lemon-project-template-glpk

comparison deps/glpk/src/glplux.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:aafb34031aeb
1 /* glplux.c */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9 * E-mail: <mao@gnu.org>.
10 *
11 * GLPK is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19 * License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23 ***********************************************************************/
24
25 #include "glplux.h"
26 #define xfault xerror
27 #define dmp_create_poolx(size) dmp_create_pool()
28
29 /*----------------------------------------------------------------------
30 // lux_create - create LU-factorization.
31 //
32 // SYNOPSIS
33 //
34 // #include "glplux.h"
35 // LUX *lux_create(int n);
36 //
37 // DESCRIPTION
38 //
39 // The routine lux_create creates LU-factorization data structure for
40 // a matrix of the order n. Initially the factorization corresponds to
41 // the unity matrix (F = V = P = Q = I, so A = I).
42 //
43 // RETURNS
44 //
45 // The routine returns a pointer to the created LU-factorization data
46 // structure, which represents the unity matrix of the order n. */
47
48 LUX *lux_create(int n)
49 { LUX *lux;
50 int k;
51 if (n < 1)
52 xfault("lux_create: n = %d; invalid parameter\n", n);
53 lux = xmalloc(sizeof(LUX));
54 lux->n = n;
55 lux->pool = dmp_create_poolx(sizeof(LUXELM));
56 lux->F_row = xcalloc(1+n, sizeof(LUXELM *));
57 lux->F_col = xcalloc(1+n, sizeof(LUXELM *));
58 lux->V_piv = xcalloc(1+n, sizeof(mpq_t));
59 lux->V_row = xcalloc(1+n, sizeof(LUXELM *));
60 lux->V_col = xcalloc(1+n, sizeof(LUXELM *));
61 lux->P_row = xcalloc(1+n, sizeof(int));
62 lux->P_col = xcalloc(1+n, sizeof(int));
63 lux->Q_row = xcalloc(1+n, sizeof(int));
64 lux->Q_col = xcalloc(1+n, sizeof(int));
65 for (k = 1; k <= n; k++)
66 { lux->F_row[k] = lux->F_col[k] = NULL;
67 mpq_init(lux->V_piv[k]);
68 mpq_set_si(lux->V_piv[k], 1, 1);
69 lux->V_row[k] = lux->V_col[k] = NULL;
70 lux->P_row[k] = lux->P_col[k] = k;
71 lux->Q_row[k] = lux->Q_col[k] = k;
72 }
73 lux->rank = n;
74 return lux;
75 }
76
77 /*----------------------------------------------------------------------
78 // initialize - initialize LU-factorization data structures.
79 //
80 // This routine initializes data structures for subsequent computing
81 // the LU-factorization of a given matrix A, which is specified by the
82 // formal routine col. On exit V = A and F = P = Q = I, where I is the
83 // unity matrix. */
84
85 static void initialize(LUX *lux, int (*col)(void *info, int j,
86 int ind[], mpq_t val[]), void *info, LUXWKA *wka)
87 { int n = lux->n;
88 DMP *pool = lux->pool;
89 LUXELM **F_row = lux->F_row;
90 LUXELM **F_col = lux->F_col;
91 mpq_t *V_piv = lux->V_piv;
92 LUXELM **V_row = lux->V_row;
93 LUXELM **V_col = lux->V_col;
94 int *P_row = lux->P_row;
95 int *P_col = lux->P_col;
96 int *Q_row = lux->Q_row;
97 int *Q_col = lux->Q_col;
98 int *R_len = wka->R_len;
99 int *R_head = wka->R_head;
100 int *R_prev = wka->R_prev;
101 int *R_next = wka->R_next;
102 int *C_len = wka->C_len;
103 int *C_head = wka->C_head;
104 int *C_prev = wka->C_prev;
105 int *C_next = wka->C_next;
106 LUXELM *fij, *vij;
107 int i, j, k, len, *ind;
108 mpq_t *val;
109 /* F := I */
110 for (i = 1; i <= n; i++)
111 { while (F_row[i] != NULL)
112 { fij = F_row[i], F_row[i] = fij->r_next;
113 mpq_clear(fij->val);
114 dmp_free_atom(pool, fij, sizeof(LUXELM));
115 }
116 }
117 for (j = 1; j <= n; j++) F_col[j] = NULL;
118 /* V := 0 */
119 for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1);
120 for (i = 1; i <= n; i++)
121 { while (V_row[i] != NULL)
122 { vij = V_row[i], V_row[i] = vij->r_next;
123 mpq_clear(vij->val);
124 dmp_free_atom(pool, vij, sizeof(LUXELM));
125 }
126 }
127 for (j = 1; j <= n; j++) V_col[j] = NULL;
128 /* V := A */
129 ind = xcalloc(1+n, sizeof(int));
130 val = xcalloc(1+n, sizeof(mpq_t));
131 for (k = 1; k <= n; k++) mpq_init(val[k]);
132 for (j = 1; j <= n; j++)
133 { /* obtain j-th column of matrix A */
134 len = col(info, j, ind, val);
135 if (!(0 <= len && len <= n))
136 xfault("lux_decomp: j = %d: len = %d; invalid column length"
137 "\n", j, len);
138 /* copy elements of j-th column to matrix V */
139 for (k = 1; k <= len; k++)
140 { /* get row index of a[i,j] */
141 i = ind[k];
142 if (!(1 <= i && i <= n))
143 xfault("lux_decomp: j = %d: i = %d; row index out of ran"
144 "ge\n", j, i);
145 /* check for duplicate indices */
146 if (V_row[i] != NULL && V_row[i]->j == j)
147 xfault("lux_decomp: j = %d: i = %d; duplicate row indice"
148 "s not allowed\n", j, i);
149 /* check for zero value */
150 if (mpq_sgn(val[k]) == 0)
151 xfault("lux_decomp: j = %d: i = %d; zero elements not al"
152 "lowed\n", j, i);
153 /* add new element v[i,j] = a[i,j] to V */
154 vij = dmp_get_atom(pool, sizeof(LUXELM));
155 vij->i = i, vij->j = j;
156 mpq_init(vij->val);
157 mpq_set(vij->val, val[k]);
158 vij->r_prev = NULL;
159 vij->r_next = V_row[i];
160 vij->c_prev = NULL;
161 vij->c_next = V_col[j];
162 if (vij->r_next != NULL) vij->r_next->r_prev = vij;
163 if (vij->c_next != NULL) vij->c_next->c_prev = vij;
164 V_row[i] = V_col[j] = vij;
165 }
166 }
167 xfree(ind);
168 for (k = 1; k <= n; k++) mpq_clear(val[k]);
169 xfree(val);
170 /* P := Q := I */
171 for (k = 1; k <= n; k++)
172 P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k;
173 /* the rank of A and V is not determined yet */
174 lux->rank = -1;
175 /* initially the entire matrix V is active */
176 /* determine its row lengths */
177 for (i = 1; i <= n; i++)
178 { len = 0;
179 for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++;
180 R_len[i] = len;
181 }
182 /* build linked lists of active rows */
183 for (len = 0; len <= n; len++) R_head[len] = 0;
184 for (i = 1; i <= n; i++)
185 { len = R_len[i];
186 R_prev[i] = 0;
187 R_next[i] = R_head[len];
188 if (R_next[i] != 0) R_prev[R_next[i]] = i;
189 R_head[len] = i;
190 }
191 /* determine its column lengths */
192 for (j = 1; j <= n; j++)
193 { len = 0;
194 for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++;
195 C_len[j] = len;
196 }
197 /* build linked lists of active columns */
198 for (len = 0; len <= n; len++) C_head[len] = 0;
199 for (j = 1; j <= n; j++)
200 { len = C_len[j];
201 C_prev[j] = 0;
202 C_next[j] = C_head[len];
203 if (C_next[j] != 0) C_prev[C_next[j]] = j;
204 C_head[len] = j;
205 }
206 return;
207 }
208
209 /*----------------------------------------------------------------------
210 // find_pivot - choose a pivot element.
211 //
212 // This routine chooses a pivot element v[p,q] in the active submatrix
213 // of matrix U = P*V*Q.
214 //
215 // It is assumed that on entry the matrix U has the following partially
216 // triangularized form:
217 //
218 // 1 k n
219 // 1 x x x x x x x x x x
220 // . x x x x x x x x x
221 // . . x x x x x x x x
222 // . . . x x x x x x x
223 // k . . . . * * * * * *
224 // . . . . * * * * * *
225 // . . . . * * * * * *
226 // . . . . * * * * * *
227 // . . . . * * * * * *
228 // n . . . . * * * * * *
229 //
230 // where rows and columns k, k+1, ..., n belong to the active submatrix
231 // (elements of the active submatrix are marked by '*').
232 //
233 // Since the matrix U = P*V*Q is not stored, the routine works with the
234 // matrix V. It is assumed that the row-wise representation corresponds
235 // to the matrix V, but the column-wise representation corresponds to
236 // the active submatrix of the matrix V, i.e. elements of the matrix V,
237 // which does not belong to the active submatrix, are missing from the
238 // column linked lists. It is also assumed that each active row of the
239 // matrix V is in the set R[len], where len is number of non-zeros in
240 // the row, and each active column of the matrix V is in the set C[len],
241 // where len is number of non-zeros in the column (in the latter case
242 // only elements of the active submatrix are counted; such elements are
243 // marked by '*' on the figure above).
244 //
245 // Due to exact arithmetic any non-zero element of the active submatrix
246 // can be chosen as a pivot. However, to keep sparsity of the matrix V
247 // the routine uses Markowitz strategy, trying to choose such element
248 // v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1),
249 // where nr[p] and nc[q] are the number of non-zero elements, resp., in
250 // p-th row and in q-th column of the active submatrix.
251 //
252 // In order to reduce the search, i.e. not to walk through all elements
253 // of the active submatrix, the routine exploits a technique proposed by
254 // I.Duff. This technique is based on using the sets R[len] and C[len]
255 // of active rows and columns.
256 //
257 // On exit the routine returns a pointer to a pivot v[p,q] chosen, or
258 // NULL, if the active submatrix is empty. */
259
260 static LUXELM *find_pivot(LUX *lux, LUXWKA *wka)
261 { int n = lux->n;
262 LUXELM **V_row = lux->V_row;
263 LUXELM **V_col = lux->V_col;
264 int *R_len = wka->R_len;
265 int *R_head = wka->R_head;
266 int *R_next = wka->R_next;
267 int *C_len = wka->C_len;
268 int *C_head = wka->C_head;
269 int *C_next = wka->C_next;
270 LUXELM *piv, *some, *vij;
271 int i, j, len, min_len, ncand, piv_lim = 5;
272 double best, cost;
273 /* nothing is chosen so far */
274 piv = NULL, best = DBL_MAX, ncand = 0;
275 /* if in the active submatrix there is a column that has the only
276 non-zero (column singleton), choose it as a pivot */
277 j = C_head[1];
278 if (j != 0)
279 { xassert(C_len[j] == 1);
280 piv = V_col[j];
281 xassert(piv != NULL && piv->c_next == NULL);
282 goto done;
283 }
284 /* if in the active submatrix there is a row that has the only
285 non-zero (row singleton), choose it as a pivot */
286 i = R_head[1];
287 if (i != 0)
288 { xassert(R_len[i] == 1);
289 piv = V_row[i];
290 xassert(piv != NULL && piv->r_next == NULL);
291 goto done;
292 }
293 /* there are no singletons in the active submatrix; walk through
294 other non-empty rows and columns */
295 for (len = 2; len <= n; len++)
296 { /* consider active columns having len non-zeros */
297 for (j = C_head[len]; j != 0; j = C_next[j])
298 { /* j-th column has len non-zeros */
299 /* find an element in the row of minimal length */
300 some = NULL, min_len = INT_MAX;
301 for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
302 { if (min_len > R_len[vij->i])
303 some = vij, min_len = R_len[vij->i];
304 /* if Markowitz cost of this element is not greater than
305 (len-1)**2, it can be chosen right now; this heuristic
306 reduces the search and works well in many cases */
307 if (min_len <= len)
308 { piv = some;
309 goto done;
310 }
311 }
312 /* j-th column has been scanned */
313 /* the minimal element found is a next pivot candidate */
314 xassert(some != NULL);
315 ncand++;
316 /* compute its Markowitz cost */
317 cost = (double)(min_len - 1) * (double)(len - 1);
318 /* choose between the current candidate and this element */
319 if (cost < best) piv = some, best = cost;
320 /* if piv_lim candidates have been considered, there is a
321 doubt that a much better candidate exists; therefore it
322 is the time to terminate the search */
323 if (ncand == piv_lim) goto done;
324 }
325 /* now consider active rows having len non-zeros */
326 for (i = R_head[len]; i != 0; i = R_next[i])
327 { /* i-th row has len non-zeros */
328 /* find an element in the column of minimal length */
329 some = NULL, min_len = INT_MAX;
330 for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
331 { if (min_len > C_len[vij->j])
332 some = vij, min_len = C_len[vij->j];
333 /* if Markowitz cost of this element is not greater than
334 (len-1)**2, it can be chosen right now; this heuristic
335 reduces the search and works well in many cases */
336 if (min_len <= len)
337 { piv = some;
338 goto done;
339 }
340 }
341 /* i-th row has been scanned */
342 /* the minimal element found is a next pivot candidate */
343 xassert(some != NULL);
344 ncand++;
345 /* compute its Markowitz cost */
346 cost = (double)(len - 1) * (double)(min_len - 1);
347 /* choose between the current candidate and this element */
348 if (cost < best) piv = some, best = cost;
349 /* if piv_lim candidates have been considered, there is a
350 doubt that a much better candidate exists; therefore it
351 is the time to terminate the search */
352 if (ncand == piv_lim) goto done;
353 }
354 }
355 done: /* bring the pivot v[p,q] to the factorizing routine */
356 return piv;
357 }
358
359 /*----------------------------------------------------------------------
360 // eliminate - perform gaussian elimination.
361 //
362 // This routine performs elementary gaussian transformations in order
363 // to eliminate subdiagonal elements in the k-th column of the matrix
364 // U = P*V*Q using the pivot element u[k,k], where k is the number of
365 // the current elimination step.
366 //
367 // The parameter piv specifies the pivot element v[p,q] = u[k,k].
368 //
369 // Each time when the routine applies the elementary transformation to
370 // a non-pivot row of the matrix V, it stores the corresponding element
371 // to the matrix F in order to keep the main equality A = F*V.
372 //
373 // The routine assumes that on entry the matrices L = P*F*inv(P) and
374 // U = P*V*Q are the following:
375 //
376 // 1 k 1 k n
377 // 1 1 . . . . . . . . . 1 x x x x x x x x x x
378 // x 1 . . . . . . . . . x x x x x x x x x
379 // x x 1 . . . . . . . . . x x x x x x x x
380 // x x x 1 . . . . . . . . . x x x x x x x
381 // k x x x x 1 . . . . . k . . . . * * * * * *
382 // x x x x _ 1 . . . . . . . . # * * * * *
383 // x x x x _ . 1 . . . . . . . # * * * * *
384 // x x x x _ . . 1 . . . . . . # * * * * *
385 // x x x x _ . . . 1 . . . . . # * * * * *
386 // n x x x x _ . . . . 1 n . . . . # * * * * *
387 //
388 // matrix L matrix U
389 //
390 // where rows and columns of the matrix U with numbers k, k+1, ..., n
391 // form the active submatrix (eliminated elements are marked by '#' and
392 // other elements of the active submatrix are marked by '*'). Note that
393 // each eliminated non-zero element u[i,k] of the matrix U gives the
394 // corresponding element l[i,k] of the matrix L (marked by '_').
395 //
396 // Actually all operations are performed on the matrix V. Should note
397 // that the row-wise representation corresponds to the matrix V, but the
398 // column-wise representation corresponds to the active submatrix of the
399 // matrix V, i.e. elements of the matrix V, which doesn't belong to the
400 // active submatrix, are missing from the column linked lists.
401 //
402 // Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal
403 // elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies
404 // the following elementary gaussian transformations:
405 //
406 // (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V),
407 //
408 // where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier.
409 //
410 // Additionally, in order to keep the main equality A = F*V, each time
411 // when the routine applies the transformation to i-th row of the matrix
412 // V, it also adds f[i,p] as a new element to the matrix F.
413 //
414 // IMPORTANT: On entry the working arrays flag and work should contain
415 // zeros. This status is provided by the routine on exit. */
416
417 static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[],
418 mpq_t work[])
419 { DMP *pool = lux->pool;
420 LUXELM **F_row = lux->F_row;
421 LUXELM **F_col = lux->F_col;
422 mpq_t *V_piv = lux->V_piv;
423 LUXELM **V_row = lux->V_row;
424 LUXELM **V_col = lux->V_col;
425 int *R_len = wka->R_len;
426 int *R_head = wka->R_head;
427 int *R_prev = wka->R_prev;
428 int *R_next = wka->R_next;
429 int *C_len = wka->C_len;
430 int *C_head = wka->C_head;
431 int *C_prev = wka->C_prev;
432 int *C_next = wka->C_next;
433 LUXELM *fip, *vij, *vpj, *viq, *next;
434 mpq_t temp;
435 int i, j, p, q;
436 mpq_init(temp);
437 /* determine row and column indices of the pivot v[p,q] */
438 xassert(piv != NULL);
439 p = piv->i, q = piv->j;
440 /* remove p-th (pivot) row from the active set; it will never
441 return there */
442 if (R_prev[p] == 0)
443 R_head[R_len[p]] = R_next[p];
444 else
445 R_next[R_prev[p]] = R_next[p];
446 if (R_next[p] == 0)
447 ;
448 else
449 R_prev[R_next[p]] = R_prev[p];
450 /* remove q-th (pivot) column from the active set; it will never
451 return there */
452 if (C_prev[q] == 0)
453 C_head[C_len[q]] = C_next[q];
454 else
455 C_next[C_prev[q]] = C_next[q];
456 if (C_next[q] == 0)
457 ;
458 else
459 C_prev[C_next[q]] = C_prev[q];
460 /* store the pivot value in a separate array */
461 mpq_set(V_piv[p], piv->val);
462 /* remove the pivot from p-th row */
463 if (piv->r_prev == NULL)
464 V_row[p] = piv->r_next;
465 else
466 piv->r_prev->r_next = piv->r_next;
467 if (piv->r_next == NULL)
468 ;
469 else
470 piv->r_next->r_prev = piv->r_prev;
471 R_len[p]--;
472 /* remove the pivot from q-th column */
473 if (piv->c_prev == NULL)
474 V_col[q] = piv->c_next;
475 else
476 piv->c_prev->c_next = piv->c_next;
477 if (piv->c_next == NULL)
478 ;
479 else
480 piv->c_next->c_prev = piv->c_prev;
481 C_len[q]--;
482 /* free the space occupied by the pivot */
483 mpq_clear(piv->val);
484 dmp_free_atom(pool, piv, sizeof(LUXELM));
485 /* walk through p-th (pivot) row, which already does not contain
486 the pivot v[p,q], and do the following... */
487 for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
488 { /* get column index of v[p,j] */
489 j = vpj->j;
490 /* store v[p,j] in the working array */
491 flag[j] = 1;
492 mpq_set(work[j], vpj->val);
493 /* remove j-th column from the active set; it will return there
494 later with a new length */
495 if (C_prev[j] == 0)
496 C_head[C_len[j]] = C_next[j];
497 else
498 C_next[C_prev[j]] = C_next[j];
499 if (C_next[j] == 0)
500 ;
501 else
502 C_prev[C_next[j]] = C_prev[j];
503 /* v[p,j] leaves the active submatrix, so remove it from j-th
504 column; however, v[p,j] is kept in p-th row */
505 if (vpj->c_prev == NULL)
506 V_col[j] = vpj->c_next;
507 else
508 vpj->c_prev->c_next = vpj->c_next;
509 if (vpj->c_next == NULL)
510 ;
511 else
512 vpj->c_next->c_prev = vpj->c_prev;
513 C_len[j]--;
514 }
515 /* now walk through q-th (pivot) column, which already does not
516 contain the pivot v[p,q], and perform gaussian elimination */
517 while (V_col[q] != NULL)
518 { /* element v[i,q] has to be eliminated */
519 viq = V_col[q];
520 /* get row index of v[i,q] */
521 i = viq->i;
522 /* remove i-th row from the active set; later it will return
523 there with a new length */
524 if (R_prev[i] == 0)
525 R_head[R_len[i]] = R_next[i];
526 else
527 R_next[R_prev[i]] = R_next[i];
528 if (R_next[i] == 0)
529 ;
530 else
531 R_prev[R_next[i]] = R_prev[i];
532 /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and
533 store it in the matrix F */
534 fip = dmp_get_atom(pool, sizeof(LUXELM));
535 fip->i = i, fip->j = p;
536 mpq_init(fip->val);
537 mpq_div(fip->val, viq->val, V_piv[p]);
538 fip->r_prev = NULL;
539 fip->r_next = F_row[i];
540 fip->c_prev = NULL;
541 fip->c_next = F_col[p];
542 if (fip->r_next != NULL) fip->r_next->r_prev = fip;
543 if (fip->c_next != NULL) fip->c_next->c_prev = fip;
544 F_row[i] = F_col[p] = fip;
545 /* v[i,q] has to be eliminated, so remove it from i-th row */
546 if (viq->r_prev == NULL)
547 V_row[i] = viq->r_next;
548 else
549 viq->r_prev->r_next = viq->r_next;
550 if (viq->r_next == NULL)
551 ;
552 else
553 viq->r_next->r_prev = viq->r_prev;
554 R_len[i]--;
555 /* and also from q-th column */
556 V_col[q] = viq->c_next;
557 C_len[q]--;
558 /* free the space occupied by v[i,q] */
559 mpq_clear(viq->val);
560 dmp_free_atom(pool, viq, sizeof(LUXELM));
561 /* perform gaussian transformation:
562 (i-th row) := (i-th row) - f[i,p] * (p-th row)
563 note that now p-th row, which is in the working array,
564 does not contain the pivot v[p,q], and i-th row does not
565 contain the element v[i,q] to be eliminated */
566 /* walk through i-th row and transform existing non-zero
567 elements */
568 for (vij = V_row[i]; vij != NULL; vij = next)
569 { next = vij->r_next;
570 /* get column index of v[i,j] */
571 j = vij->j;
572 /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
573 if (flag[j])
574 { /* v[p,j] != 0 */
575 flag[j] = 0;
576 mpq_mul(temp, fip->val, work[j]);
577 mpq_sub(vij->val, vij->val, temp);
578 if (mpq_sgn(vij->val) == 0)
579 { /* new v[i,j] is zero, so remove it from the active
580 submatrix */
581 /* remove v[i,j] from i-th row */
582 if (vij->r_prev == NULL)
583 V_row[i] = vij->r_next;
584 else
585 vij->r_prev->r_next = vij->r_next;
586 if (vij->r_next == NULL)
587 ;
588 else
589 vij->r_next->r_prev = vij->r_prev;
590 R_len[i]--;
591 /* remove v[i,j] from j-th column */
592 if (vij->c_prev == NULL)
593 V_col[j] = vij->c_next;
594 else
595 vij->c_prev->c_next = vij->c_next;
596 if (vij->c_next == NULL)
597 ;
598 else
599 vij->c_next->c_prev = vij->c_prev;
600 C_len[j]--;
601 /* free the space occupied by v[i,j] */
602 mpq_clear(vij->val);
603 dmp_free_atom(pool, vij, sizeof(LUXELM));
604 }
605 }
606 }
607 /* now flag is the pattern of the set v[p,*] \ v[i,*] */
608 /* walk through p-th (pivot) row and create new elements in
609 i-th row, which appear due to fill-in */
610 for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
611 { j = vpj->j;
612 if (flag[j])
613 { /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and
614 add it to i-th row and j-th column */
615 vij = dmp_get_atom(pool, sizeof(LUXELM));
616 vij->i = i, vij->j = j;
617 mpq_init(vij->val);
618 mpq_mul(vij->val, fip->val, work[j]);
619 mpq_neg(vij->val, vij->val);
620 vij->r_prev = NULL;
621 vij->r_next = V_row[i];
622 vij->c_prev = NULL;
623 vij->c_next = V_col[j];
624 if (vij->r_next != NULL) vij->r_next->r_prev = vij;
625 if (vij->c_next != NULL) vij->c_next->c_prev = vij;
626 V_row[i] = V_col[j] = vij;
627 R_len[i]++, C_len[j]++;
628 }
629 else
630 { /* there is no fill-in, because v[i,j] already exists in
631 i-th row; restore the flag, which was reset before */
632 flag[j] = 1;
633 }
634 }
635 /* now i-th row has been completely transformed and can return
636 to the active set with a new length */
637 R_prev[i] = 0;
638 R_next[i] = R_head[R_len[i]];
639 if (R_next[i] != 0) R_prev[R_next[i]] = i;
640 R_head[R_len[i]] = i;
641 }
642 /* at this point q-th (pivot) column must be empty */
643 xassert(C_len[q] == 0);
644 /* walk through p-th (pivot) row again and do the following... */
645 for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
646 { /* get column index of v[p,j] */
647 j = vpj->j;
648 /* erase v[p,j] from the working array */
649 flag[j] = 0;
650 mpq_set_si(work[j], 0, 1);
651 /* now j-th column has been completely transformed, so it can
652 return to the active list with a new length */
653 C_prev[j] = 0;
654 C_next[j] = C_head[C_len[j]];
655 if (C_next[j] != 0) C_prev[C_next[j]] = j;
656 C_head[C_len[j]] = j;
657 }
658 mpq_clear(temp);
659 /* return to the factorizing routine */
660 return;
661 }
662
663 /*----------------------------------------------------------------------
664 // lux_decomp - compute LU-factorization.
665 //
666 // SYNOPSIS
667 //
668 // #include "glplux.h"
669 // int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
670 // mpq_t val[]), void *info);
671 //
672 // DESCRIPTION
673 //
674 // The routine lux_decomp computes LU-factorization of a given square
675 // matrix A.
676 //
677 // The parameter lux specifies LU-factorization data structure built by
678 // means of the routine lux_create.
679 //
680 // The formal routine col specifies the original matrix A. In order to
681 // obtain j-th column of the matrix A the routine lux_decomp calls the
682 // routine col with the parameter j (1 <= j <= n, where n is the order
683 // of A). In response the routine col should store row indices and
684 // numerical values of non-zero elements of j-th column of A to the
685 // locations ind[1], ..., ind[len] and val[1], ..., val[len], resp.,
686 // where len is the number of non-zeros in j-th column, which should be
687 // returned on exit. Neiter zero nor duplicate elements are allowed.
688 //
689 // The parameter info is a transit pointer passed to the formal routine
690 // col; it can be used for various purposes.
691 //
692 // RETURNS
693 //
694 // The routine lux_decomp returns the singularity flag. Zero flag means
695 // that the original matrix A is non-singular while non-zero flag means
696 // that A is (exactly!) singular.
697 //
698 // Note that LU-factorization is valid in both cases, however, in case
699 // of singularity some rows of the matrix V (including pivot elements)
700 // will be empty.
701 //
702 // REPAIRING SINGULAR MATRIX
703 //
704 // If the routine lux_decomp returns non-zero flag, it provides all
705 // necessary information that can be used for "repairing" the matrix A,
706 // where "repairing" means replacing linearly dependent columns of the
707 // matrix A by appropriate columns of the unity matrix. This feature is
708 // needed when the routine lux_decomp is used for reinverting the basis
709 // matrix within the simplex method procedure.
710 //
711 // On exit linearly dependent columns of the matrix U have the numbers
712 // rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A
713 // stored by the routine to the member lux->rank. The correspondence
714 // between columns of A and U is the same as between columns of V and U.
715 // Thus, linearly dependent columns of the matrix A have the numbers
716 // Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array
717 // representing the permutation matrix Q in column-like format. It is
718 // understood that each j-th linearly dependent column of the matrix U
719 // should be replaced by the unity vector, where all elements are zero
720 // except the unity diagonal element u[j,j]. On the other hand j-th row
721 // of the matrix U corresponds to the row of the matrix V (and therefore
722 // of the matrix A) with the number P_row[j], where P_row is an array
723 // representing the permutation matrix P in row-like format. Thus, each
724 // j-th linearly dependent column of the matrix U should be replaced by
725 // a column of the unity matrix with the number P_row[j].
726 //
727 // The code that repairs the matrix A may look like follows:
728 //
729 // for (j = rank+1; j <= n; j++)
730 // { replace column Q_col[j] of the matrix A by column P_row[j] of
731 // the unity matrix;
732 // }
733 //
734 // where rank, P_row, and Q_col are members of the structure LUX. */
735
736 int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
737 mpq_t val[]), void *info)
738 { int n = lux->n;
739 LUXELM **V_row = lux->V_row;
740 LUXELM **V_col = lux->V_col;
741 int *P_row = lux->P_row;
742 int *P_col = lux->P_col;
743 int *Q_row = lux->Q_row;
744 int *Q_col = lux->Q_col;
745 LUXELM *piv, *vij;
746 LUXWKA *wka;
747 int i, j, k, p, q, t, *flag;
748 mpq_t *work;
749 /* allocate working area */
750 wka = xmalloc(sizeof(LUXWKA));
751 wka->R_len = xcalloc(1+n, sizeof(int));
752 wka->R_head = xcalloc(1+n, sizeof(int));
753 wka->R_prev = xcalloc(1+n, sizeof(int));
754 wka->R_next = xcalloc(1+n, sizeof(int));
755 wka->C_len = xcalloc(1+n, sizeof(int));
756 wka->C_head = xcalloc(1+n, sizeof(int));
757 wka->C_prev = xcalloc(1+n, sizeof(int));
758 wka->C_next = xcalloc(1+n, sizeof(int));
759 /* initialize LU-factorization data structures */
760 initialize(lux, col, info, wka);
761 /* allocate working arrays */
762 flag = xcalloc(1+n, sizeof(int));
763 work = xcalloc(1+n, sizeof(mpq_t));
764 for (k = 1; k <= n; k++)
765 { flag[k] = 0;
766 mpq_init(work[k]);
767 }
768 /* main elimination loop */
769 for (k = 1; k <= n; k++)
770 { /* choose a pivot element v[p,q] */
771 piv = find_pivot(lux, wka);
772 if (piv == NULL)
773 { /* no pivot can be chosen, because the active submatrix is
774 empty */
775 break;
776 }
777 /* determine row and column indices of the pivot element */
778 p = piv->i, q = piv->j;
779 /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th
780 rows and k-th and j'-th columns of the matrix U = P*V*Q to
781 move the element u[i',j'] to the position u[k,k] */
782 i = P_col[p], j = Q_row[q];
783 xassert(k <= i && i <= n && k <= j && j <= n);
784 /* permute k-th and i-th rows of the matrix U */
785 t = P_row[k];
786 P_row[i] = t, P_col[t] = i;
787 P_row[k] = p, P_col[p] = k;
788 /* permute k-th and j-th columns of the matrix U */
789 t = Q_col[k];
790 Q_col[j] = t, Q_row[t] = j;
791 Q_col[k] = q, Q_row[q] = k;
792 /* eliminate subdiagonal elements of k-th column of the matrix
793 U = P*V*Q using the pivot element u[k,k] = v[p,q] */
794 eliminate(lux, wka, piv, flag, work);
795 }
796 /* determine the rank of A (and V) */
797 lux->rank = k - 1;
798 /* free working arrays */
799 xfree(flag);
800 for (k = 1; k <= n; k++) mpq_clear(work[k]);
801 xfree(work);
802 /* build column lists of the matrix V using its row lists */
803 for (j = 1; j <= n; j++)
804 xassert(V_col[j] == NULL);
805 for (i = 1; i <= n; i++)
806 { for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
807 { j = vij->j;
808 vij->c_prev = NULL;
809 vij->c_next = V_col[j];
810 if (vij->c_next != NULL) vij->c_next->c_prev = vij;
811 V_col[j] = vij;
812 }
813 }
814 /* free working area */
815 xfree(wka->R_len);
816 xfree(wka->R_head);
817 xfree(wka->R_prev);
818 xfree(wka->R_next);
819 xfree(wka->C_len);
820 xfree(wka->C_head);
821 xfree(wka->C_prev);
822 xfree(wka->C_next);
823 xfree(wka);
824 /* return to the calling program */
825 return (lux->rank < n);
826 }
827
828 /*----------------------------------------------------------------------
829 // lux_f_solve - solve system F*x = b or F'*x = b.
830 //
831 // SYNOPSIS
832 //
833 // #include "glplux.h"
834 // void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
835 //
836 // DESCRIPTION
837 //
838 // The routine lux_f_solve solves either the system F*x = b (if the
839 // flag tr is zero) or the system F'*x = b (if the flag tr is non-zero),
840 // where the matrix F is a component of LU-factorization specified by
841 // the parameter lux, F' is a matrix transposed to F.
842 //
843 // On entry the array x should contain elements of the right-hand side
844 // vector b in locations x[1], ..., x[n], where n is the order of the
845 // matrix F. On exit this array will contain elements of the solution
846 // vector x in the same locations. */
847
848 void lux_f_solve(LUX *lux, int tr, mpq_t x[])
849 { int n = lux->n;
850 LUXELM **F_row = lux->F_row;
851 LUXELM **F_col = lux->F_col;
852 int *P_row = lux->P_row;
853 LUXELM *fik, *fkj;
854 int i, j, k;
855 mpq_t temp;
856 mpq_init(temp);
857 if (!tr)
858 { /* solve the system F*x = b */
859 for (j = 1; j <= n; j++)
860 { k = P_row[j];
861 if (mpq_sgn(x[k]) != 0)
862 { for (fik = F_col[k]; fik != NULL; fik = fik->c_next)
863 { mpq_mul(temp, fik->val, x[k]);
864 mpq_sub(x[fik->i], x[fik->i], temp);
865 }
866 }
867 }
868 }
869 else
870 { /* solve the system F'*x = b */
871 for (i = n; i >= 1; i--)
872 { k = P_row[i];
873 if (mpq_sgn(x[k]) != 0)
874 { for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next)
875 { mpq_mul(temp, fkj->val, x[k]);
876 mpq_sub(x[fkj->j], x[fkj->j], temp);
877 }
878 }
879 }
880 }
881 mpq_clear(temp);
882 return;
883 }
884
885 /*----------------------------------------------------------------------
886 // lux_v_solve - solve system V*x = b or V'*x = b.
887 //
888 // SYNOPSIS
889 //
890 // #include "glplux.h"
891 // void lux_v_solve(LUX *lux, int tr, double x[]);
892 //
893 // DESCRIPTION
894 //
895 // The routine lux_v_solve solves either the system V*x = b (if the
896 // flag tr is zero) or the system V'*x = b (if the flag tr is non-zero),
897 // where the matrix V is a component of LU-factorization specified by
898 // the parameter lux, V' is a matrix transposed to V.
899 //
900 // On entry the array x should contain elements of the right-hand side
901 // vector b in locations x[1], ..., x[n], where n is the order of the
902 // matrix V. On exit this array will contain elements of the solution
903 // vector x in the same locations. */
904
905 void lux_v_solve(LUX *lux, int tr, mpq_t x[])
906 { int n = lux->n;
907 mpq_t *V_piv = lux->V_piv;
908 LUXELM **V_row = lux->V_row;
909 LUXELM **V_col = lux->V_col;
910 int *P_row = lux->P_row;
911 int *Q_col = lux->Q_col;
912 LUXELM *vij;
913 int i, j, k;
914 mpq_t *b, temp;
915 b = xcalloc(1+n, sizeof(mpq_t));
916 for (k = 1; k <= n; k++)
917 mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1);
918 mpq_init(temp);
919 if (!tr)
920 { /* solve the system V*x = b */
921 for (k = n; k >= 1; k--)
922 { i = P_row[k], j = Q_col[k];
923 if (mpq_sgn(b[i]) != 0)
924 { mpq_set(x[j], b[i]);
925 mpq_div(x[j], x[j], V_piv[i]);
926 for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
927 { mpq_mul(temp, vij->val, x[j]);
928 mpq_sub(b[vij->i], b[vij->i], temp);
929 }
930 }
931 }
932 }
933 else
934 { /* solve the system V'*x = b */
935 for (k = 1; k <= n; k++)
936 { i = P_row[k], j = Q_col[k];
937 if (mpq_sgn(b[j]) != 0)
938 { mpq_set(x[i], b[j]);
939 mpq_div(x[i], x[i], V_piv[i]);
940 for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
941 { mpq_mul(temp, vij->val, x[i]);
942 mpq_sub(b[vij->j], b[vij->j], temp);
943 }
944 }
945 }
946 }
947 for (k = 1; k <= n; k++) mpq_clear(b[k]);
948 mpq_clear(temp);
949 xfree(b);
950 return;
951 }
952
953 /*----------------------------------------------------------------------
954 // lux_solve - solve system A*x = b or A'*x = b.
955 //
956 // SYNOPSIS
957 //
958 // #include "glplux.h"
959 // void lux_solve(LUX *lux, int tr, mpq_t x[]);
960 //
961 // DESCRIPTION
962 //
963 // The routine lux_solve solves either the system A*x = b (if the flag
964 // tr is zero) or the system A'*x = b (if the flag tr is non-zero),
965 // where the parameter lux specifies LU-factorization of the matrix A,
966 // A' is a matrix transposed to A.
967 //
968 // On entry the array x should contain elements of the right-hand side
969 // vector b in locations x[1], ..., x[n], where n is the order of the
970 // matrix A. On exit this array will contain elements of the solution
971 // vector x in the same locations. */
972
973 void lux_solve(LUX *lux, int tr, mpq_t x[])
974 { if (lux->rank < lux->n)
975 xfault("lux_solve: LU-factorization has incomplete rank\n");
976 if (!tr)
977 { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
978 lux_f_solve(lux, 0, x);
979 lux_v_solve(lux, 0, x);
980 }
981 else
982 { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
983 lux_v_solve(lux, 1, x);
984 lux_f_solve(lux, 1, x);
985 }
986 return;
987 }
988
989 /*----------------------------------------------------------------------
990 // lux_delete - delete LU-factorization.
991 //
992 // SYNOPSIS
993 //
994 // #include "glplux.h"
995 // void lux_delete(LUX *lux);
996 //
997 // DESCRIPTION
998 //
999 // The routine lux_delete deletes LU-factorization data structure,
1000 // which the parameter lux points to, freeing all the memory allocated
1001 // to this object. */
1002
1003 void lux_delete(LUX *lux)
1004 { int n = lux->n;
1005 LUXELM *fij, *vij;
1006 int i;
1007 for (i = 1; i <= n; i++)
1008 { for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next)
1009 mpq_clear(fij->val);
1010 mpq_clear(lux->V_piv[i]);
1011 for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next)
1012 mpq_clear(vij->val);
1013 }
1014 dmp_delete_pool(lux->pool);
1015 xfree(lux->F_row);
1016 xfree(lux->F_col);
1017 xfree(lux->V_piv);
1018 xfree(lux->V_row);
1019 xfree(lux->V_col);
1020 xfree(lux->P_row);
1021 xfree(lux->P_col);
1022 xfree(lux->Q_row);
1023 xfree(lux->Q_col);
1024 xfree(lux);
1025 return;
1026 }
1027
1028 /* eof */