COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glplux.c @ 10:5545663ca997

subpack-glpk
Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago

Import GLPK 4.47

File size: 38.0 KB
Line 
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
48LUX *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
85static 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
260static 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      }
355done: /* 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
417static 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
736int 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
848void 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
905void 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
973void 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
1003void 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 */
Note: See TracBrowser for help on using the repository browser.