COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpluf.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 69.9 KB
Line 
1/* glpluf.c (LU-factorization) */
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 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 "glpenv.h"
26#include "glpluf.h"
27#define xfault xerror
28
29/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
30
31#define N_MAX 100000000 /* = 100*10^6 */
32/* maximal order of the original matrix */
33
34/***********************************************************************
35*  NAME
36*
37*  luf_create_it - create LU-factorization
38*
39*  SYNOPSIS
40*
41*  #include "glpluf.h"
42*  LUF *luf_create_it(void);
43*
44*  DESCRIPTION
45*
46*  The routine luf_create_it creates a program object, which represents
47*  LU-factorization of a square matrix.
48*
49*  RETURNS
50*
51*  The routine luf_create_it returns a pointer to the object created. */
52
53LUF *luf_create_it(void)
54{     LUF *luf;
55      luf = xmalloc(sizeof(LUF));
56      luf->n_max = luf->n = 0;
57      luf->valid = 0;
58      luf->fr_ptr = luf->fr_len = NULL;
59      luf->fc_ptr = luf->fc_len = NULL;
60      luf->vr_ptr = luf->vr_len = luf->vr_cap = NULL;
61      luf->vr_piv = NULL;
62      luf->vc_ptr = luf->vc_len = luf->vc_cap = NULL;
63      luf->pp_row = luf->pp_col = NULL;
64      luf->qq_row = luf->qq_col = NULL;
65      luf->sv_size = 0;
66      luf->sv_beg = luf->sv_end = 0;
67      luf->sv_ind = NULL;
68      luf->sv_val = NULL;
69      luf->sv_head = luf->sv_tail = 0;
70      luf->sv_prev = luf->sv_next = NULL;
71      luf->vr_max = NULL;
72      luf->rs_head = luf->rs_prev = luf->rs_next = NULL;
73      luf->cs_head = luf->cs_prev = luf->cs_next = NULL;
74      luf->flag = NULL;
75      luf->work = NULL;
76      luf->new_sva = 0;
77      luf->piv_tol = 0.10;
78      luf->piv_lim = 4;
79      luf->suhl = 1;
80      luf->eps_tol = 1e-15;
81      luf->max_gro = 1e+10;
82      luf->nnz_a = luf->nnz_f = luf->nnz_v = 0;
83      luf->max_a = luf->big_v = 0.0;
84      luf->rank = 0;
85      return luf;
86}
87
88/***********************************************************************
89*  NAME
90*
91*  luf_defrag_sva - defragment the sparse vector area
92*
93*  SYNOPSIS
94*
95*  #include "glpluf.h"
96*  void luf_defrag_sva(LUF *luf);
97*
98*  DESCRIPTION
99*
100*  The routine luf_defrag_sva defragments the sparse vector area (SVA)
101*  gathering all unused locations in one continuous extent. In order to
102*  do that the routine moves all unused locations from the left part of
103*  SVA (which contains rows and columns of the matrix V) to the middle
104*  part (which contains free locations). This is attained by relocating
105*  elements of rows and columns of the matrix V toward the beginning of
106*  the left part.
107*
108*  NOTE that this "garbage collection" involves changing row and column
109*  pointers of the matrix V. */
110
111void luf_defrag_sva(LUF *luf)
112{     int n = luf->n;
113      int *vr_ptr = luf->vr_ptr;
114      int *vr_len = luf->vr_len;
115      int *vr_cap = luf->vr_cap;
116      int *vc_ptr = luf->vc_ptr;
117      int *vc_len = luf->vc_len;
118      int *vc_cap = luf->vc_cap;
119      int *sv_ind = luf->sv_ind;
120      double *sv_val = luf->sv_val;
121      int *sv_next = luf->sv_next;
122      int sv_beg = 1;
123      int i, j, k;
124      /* skip rows and columns, which do not need to be relocated */
125      for (k = luf->sv_head; k != 0; k = sv_next[k])
126      {  if (k <= n)
127         {  /* i-th row of the matrix V */
128            i = k;
129            if (vr_ptr[i] != sv_beg) break;
130            vr_cap[i] = vr_len[i];
131            sv_beg += vr_cap[i];
132         }
133         else
134         {  /* j-th column of the matrix V */
135            j = k - n;
136            if (vc_ptr[j] != sv_beg) break;
137            vc_cap[j] = vc_len[j];
138            sv_beg += vc_cap[j];
139         }
140      }
141      /* relocate other rows and columns in order to gather all unused
142         locations in one continuous extent */
143      for (k = k; k != 0; k = sv_next[k])
144      {  if (k <= n)
145         {  /* i-th row of the matrix V */
146            i = k;
147            memmove(&sv_ind[sv_beg], &sv_ind[vr_ptr[i]],
148               vr_len[i] * sizeof(int));
149            memmove(&sv_val[sv_beg], &sv_val[vr_ptr[i]],
150               vr_len[i] * sizeof(double));
151            vr_ptr[i] = sv_beg;
152            vr_cap[i] = vr_len[i];
153            sv_beg += vr_cap[i];
154         }
155         else
156         {  /* j-th column of the matrix V */
157            j = k - n;
158            memmove(&sv_ind[sv_beg], &sv_ind[vc_ptr[j]],
159               vc_len[j] * sizeof(int));
160            memmove(&sv_val[sv_beg], &sv_val[vc_ptr[j]],
161               vc_len[j] * sizeof(double));
162            vc_ptr[j] = sv_beg;
163            vc_cap[j] = vc_len[j];
164            sv_beg += vc_cap[j];
165         }
166      }
167      /* set new pointer to the beginning of the free part */
168      luf->sv_beg = sv_beg;
169      return;
170}
171
172/***********************************************************************
173*  NAME
174*
175*  luf_enlarge_row - enlarge row capacity
176*
177*  SYNOPSIS
178*
179*  #include "glpluf.h"
180*  int luf_enlarge_row(LUF *luf, int i, int cap);
181*
182*  DESCRIPTION
183*
184*  The routine luf_enlarge_row enlarges capacity of the i-th row of the
185*  matrix V to cap locations (assuming that its current capacity is less
186*  than cap). In order to do that the routine relocates elements of the
187*  i-th row to the end of the left part of SVA (which contains rows and
188*  columns of the matrix V) and then expands the left part by allocating
189*  cap free locations from the free part. If there are less than cap
190*  free locations, the routine defragments the sparse vector area.
191*
192*  Due to "garbage collection" this operation may change row and column
193*  pointers of the matrix V.
194*
195*  RETURNS
196*
197*  If no error occured, the routine returns zero. Otherwise, in case of
198*  overflow of the sparse vector area, the routine returns non-zero. */
199
200int luf_enlarge_row(LUF *luf, int i, int cap)
201{     int n = luf->n;
202      int *vr_ptr = luf->vr_ptr;
203      int *vr_len = luf->vr_len;
204      int *vr_cap = luf->vr_cap;
205      int *vc_cap = luf->vc_cap;
206      int *sv_ind = luf->sv_ind;
207      double *sv_val = luf->sv_val;
208      int *sv_prev = luf->sv_prev;
209      int *sv_next = luf->sv_next;
210      int ret = 0;
211      int cur, k, kk;
212      xassert(1 <= i && i <= n);
213      xassert(vr_cap[i] < cap);
214      /* if there are less than cap free locations, defragment SVA */
215      if (luf->sv_end - luf->sv_beg < cap)
216      {  luf_defrag_sva(luf);
217         if (luf->sv_end - luf->sv_beg < cap)
218         {  ret = 1;
219            goto done;
220         }
221      }
222      /* save current capacity of the i-th row */
223      cur = vr_cap[i];
224      /* copy existing elements to the beginning of the free part */
225      memmove(&sv_ind[luf->sv_beg], &sv_ind[vr_ptr[i]],
226         vr_len[i] * sizeof(int));
227      memmove(&sv_val[luf->sv_beg], &sv_val[vr_ptr[i]],
228         vr_len[i] * sizeof(double));
229      /* set new pointer and new capacity of the i-th row */
230      vr_ptr[i] = luf->sv_beg;
231      vr_cap[i] = cap;
232      /* set new pointer to the beginning of the free part */
233      luf->sv_beg += cap;
234      /* now the i-th row starts in the rightmost location among other
235         rows and columns of the matrix V, so its node should be moved
236         to the end of the row/column linked list */
237      k = i;
238      /* remove the i-th row node from the linked list */
239      if (sv_prev[k] == 0)
240         luf->sv_head = sv_next[k];
241      else
242      {  /* capacity of the previous row/column can be increased at the
243            expense of old locations of the i-th row */
244         kk = sv_prev[k];
245         if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
246         sv_next[sv_prev[k]] = sv_next[k];
247      }
248      if (sv_next[k] == 0)
249         luf->sv_tail = sv_prev[k];
250      else
251         sv_prev[sv_next[k]] = sv_prev[k];
252      /* insert the i-th row node to the end of the linked list */
253      sv_prev[k] = luf->sv_tail;
254      sv_next[k] = 0;
255      if (sv_prev[k] == 0)
256         luf->sv_head = k;
257      else
258         sv_next[sv_prev[k]] = k;
259      luf->sv_tail = k;
260done: return ret;
261}
262
263/***********************************************************************
264*  NAME
265*
266*  luf_enlarge_col - enlarge column capacity
267*
268*  SYNOPSIS
269*
270*  #include "glpluf.h"
271*  int luf_enlarge_col(LUF *luf, int j, int cap);
272*
273*  DESCRIPTION
274*
275*  The routine luf_enlarge_col enlarges capacity of the j-th column of
276*  the matrix V to cap locations (assuming that its current capacity is
277*  less than cap). In order to do that the routine relocates elements
278*  of the j-th column to the end of the left part of SVA (which contains
279*  rows and columns of the matrix V) and then expands the left part by
280*  allocating cap free locations from the free part. If there are less
281*  than cap free locations, the routine defragments the sparse vector
282*  area.
283*
284*  Due to "garbage collection" this operation may change row and column
285*  pointers of the matrix V.
286*
287*  RETURNS
288*
289*  If no error occured, the routine returns zero. Otherwise, in case of
290*  overflow of the sparse vector area, the routine returns non-zero. */
291
292int luf_enlarge_col(LUF *luf, int j, int cap)
293{     int n = luf->n;
294      int *vr_cap = luf->vr_cap;
295      int *vc_ptr = luf->vc_ptr;
296      int *vc_len = luf->vc_len;
297      int *vc_cap = luf->vc_cap;
298      int *sv_ind = luf->sv_ind;
299      double *sv_val = luf->sv_val;
300      int *sv_prev = luf->sv_prev;
301      int *sv_next = luf->sv_next;
302      int ret = 0;
303      int cur, k, kk;
304      xassert(1 <= j && j <= n);
305      xassert(vc_cap[j] < cap);
306      /* if there are less than cap free locations, defragment SVA */
307      if (luf->sv_end - luf->sv_beg < cap)
308      {  luf_defrag_sva(luf);
309         if (luf->sv_end - luf->sv_beg < cap)
310         {  ret = 1;
311            goto done;
312         }
313      }
314      /* save current capacity of the j-th column */
315      cur = vc_cap[j];
316      /* copy existing elements to the beginning of the free part */
317      memmove(&sv_ind[luf->sv_beg], &sv_ind[vc_ptr[j]],
318         vc_len[j] * sizeof(int));
319      memmove(&sv_val[luf->sv_beg], &sv_val[vc_ptr[j]],
320         vc_len[j] * sizeof(double));
321      /* set new pointer and new capacity of the j-th column */
322      vc_ptr[j] = luf->sv_beg;
323      vc_cap[j] = cap;
324      /* set new pointer to the beginning of the free part */
325      luf->sv_beg += cap;
326      /* now the j-th column starts in the rightmost location among
327         other rows and columns of the matrix V, so its node should be
328         moved to the end of the row/column linked list */
329      k = n + j;
330      /* remove the j-th column node from the linked list */
331      if (sv_prev[k] == 0)
332         luf->sv_head = sv_next[k];
333      else
334      {  /* capacity of the previous row/column can be increased at the
335            expense of old locations of the j-th column */
336         kk = sv_prev[k];
337         if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
338         sv_next[sv_prev[k]] = sv_next[k];
339      }
340      if (sv_next[k] == 0)
341         luf->sv_tail = sv_prev[k];
342      else
343         sv_prev[sv_next[k]] = sv_prev[k];
344      /* insert the j-th column node to the end of the linked list */
345      sv_prev[k] = luf->sv_tail;
346      sv_next[k] = 0;
347      if (sv_prev[k] == 0)
348         luf->sv_head = k;
349      else
350         sv_next[sv_prev[k]] = k;
351      luf->sv_tail = k;
352done: return ret;
353}
354
355/***********************************************************************
356*  reallocate - reallocate LU-factorization arrays
357*
358*  This routine reallocates arrays, whose size depends of n, the order
359*  of the matrix A to be factorized. */
360
361static void reallocate(LUF *luf, int n)
362{     int n_max = luf->n_max;
363      luf->n = n;
364      if (n <= n_max) goto done;
365      if (luf->fr_ptr != NULL) xfree(luf->fr_ptr);
366      if (luf->fr_len != NULL) xfree(luf->fr_len);
367      if (luf->fc_ptr != NULL) xfree(luf->fc_ptr);
368      if (luf->fc_len != NULL) xfree(luf->fc_len);
369      if (luf->vr_ptr != NULL) xfree(luf->vr_ptr);
370      if (luf->vr_len != NULL) xfree(luf->vr_len);
371      if (luf->vr_cap != NULL) xfree(luf->vr_cap);
372      if (luf->vr_piv != NULL) xfree(luf->vr_piv);
373      if (luf->vc_ptr != NULL) xfree(luf->vc_ptr);
374      if (luf->vc_len != NULL) xfree(luf->vc_len);
375      if (luf->vc_cap != NULL) xfree(luf->vc_cap);
376      if (luf->pp_row != NULL) xfree(luf->pp_row);
377      if (luf->pp_col != NULL) xfree(luf->pp_col);
378      if (luf->qq_row != NULL) xfree(luf->qq_row);
379      if (luf->qq_col != NULL) xfree(luf->qq_col);
380      if (luf->sv_prev != NULL) xfree(luf->sv_prev);
381      if (luf->sv_next != NULL) xfree(luf->sv_next);
382      if (luf->vr_max != NULL) xfree(luf->vr_max);
383      if (luf->rs_head != NULL) xfree(luf->rs_head);
384      if (luf->rs_prev != NULL) xfree(luf->rs_prev);
385      if (luf->rs_next != NULL) xfree(luf->rs_next);
386      if (luf->cs_head != NULL) xfree(luf->cs_head);
387      if (luf->cs_prev != NULL) xfree(luf->cs_prev);
388      if (luf->cs_next != NULL) xfree(luf->cs_next);
389      if (luf->flag != NULL) xfree(luf->flag);
390      if (luf->work != NULL) xfree(luf->work);
391      luf->n_max = n_max = n + 100;
392      luf->fr_ptr = xcalloc(1+n_max, sizeof(int));
393      luf->fr_len = xcalloc(1+n_max, sizeof(int));
394      luf->fc_ptr = xcalloc(1+n_max, sizeof(int));
395      luf->fc_len = xcalloc(1+n_max, sizeof(int));
396      luf->vr_ptr = xcalloc(1+n_max, sizeof(int));
397      luf->vr_len = xcalloc(1+n_max, sizeof(int));
398      luf->vr_cap = xcalloc(1+n_max, sizeof(int));
399      luf->vr_piv = xcalloc(1+n_max, sizeof(double));
400      luf->vc_ptr = xcalloc(1+n_max, sizeof(int));
401      luf->vc_len = xcalloc(1+n_max, sizeof(int));
402      luf->vc_cap = xcalloc(1+n_max, sizeof(int));
403      luf->pp_row = xcalloc(1+n_max, sizeof(int));
404      luf->pp_col = xcalloc(1+n_max, sizeof(int));
405      luf->qq_row = xcalloc(1+n_max, sizeof(int));
406      luf->qq_col = xcalloc(1+n_max, sizeof(int));
407      luf->sv_prev = xcalloc(1+n_max+n_max, sizeof(int));
408      luf->sv_next = xcalloc(1+n_max+n_max, sizeof(int));
409      luf->vr_max = xcalloc(1+n_max, sizeof(double));
410      luf->rs_head = xcalloc(1+n_max, sizeof(int));
411      luf->rs_prev = xcalloc(1+n_max, sizeof(int));
412      luf->rs_next = xcalloc(1+n_max, sizeof(int));
413      luf->cs_head = xcalloc(1+n_max, sizeof(int));
414      luf->cs_prev = xcalloc(1+n_max, sizeof(int));
415      luf->cs_next = xcalloc(1+n_max, sizeof(int));
416      luf->flag = xcalloc(1+n_max, sizeof(int));
417      luf->work = xcalloc(1+n_max, sizeof(double));
418done: return;
419}
420
421/***********************************************************************
422*  initialize - initialize LU-factorization data structures
423*
424*  This routine initializes data structures for subsequent computing
425*  the LU-factorization of a given matrix A, which is specified by the
426*  formal routine col. On exit V = A and F = P = Q = I, where I is the
427*  unity matrix. (Row-wise representation of the matrix F is not used
428*  at the factorization stage and therefore is not initialized.)
429*
430*  If no error occured, the routine returns zero. Otherwise, in case of
431*  overflow of the sparse vector area, the routine returns non-zero. */
432
433static int initialize(LUF *luf, int (*col)(void *info, int j, int rn[],
434      double aj[]), void *info)
435{     int n = luf->n;
436      int *fc_ptr = luf->fc_ptr;
437      int *fc_len = luf->fc_len;
438      int *vr_ptr = luf->vr_ptr;
439      int *vr_len = luf->vr_len;
440      int *vr_cap = luf->vr_cap;
441      int *vc_ptr = luf->vc_ptr;
442      int *vc_len = luf->vc_len;
443      int *vc_cap = luf->vc_cap;
444      int *pp_row = luf->pp_row;
445      int *pp_col = luf->pp_col;
446      int *qq_row = luf->qq_row;
447      int *qq_col = luf->qq_col;
448      int *sv_ind = luf->sv_ind;
449      double *sv_val = luf->sv_val;
450      int *sv_prev = luf->sv_prev;
451      int *sv_next = luf->sv_next;
452      double *vr_max = luf->vr_max;
453      int *rs_head = luf->rs_head;
454      int *rs_prev = luf->rs_prev;
455      int *rs_next = luf->rs_next;
456      int *cs_head = luf->cs_head;
457      int *cs_prev = luf->cs_prev;
458      int *cs_next = luf->cs_next;
459      int *flag = luf->flag;
460      double *work = luf->work;
461      int ret = 0;
462      int i, i_ptr, j, j_beg, j_end, k, len, nnz, sv_beg, sv_end, ptr;
463      double big, val;
464      /* free all locations of the sparse vector area */
465      sv_beg = 1;
466      sv_end = luf->sv_size + 1;
467      /* (row-wise representation of the matrix F is not initialized,
468         because it is not used at the factorization stage) */
469      /* build the matrix F in column-wise format (initially F = I) */
470      for (j = 1; j <= n; j++)
471      {  fc_ptr[j] = sv_end;
472         fc_len[j] = 0;
473      }
474      /* clear rows of the matrix V; clear the flag array */
475      for (i = 1; i <= n; i++)
476         vr_len[i] = vr_cap[i] = 0, flag[i] = 0;
477      /* build the matrix V in column-wise format (initially V = A);
478         count non-zeros in rows of this matrix; count total number of
479         non-zeros; compute largest of absolute values of elements */
480      nnz = 0;
481      big = 0.0;
482      for (j = 1; j <= n; j++)
483      {  int *rn = pp_row;
484         double *aj = work;
485         /* obtain j-th column of the matrix A */
486         len = col(info, j, rn, aj);
487         if (!(0 <= len && len <= n))
488            xfault("luf_factorize: j = %d; len = %d; invalid column len"
489               "gth\n", j, len);
490         /* check for free locations */
491         if (sv_end - sv_beg < len)
492         {  /* overflow of the sparse vector area */
493            ret = 1;
494            goto done;
495         }
496         /* set pointer to the j-th column */
497         vc_ptr[j] = sv_beg;
498         /* set length of the j-th column */
499         vc_len[j] = vc_cap[j] = len;
500         /* count total number of non-zeros */
501         nnz += len;
502         /* walk through elements of the j-th column */
503         for (ptr = 1; ptr <= len; ptr++)
504         {  /* get row index and numerical value of a[i,j] */
505            i = rn[ptr];
506            val = aj[ptr];
507            if (!(1 <= i && i <= n))
508               xfault("luf_factorize: i = %d; j = %d; invalid row index"
509                  "\n", i, j);
510            if (flag[i])
511               xfault("luf_factorize: i = %d; j = %d; duplicate element"
512                  " not allowed\n", i, j);
513            if (val == 0.0)
514               xfault("luf_factorize: i = %d; j = %d; zero element not "
515                  "allowed\n", i, j);
516            /* add new element v[i,j] = a[i,j] to j-th column */
517            sv_ind[sv_beg] = i;
518            sv_val[sv_beg] = val;
519            sv_beg++;
520            /* big := max(big, |a[i,j]|) */
521            if (val < 0.0) val = - val;
522            if (big < val) big = val;
523            /* mark non-zero in the i-th position of the j-th column */
524            flag[i] = 1;
525            /* increase length of the i-th row */
526            vr_cap[i]++;
527         }
528         /* reset all non-zero marks */
529         for (ptr = 1; ptr <= len; ptr++) flag[rn[ptr]] = 0;
530      }
531      /* allocate rows of the matrix V */
532      for (i = 1; i <= n; i++)
533      {  /* get length of the i-th row */
534         len = vr_cap[i];
535         /* check for free locations */
536         if (sv_end - sv_beg < len)
537         {  /* overflow of the sparse vector area */
538            ret = 1;
539            goto done;
540         }
541         /* set pointer to the i-th row */
542         vr_ptr[i] = sv_beg;
543         /* reserve locations for the i-th row */
544         sv_beg += len;
545      }
546      /* build the matrix V in row-wise format using representation of
547         this matrix in column-wise format */
548      for (j = 1; j <= n; j++)
549      {  /* walk through elements of the j-th column */
550         j_beg = vc_ptr[j];
551         j_end = j_beg + vc_len[j] - 1;
552         for (k = j_beg; k <= j_end; k++)
553         {  /* get row index and numerical value of v[i,j] */
554            i = sv_ind[k];
555            val = sv_val[k];
556            /* store element in the i-th row */
557            i_ptr = vr_ptr[i] + vr_len[i];
558            sv_ind[i_ptr] = j;
559            sv_val[i_ptr] = val;
560            /* increase count of the i-th row */
561            vr_len[i]++;
562         }
563      }
564      /* initialize the matrices P and Q (initially P = Q = I) */
565      for (k = 1; k <= n; k++)
566         pp_row[k] = pp_col[k] = qq_row[k] = qq_col[k] = k;
567      /* set sva partitioning pointers */
568      luf->sv_beg = sv_beg;
569      luf->sv_end = sv_end;
570      /* the initial physical order of rows and columns of the matrix V
571         is n+1, ..., n+n, 1, ..., n (firstly columns, then rows) */
572      luf->sv_head = n+1;
573      luf->sv_tail = n;
574      for (i = 1; i <= n; i++)
575      {  sv_prev[i] = i-1;
576         sv_next[i] = i+1;
577      }
578      sv_prev[1] = n+n;
579      sv_next[n] = 0;
580      for (j = 1; j <= n; j++)
581      {  sv_prev[n+j] = n+j-1;
582         sv_next[n+j] = n+j+1;
583      }
584      sv_prev[n+1] = 0;
585      sv_next[n+n] = 1;
586      /* clear working arrays */
587      for (k = 1; k <= n; k++)
588      {  flag[k] = 0;
589         work[k] = 0.0;
590      }
591      /* initialize some statistics */
592      luf->nnz_a = nnz;
593      luf->nnz_f = 0;
594      luf->nnz_v = nnz;
595      luf->max_a = big;
596      luf->big_v = big;
597      luf->rank = -1;
598      /* initially the active submatrix is the entire matrix V */
599      /* largest of absolute values of elements in each active row is
600         unknown yet */
601      for (i = 1; i <= n; i++) vr_max[i] = -1.0;
602      /* build linked lists of active rows */
603      for (len = 0; len <= n; len++) rs_head[len] = 0;
604      for (i = 1; i <= n; i++)
605      {  len = vr_len[i];
606         rs_prev[i] = 0;
607         rs_next[i] = rs_head[len];
608         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
609         rs_head[len] = i;
610      }
611      /* build linked lists of active columns */
612      for (len = 0; len <= n; len++) cs_head[len] = 0;
613      for (j = 1; j <= n; j++)
614      {  len = vc_len[j];
615         cs_prev[j] = 0;
616         cs_next[j] = cs_head[len];
617         if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
618         cs_head[len] = j;
619      }
620done: /* return to the factorizing routine */
621      return ret;
622}
623
624/***********************************************************************
625*  find_pivot - choose a pivot element
626*
627*  This routine chooses a pivot element in the active submatrix of the
628*  matrix U = P*V*Q.
629*
630*  It is assumed that on entry the matrix U has the following partially
631*  triangularized form:
632*
633*        1       k         n
634*     1  x x x x x x x x x x
635*        . x x x x x x x x x
636*        . . x x x x x x x x
637*        . . . x x x x x x x
638*     k  . . . . * * * * * *
639*        . . . . * * * * * *
640*        . . . . * * * * * *
641*        . . . . * * * * * *
642*        . . . . * * * * * *
643*     n  . . . . * * * * * *
644*
645*  where rows and columns k, k+1, ..., n belong to the active submatrix
646*  (elements of the active submatrix are marked by '*').
647*
648*  Since the matrix U = P*V*Q is not stored, the routine works with the
649*  matrix V. It is assumed that the row-wise representation corresponds
650*  to the matrix V, but the column-wise representation corresponds to
651*  the active submatrix of the matrix V, i.e. elements of the matrix V,
652*  which doesn't belong to the active submatrix, are missing from the
653*  column linked lists. It is also assumed that each active row of the
654*  matrix V is in the set R[len], where len is number of non-zeros in
655*  the row, and each active column of the matrix V is in the set C[len],
656*  where len is number of non-zeros in the column (in the latter case
657*  only elements of the active submatrix are counted; such elements are
658*  marked by '*' on the figure above).
659*
660*  For the reason of numerical stability the routine applies so called
661*  threshold pivoting proposed by J.Reid. It is assumed that an element
662*  v[i,j] can be selected as a pivot candidate if it is not very small
663*  (in absolute value) among other elements in the same row, i.e. if it
664*  satisfies to the stability condition |v[i,j]| >= tol * max|v[i,*]|,
665*  where 0 < tol < 1 is a given tolerance.
666*
667*  In order to keep sparsity of the matrix V the routine uses Markowitz
668*  strategy, trying to choose such element v[p,q], which satisfies to
669*  the stability condition (see above) and has smallest Markowitz cost
670*  (nr[p]-1) * (nc[q]-1), where nr[p] and nc[q] are numbers of non-zero
671*  elements, respectively, in the p-th row and in the q-th column of the
672*  active submatrix.
673*
674*  In order to reduce the search, i.e. not to walk through all elements
675*  of the active submatrix, the routine exploits a technique proposed by
676*  I.Duff. This technique is based on using the sets R[len] and C[len]
677*  of active rows and columns.
678*
679*  If the pivot element v[p,q] has been chosen, the routine stores its
680*  indices to the locations *p and *q and returns zero. Otherwise, if
681*  the active submatrix is empty and therefore the pivot element can't
682*  be chosen, the routine returns non-zero. */
683
684static int find_pivot(LUF *luf, int *_p, int *_q)
685{     int n = luf->n;
686      int *vr_ptr = luf->vr_ptr;
687      int *vr_len = luf->vr_len;
688      int *vc_ptr = luf->vc_ptr;
689      int *vc_len = luf->vc_len;
690      int *sv_ind = luf->sv_ind;
691      double *sv_val = luf->sv_val;
692      double *vr_max = luf->vr_max;
693      int *rs_head = luf->rs_head;
694      int *rs_next = luf->rs_next;
695      int *cs_head = luf->cs_head;
696      int *cs_prev = luf->cs_prev;
697      int *cs_next = luf->cs_next;
698      double piv_tol = luf->piv_tol;
699      int piv_lim = luf->piv_lim;
700      int suhl = luf->suhl;
701      int p, q, len, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr,
702         ncand, next_j, min_p, min_q, min_len;
703      double best, cost, big, temp;
704      /* initially no pivot candidates have been found so far */
705      p = q = 0, best = DBL_MAX, ncand = 0;
706      /* if in the active submatrix there is a column that has the only
707         non-zero (column singleton), choose it as pivot */
708      j = cs_head[1];
709      if (j != 0)
710      {  xassert(vc_len[j] == 1);
711         p = sv_ind[vc_ptr[j]], q = j;
712         goto done;
713      }
714      /* if in the active submatrix there is a row that has the only
715         non-zero (row singleton), choose it as pivot */
716      i = rs_head[1];
717      if (i != 0)
718      {  xassert(vr_len[i] == 1);
719         p = i, q = sv_ind[vr_ptr[i]];
720         goto done;
721      }
722      /* there are no singletons in the active submatrix; walk through
723         other non-empty rows and columns */
724      for (len = 2; len <= n; len++)
725      {  /* consider active columns that have len non-zeros */
726         for (j = cs_head[len]; j != 0; j = next_j)
727         {  /* the j-th column has len non-zeros */
728            j_beg = vc_ptr[j];
729            j_end = j_beg + vc_len[j] - 1;
730            /* save pointer to the next column with the same length */
731            next_j = cs_next[j];
732            /* find an element in the j-th column, which is placed in a
733               row with minimal number of non-zeros and satisfies to the
734               stability condition (such element may not exist) */
735            min_p = min_q = 0, min_len = INT_MAX;
736            for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
737            {  /* get row index of v[i,j] */
738               i = sv_ind[j_ptr];
739               i_beg = vr_ptr[i];
740               i_end = i_beg + vr_len[i] - 1;
741               /* if the i-th row is not shorter than that one, where
742                  minimal element is currently placed, skip v[i,j] */
743               if (vr_len[i] >= min_len) continue;
744               /* determine the largest of absolute values of elements
745                  in the i-th row */
746               big = vr_max[i];
747               if (big < 0.0)
748               {  /* the largest value is unknown yet; compute it */
749                  for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
750                  {  temp = sv_val[i_ptr];
751                     if (temp < 0.0) temp = - temp;
752                     if (big < temp) big = temp;
753                  }
754                  vr_max[i] = big;
755               }
756               /* find v[i,j] in the i-th row */
757               for (i_ptr = vr_ptr[i]; sv_ind[i_ptr] != j; i_ptr++);
758               xassert(i_ptr <= i_end);
759               /* if v[i,j] doesn't satisfy to the stability condition,
760                  skip it */
761               temp = sv_val[i_ptr];
762               if (temp < 0.0) temp = - temp;
763               if (temp < piv_tol * big) continue;
764               /* v[i,j] is better than the current minimal element */
765               min_p = i, min_q = j, min_len = vr_len[i];
766               /* if Markowitz cost of the current minimal element is
767                  not greater than (len-1)**2, it can be chosen right
768                  now; this heuristic reduces the search and works well
769                  in many cases */
770               if (min_len <= len)
771               {  p = min_p, q = min_q;
772                  goto done;
773               }
774            }
775            /* the j-th column has been scanned */
776            if (min_p != 0)
777            {  /* the minimal element is a next pivot candidate */
778               ncand++;
779               /* compute its Markowitz cost */
780               cost = (double)(min_len - 1) * (double)(len - 1);
781               /* choose between the minimal element and the current
782                  candidate */
783               if (cost < best) p = min_p, q = min_q, best = cost;
784               /* if piv_lim candidates have been considered, there are
785                  doubts that a much better candidate exists; therefore
786                  it's time to terminate the search */
787               if (ncand == piv_lim) goto done;
788            }
789            else
790            {  /* the j-th column has no elements, which satisfy to the
791                  stability condition; Uwe Suhl suggests to exclude such
792                  column from the further consideration until it becomes
793                  a column singleton; in hard cases this significantly
794                  reduces a time needed for pivot searching */
795               if (suhl)
796               {  /* remove the j-th column from the active set */
797                  if (cs_prev[j] == 0)
798                     cs_head[len] = cs_next[j];
799                  else
800                     cs_next[cs_prev[j]] = cs_next[j];
801                  if (cs_next[j] == 0)
802                     /* nop */;
803                  else
804                     cs_prev[cs_next[j]] = cs_prev[j];
805                  /* the following assignment is used to avoid an error
806                     when the routine eliminate (see below) will try to
807                     remove the j-th column from the active set */
808                  cs_prev[j] = cs_next[j] = j;
809               }
810            }
811         }
812         /* consider active rows that have len non-zeros */
813         for (i = rs_head[len]; i != 0; i = rs_next[i])
814         {  /* the i-th row has len non-zeros */
815            i_beg = vr_ptr[i];
816            i_end = i_beg + vr_len[i] - 1;
817            /* determine the largest of absolute values of elements in
818               the i-th row */
819            big = vr_max[i];
820            if (big < 0.0)
821            {  /* the largest value is unknown yet; compute it */
822               for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
823               {  temp = sv_val[i_ptr];
824                  if (temp < 0.0) temp = - temp;
825                  if (big < temp) big = temp;
826               }
827               vr_max[i] = big;
828            }
829            /* find an element in the i-th row, which is placed in a
830               column with minimal number of non-zeros and satisfies to
831               the stability condition (such element always exists) */
832            min_p = min_q = 0, min_len = INT_MAX;
833            for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
834            {  /* get column index of v[i,j] */
835               j = sv_ind[i_ptr];
836               /* if the j-th column is not shorter than that one, where
837                  minimal element is currently placed, skip v[i,j] */
838               if (vc_len[j] >= min_len) continue;
839               /* if v[i,j] doesn't satisfy to the stability condition,
840                  skip it */
841               temp = sv_val[i_ptr];
842               if (temp < 0.0) temp = - temp;
843               if (temp < piv_tol * big) continue;
844               /* v[i,j] is better than the current minimal element */
845               min_p = i, min_q = j, min_len = vc_len[j];
846               /* if Markowitz cost of the current minimal element is
847                  not greater than (len-1)**2, it can be chosen right
848                  now; this heuristic reduces the search and works well
849                  in many cases */
850               if (min_len <= len)
851               {  p = min_p, q = min_q;
852                  goto done;
853               }
854            }
855            /* the i-th row has been scanned */
856            if (min_p != 0)
857            {  /* the minimal element is a next pivot candidate */
858               ncand++;
859               /* compute its Markowitz cost */
860               cost = (double)(len - 1) * (double)(min_len - 1);
861               /* choose between the minimal element and the current
862                  candidate */
863               if (cost < best) p = min_p, q = min_q, best = cost;
864               /* if piv_lim candidates have been considered, there are
865                  doubts that a much better candidate exists; therefore
866                  it's time to terminate the search */
867               if (ncand == piv_lim) goto done;
868            }
869            else
870            {  /* this can't be because this can never be */
871               xassert(min_p != min_p);
872            }
873         }
874      }
875done: /* bring the pivot to the factorizing routine */
876      *_p = p, *_q = q;
877      return (p == 0);
878}
879
880/***********************************************************************
881*  eliminate - perform gaussian elimination.
882*
883*  This routine performs elementary gaussian transformations in order
884*  to eliminate subdiagonal elements in the k-th column of the matrix
885*  U = P*V*Q using the pivot element u[k,k], where k is the number of
886*  the current elimination step.
887*
888*  The parameters p and q are, respectively, row and column indices of
889*  the element v[p,q], which corresponds to the element u[k,k].
890*
891*  Each time when the routine applies the elementary transformation to
892*  a non-pivot row of the matrix V, it stores the corresponding element
893*  to the matrix F in order to keep the main equality A = F*V.
894*
895*  The routine assumes that on entry the matrices L = P*F*inv(P) and
896*  U = P*V*Q are the following:
897*
898*        1       k                  1       k         n
899*     1  1 . . . . . . . . .     1  x x x x x x x x x x
900*        x 1 . . . . . . . .        . x x x x x x x x x
901*        x x 1 . . . . . . .        . . x x x x x x x x
902*        x x x 1 . . . . . .        . . . x x x x x x x
903*     k  x x x x 1 . . . . .     k  . . . . * * * * * *
904*        x x x x _ 1 . . . .        . . . . # * * * * *
905*        x x x x _ . 1 . . .        . . . . # * * * * *
906*        x x x x _ . . 1 . .        . . . . # * * * * *
907*        x x x x _ . . . 1 .        . . . . # * * * * *
908*     n  x x x x _ . . . . 1     n  . . . . # * * * * *
909*
910*             matrix L                   matrix U
911*
912*  where rows and columns of the matrix U with numbers k, k+1, ..., n
913*  form the active submatrix (eliminated elements are marked by '#' and
914*  other elements of the active submatrix are marked by '*'). Note that
915*  each eliminated non-zero element u[i,k] of the matrix U gives the
916*  corresponding element l[i,k] of the matrix L (marked by '_').
917*
918*  Actually all operations are performed on the matrix V. Should note
919*  that the row-wise representation corresponds to the matrix V, but the
920*  column-wise representation corresponds to the active submatrix of the
921*  matrix V, i.e. elements of the matrix V, which doesn't belong to the
922*  active submatrix, are missing from the column linked lists.
923*
924*  Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal
925*  elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies
926*  the following elementary gaussian transformations:
927*
928*     (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V),
929*
930*  where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier.
931*
932*  Additionally, in order to keep the main equality A = F*V, each time
933*  when the routine applies the transformation to i-th row of the matrix
934*  V, it also adds f[i,p] as a new element to the matrix F.
935*
936*  IMPORTANT: On entry the working arrays flag and work should contain
937*  zeros. This status is provided by the routine on exit.
938*
939*  If no error occured, the routine returns zero. Otherwise, in case of
940*  overflow of the sparse vector area, the routine returns non-zero. */
941
942static int eliminate(LUF *luf, int p, int q)
943{     int n = luf->n;
944      int *fc_ptr = luf->fc_ptr;
945      int *fc_len = luf->fc_len;
946      int *vr_ptr = luf->vr_ptr;
947      int *vr_len = luf->vr_len;
948      int *vr_cap = luf->vr_cap;
949      double *vr_piv = luf->vr_piv;
950      int *vc_ptr = luf->vc_ptr;
951      int *vc_len = luf->vc_len;
952      int *vc_cap = luf->vc_cap;
953      int *sv_ind = luf->sv_ind;
954      double *sv_val = luf->sv_val;
955      int *sv_prev = luf->sv_prev;
956      int *sv_next = luf->sv_next;
957      double *vr_max = luf->vr_max;
958      int *rs_head = luf->rs_head;
959      int *rs_prev = luf->rs_prev;
960      int *rs_next = luf->rs_next;
961      int *cs_head = luf->cs_head;
962      int *cs_prev = luf->cs_prev;
963      int *cs_next = luf->cs_next;
964      int *flag = luf->flag;
965      double *work = luf->work;
966      double eps_tol = luf->eps_tol;
967      /* at this stage the row-wise representation of the matrix F is
968         not used, so fr_len can be used as a working array */
969      int *ndx = luf->fr_len;
970      int ret = 0;
971      int len, fill, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, k,
972         p_beg, p_end, p_ptr, q_beg, q_end, q_ptr;
973      double fip, val, vpq, temp;
974      xassert(1 <= p && p <= n);
975      xassert(1 <= q && q <= n);
976      /* remove the p-th (pivot) row from the active set; this row will
977         never return there */
978      if (rs_prev[p] == 0)
979         rs_head[vr_len[p]] = rs_next[p];
980      else
981         rs_next[rs_prev[p]] = rs_next[p];
982      if (rs_next[p] == 0)
983         ;
984      else
985         rs_prev[rs_next[p]] = rs_prev[p];
986      /* remove the q-th (pivot) column from the active set; this column
987         will never return there */
988      if (cs_prev[q] == 0)
989         cs_head[vc_len[q]] = cs_next[q];
990      else
991         cs_next[cs_prev[q]] = cs_next[q];
992      if (cs_next[q] == 0)
993         ;
994      else
995         cs_prev[cs_next[q]] = cs_prev[q];
996      /* find the pivot v[p,q] = u[k,k] in the p-th row */
997      p_beg = vr_ptr[p];
998      p_end = p_beg + vr_len[p] - 1;
999      for (p_ptr = p_beg; sv_ind[p_ptr] != q; p_ptr++) /* nop */;
1000      xassert(p_ptr <= p_end);
1001      /* store value of the pivot */
1002      vpq = (vr_piv[p] = sv_val[p_ptr]);
1003      /* remove the pivot from the p-th row */
1004      sv_ind[p_ptr] = sv_ind[p_end];
1005      sv_val[p_ptr] = sv_val[p_end];
1006      vr_len[p]--;
1007      p_end--;
1008      /* find the pivot v[p,q] = u[k,k] in the q-th column */
1009      q_beg = vc_ptr[q];
1010      q_end = q_beg + vc_len[q] - 1;
1011      for (q_ptr = q_beg; sv_ind[q_ptr] != p; q_ptr++) /* nop */;
1012      xassert(q_ptr <= q_end);
1013      /* remove the pivot from the q-th column */
1014      sv_ind[q_ptr] = sv_ind[q_end];
1015      vc_len[q]--;
1016      q_end--;
1017      /* walk through the p-th (pivot) row, which doesn't contain the
1018         pivot v[p,q] already, and do the following... */
1019      for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
1020      {  /* get column index of v[p,j] */
1021         j = sv_ind[p_ptr];
1022         /* store v[p,j] to the working array */
1023         flag[j] = 1;
1024         work[j] = sv_val[p_ptr];
1025         /* remove the j-th column from the active set; this column will
1026            return there later with new length */
1027         if (cs_prev[j] == 0)
1028            cs_head[vc_len[j]] = cs_next[j];
1029         else
1030            cs_next[cs_prev[j]] = cs_next[j];
1031         if (cs_next[j] == 0)
1032            ;
1033         else
1034            cs_prev[cs_next[j]] = cs_prev[j];
1035         /* find v[p,j] in the j-th column */
1036         j_beg = vc_ptr[j];
1037         j_end = j_beg + vc_len[j] - 1;
1038         for (j_ptr = j_beg; sv_ind[j_ptr] != p; j_ptr++) /* nop */;
1039         xassert(j_ptr <= j_end);
1040         /* since v[p,j] leaves the active submatrix, remove it from the
1041            j-th column; however, v[p,j] is kept in the p-th row */
1042         sv_ind[j_ptr] = sv_ind[j_end];
1043         vc_len[j]--;
1044      }
1045      /* walk through the q-th (pivot) column, which doesn't contain the
1046         pivot v[p,q] already, and perform gaussian elimination */
1047      while (q_beg <= q_end)
1048      {  /* element v[i,q] should be eliminated */
1049         /* get row index of v[i,q] */
1050         i = sv_ind[q_beg];
1051         /* remove the i-th row from the active set; later this row will
1052            return there with new length */
1053         if (rs_prev[i] == 0)
1054            rs_head[vr_len[i]] = rs_next[i];
1055         else
1056            rs_next[rs_prev[i]] = rs_next[i];
1057         if (rs_next[i] == 0)
1058            ;
1059         else
1060            rs_prev[rs_next[i]] = rs_prev[i];
1061         /* find v[i,q] in the i-th row */
1062         i_beg = vr_ptr[i];
1063         i_end = i_beg + vr_len[i] - 1;
1064         for (i_ptr = i_beg; sv_ind[i_ptr] != q; i_ptr++) /* nop */;
1065         xassert(i_ptr <= i_end);
1066         /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] */
1067         fip = sv_val[i_ptr] / vpq;
1068         /* since v[i,q] should be eliminated, remove it from the i-th
1069            row */
1070         sv_ind[i_ptr] = sv_ind[i_end];
1071         sv_val[i_ptr] = sv_val[i_end];
1072         vr_len[i]--;
1073         i_end--;
1074         /* and from the q-th column */
1075         sv_ind[q_beg] = sv_ind[q_end];
1076         vc_len[q]--;
1077         q_end--;
1078         /* perform gaussian transformation:
1079            (i-th row) := (i-th row) - f[i,p] * (p-th row)
1080            note that now the p-th row, which is in the working array,
1081            doesn't contain the pivot v[p,q], and the i-th row doesn't
1082            contain the eliminated element v[i,q] */
1083         /* walk through the i-th row and transform existing non-zero
1084            elements */
1085         fill = vr_len[p];
1086         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
1087         {  /* get column index of v[i,j] */
1088            j = sv_ind[i_ptr];
1089            /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
1090            if (flag[j])
1091            {  /* v[p,j] != 0 */
1092               temp = (sv_val[i_ptr] -= fip * work[j]);
1093               if (temp < 0.0) temp = - temp;
1094               flag[j] = 0;
1095               fill--; /* since both v[i,j] and v[p,j] exist */
1096               if (temp == 0.0 || temp < eps_tol)
1097               {  /* new v[i,j] is closer to zero; replace it by exact
1098                     zero, i.e. remove it from the active submatrix */
1099                  /* remove v[i,j] from the i-th row */
1100                  sv_ind[i_ptr] = sv_ind[i_end];
1101                  sv_val[i_ptr] = sv_val[i_end];
1102                  vr_len[i]--;
1103                  i_ptr--;
1104                  i_end--;
1105                  /* find v[i,j] in the j-th column */
1106                  j_beg = vc_ptr[j];
1107                  j_end = j_beg + vc_len[j] - 1;
1108                  for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++);
1109                  xassert(j_ptr <= j_end);
1110                  /* remove v[i,j] from the j-th column */
1111                  sv_ind[j_ptr] = sv_ind[j_end];
1112                  vc_len[j]--;
1113               }
1114               else
1115               {  /* v_big := max(v_big, |v[i,j]|) */
1116                  if (luf->big_v < temp) luf->big_v = temp;
1117               }
1118            }
1119         }
1120         /* now flag is the pattern of the set v[p,*] \ v[i,*], and fill
1121            is number of non-zeros in this set; therefore up to fill new
1122            non-zeros may appear in the i-th row */
1123         if (vr_len[i] + fill > vr_cap[i])
1124         {  /* enlarge the i-th row */
1125            if (luf_enlarge_row(luf, i, vr_len[i] + fill))
1126            {  /* overflow of the sparse vector area */
1127               ret = 1;
1128               goto done;
1129            }
1130            /* defragmentation may change row and column pointers of the
1131               matrix V */
1132            p_beg = vr_ptr[p];
1133            p_end = p_beg + vr_len[p] - 1;
1134            q_beg = vc_ptr[q];
1135            q_end = q_beg + vc_len[q] - 1;
1136         }
1137         /* walk through the p-th (pivot) row and create new elements
1138            of the i-th row that appear due to fill-in; column indices
1139            of these new elements are accumulated in the array ndx */
1140         len = 0;
1141         for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
1142         {  /* get column index of v[p,j], which may cause fill-in */
1143            j = sv_ind[p_ptr];
1144            if (flag[j])
1145            {  /* compute new non-zero v[i,j] = 0 - f[i,p] * v[p,j] */
1146               temp = (val = - fip * work[j]);
1147               if (temp < 0.0) temp = - temp;
1148               if (temp == 0.0 || temp < eps_tol)
1149                  /* if v[i,j] is closer to zero; just ignore it */;
1150               else
1151               {  /* add v[i,j] to the i-th row */
1152                  i_ptr = vr_ptr[i] + vr_len[i];
1153                  sv_ind[i_ptr] = j;
1154                  sv_val[i_ptr] = val;
1155                  vr_len[i]++;
1156                  /* remember column index of v[i,j] */
1157                  ndx[++len] = j;
1158                  /* big_v := max(big_v, |v[i,j]|) */
1159                  if (luf->big_v < temp) luf->big_v = temp;
1160               }
1161            }
1162            else
1163            {  /* there is no fill-in, because v[i,j] already exists in
1164                  the i-th row; restore the flag of the element v[p,j],
1165                  which was reset before */
1166               flag[j] = 1;
1167            }
1168         }
1169         /* add new non-zeros v[i,j] to the corresponding columns */
1170         for (k = 1; k <= len; k++)
1171         {  /* get column index of new non-zero v[i,j] */
1172            j = ndx[k];
1173            /* one free location is needed in the j-th column */
1174            if (vc_len[j] + 1 > vc_cap[j])
1175            {  /* enlarge the j-th column */
1176               if (luf_enlarge_col(luf, j, vc_len[j] + 10))
1177               {  /* overflow of the sparse vector area */
1178                  ret = 1;
1179                  goto done;
1180               }
1181               /* defragmentation may change row and column pointers of
1182                  the matrix V */
1183               p_beg = vr_ptr[p];
1184               p_end = p_beg + vr_len[p] - 1;
1185               q_beg = vc_ptr[q];
1186               q_end = q_beg + vc_len[q] - 1;
1187            }
1188            /* add new non-zero v[i,j] to the j-th column */
1189            j_ptr = vc_ptr[j] + vc_len[j];
1190            sv_ind[j_ptr] = i;
1191            vc_len[j]++;
1192         }
1193         /* now the i-th row has been completely transformed, therefore
1194            it can return to the active set with new length */
1195         rs_prev[i] = 0;
1196         rs_next[i] = rs_head[vr_len[i]];
1197         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
1198         rs_head[vr_len[i]] = i;
1199         /* the largest of absolute values of elements in the i-th row
1200            is currently unknown */
1201         vr_max[i] = -1.0;
1202         /* at least one free location is needed to store the gaussian
1203            multiplier */
1204         if (luf->sv_end - luf->sv_beg < 1)
1205         {  /* there are no free locations at all; defragment SVA */
1206            luf_defrag_sva(luf);
1207            if (luf->sv_end - luf->sv_beg < 1)
1208            {  /* overflow of the sparse vector area */
1209               ret = 1;
1210               goto done;
1211            }
1212            /* defragmentation may change row and column pointers of the
1213               matrix V */
1214            p_beg = vr_ptr[p];
1215            p_end = p_beg + vr_len[p] - 1;
1216            q_beg = vc_ptr[q];
1217            q_end = q_beg + vc_len[q] - 1;
1218         }
1219         /* add the element f[i,p], which is the gaussian multiplier,
1220            to the matrix F */
1221         luf->sv_end--;
1222         sv_ind[luf->sv_end] = i;
1223         sv_val[luf->sv_end] = fip;
1224         fc_len[p]++;
1225         /* end of elimination loop */
1226      }
1227      /* at this point the q-th (pivot) column should be empty */
1228      xassert(vc_len[q] == 0);
1229      /* reset capacity of the q-th column */
1230      vc_cap[q] = 0;
1231      /* remove node of the q-th column from the addressing list */
1232      k = n + q;
1233      if (sv_prev[k] == 0)
1234         luf->sv_head = sv_next[k];
1235      else
1236         sv_next[sv_prev[k]] = sv_next[k];
1237      if (sv_next[k] == 0)
1238         luf->sv_tail = sv_prev[k];
1239      else
1240         sv_prev[sv_next[k]] = sv_prev[k];
1241      /* the p-th column of the matrix F has been completely built; set
1242         its pointer */
1243      fc_ptr[p] = luf->sv_end;
1244      /* walk through the p-th (pivot) row and do the following... */
1245      for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
1246      {  /* get column index of v[p,j] */
1247         j = sv_ind[p_ptr];
1248         /* erase v[p,j] from the working array */
1249         flag[j] = 0;
1250         work[j] = 0.0;
1251         /* the j-th column has been completely transformed, therefore
1252            it can return to the active set with new length; however
1253            the special case c_prev[j] = c_next[j] = j means that the
1254            routine find_pivot excluded the j-th column from the active
1255            set due to Uwe Suhl's rule, and therefore in this case the
1256            column can return to the active set only if it is a column
1257            singleton */
1258         if (!(vc_len[j] != 1 && cs_prev[j] == j && cs_next[j] == j))
1259         {  cs_prev[j] = 0;
1260            cs_next[j] = cs_head[vc_len[j]];
1261            if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
1262            cs_head[vc_len[j]] = j;
1263         }
1264      }
1265done: /* return to the factorizing routine */
1266      return ret;
1267}
1268
1269/***********************************************************************
1270*  build_v_cols - build the matrix V in column-wise format
1271*
1272*  This routine builds the column-wise representation of the matrix V
1273*  using its row-wise representation.
1274*
1275*  If no error occured, the routine returns zero. Otherwise, in case of
1276*  overflow of the sparse vector area, the routine returns non-zero. */
1277
1278static int build_v_cols(LUF *luf)
1279{     int n = luf->n;
1280      int *vr_ptr = luf->vr_ptr;
1281      int *vr_len = luf->vr_len;
1282      int *vc_ptr = luf->vc_ptr;
1283      int *vc_len = luf->vc_len;
1284      int *vc_cap = luf->vc_cap;
1285      int *sv_ind = luf->sv_ind;
1286      double *sv_val = luf->sv_val;
1287      int *sv_prev = luf->sv_prev;
1288      int *sv_next = luf->sv_next;
1289      int ret = 0;
1290      int i, i_beg, i_end, i_ptr, j, j_ptr, k, nnz;
1291      /* it is assumed that on entry all columns of the matrix V are
1292         empty, i.e. vc_len[j] = vc_cap[j] = 0 for all j = 1, ..., n,
1293         and have been removed from the addressing list */
1294      /* count non-zeros in columns of the matrix V; count total number
1295         of non-zeros in this matrix */
1296      nnz = 0;
1297      for (i = 1; i <= n; i++)
1298      {  /* walk through elements of the i-th row and count non-zeros
1299            in the corresponding columns */
1300         i_beg = vr_ptr[i];
1301         i_end = i_beg + vr_len[i] - 1;
1302         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
1303            vc_cap[sv_ind[i_ptr]]++;
1304         /* count total number of non-zeros */
1305         nnz += vr_len[i];
1306      }
1307      /* store total number of non-zeros */
1308      luf->nnz_v = nnz;
1309      /* check for free locations */
1310      if (luf->sv_end - luf->sv_beg < nnz)
1311      {  /* overflow of the sparse vector area */
1312         ret = 1;
1313         goto done;
1314      }
1315      /* allocate columns of the matrix V */
1316      for (j = 1; j <= n; j++)
1317      {  /* set pointer to the j-th column */
1318         vc_ptr[j] = luf->sv_beg;
1319         /* reserve locations for the j-th column */
1320         luf->sv_beg += vc_cap[j];
1321      }
1322      /* build the matrix V in column-wise format using this matrix in
1323         row-wise format */
1324      for (i = 1; i <= n; i++)
1325      {  /* walk through elements of the i-th row */
1326         i_beg = vr_ptr[i];
1327         i_end = i_beg + vr_len[i] - 1;
1328         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
1329         {  /* get column index */
1330            j = sv_ind[i_ptr];
1331            /* store element in the j-th column */
1332            j_ptr = vc_ptr[j] + vc_len[j];
1333            sv_ind[j_ptr] = i;
1334            sv_val[j_ptr] = sv_val[i_ptr];
1335            /* increase length of the j-th column */
1336            vc_len[j]++;
1337         }
1338      }
1339      /* now columns are placed in the sparse vector area behind rows
1340         in the order n+1, n+2, ..., n+n; so insert column nodes in the
1341         addressing list using this order */
1342      for (k = n+1; k <= n+n; k++)
1343      {  sv_prev[k] = k-1;
1344         sv_next[k] = k+1;
1345      }
1346      sv_prev[n+1] = luf->sv_tail;
1347      sv_next[luf->sv_tail] = n+1;
1348      sv_next[n+n] = 0;
1349      luf->sv_tail = n+n;
1350done: /* return to the factorizing routine */
1351      return ret;
1352}
1353
1354/***********************************************************************
1355*  build_f_rows - build the matrix F in row-wise format
1356*
1357*  This routine builds the row-wise representation of the matrix F using
1358*  its column-wise representation.
1359*
1360*  If no error occured, the routine returns zero. Otherwise, in case of
1361*  overflow of the sparse vector area, the routine returns non-zero. */
1362
1363static int build_f_rows(LUF *luf)
1364{     int n = luf->n;
1365      int *fr_ptr = luf->fr_ptr;
1366      int *fr_len = luf->fr_len;
1367      int *fc_ptr = luf->fc_ptr;
1368      int *fc_len = luf->fc_len;
1369      int *sv_ind = luf->sv_ind;
1370      double *sv_val = luf->sv_val;
1371      int ret = 0;
1372      int i, j, j_beg, j_end, j_ptr, ptr, nnz;
1373      /* clear rows of the matrix F */
1374      for (i = 1; i <= n; i++) fr_len[i] = 0;
1375      /* count non-zeros in rows of the matrix F; count total number of
1376         non-zeros in this matrix */
1377      nnz = 0;
1378      for (j = 1; j <= n; j++)
1379      {  /* walk through elements of the j-th column and count non-zeros
1380            in the corresponding rows */
1381         j_beg = fc_ptr[j];
1382         j_end = j_beg + fc_len[j] - 1;
1383         for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
1384            fr_len[sv_ind[j_ptr]]++;
1385         /* increase total number of non-zeros */
1386         nnz += fc_len[j];
1387      }
1388      /* store total number of non-zeros */
1389      luf->nnz_f = nnz;
1390      /* check for free locations */
1391      if (luf->sv_end - luf->sv_beg < nnz)
1392      {  /* overflow of the sparse vector area */
1393         ret = 1;
1394         goto done;
1395      }
1396      /* allocate rows of the matrix F */
1397      for (i = 1; i <= n; i++)
1398      {  /* set pointer to the end of the i-th row; later this pointer
1399            will be set to the beginning of the i-th row */
1400         fr_ptr[i] = luf->sv_end;
1401         /* reserve locations for the i-th row */
1402         luf->sv_end -= fr_len[i];
1403      }
1404      /* build the matrix F in row-wise format using this matrix in
1405         column-wise format */
1406      for (j = 1; j <= n; j++)
1407      {  /* walk through elements of the j-th column */
1408         j_beg = fc_ptr[j];
1409         j_end = j_beg + fc_len[j] - 1;
1410         for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
1411         {  /* get row index */
1412            i = sv_ind[j_ptr];
1413            /* store element in the i-th row */
1414            ptr = --fr_ptr[i];
1415            sv_ind[ptr] = j;
1416            sv_val[ptr] = sv_val[j_ptr];
1417         }
1418      }
1419done: /* return to the factorizing routine */
1420      return ret;
1421}
1422
1423/***********************************************************************
1424*  NAME
1425*
1426*  luf_factorize - compute LU-factorization
1427*
1428*  SYNOPSIS
1429*
1430*  #include "glpluf.h"
1431*  int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
1432*     int ind[], double val[]), void *info);
1433*
1434*  DESCRIPTION
1435*
1436*  The routine luf_factorize computes LU-factorization of a specified
1437*  square matrix A.
1438*
1439*  The parameter luf specifies LU-factorization program object created
1440*  by the routine luf_create_it.
1441*
1442*  The parameter n specifies the order of A, n > 0.
1443*
1444*  The formal routine col specifies the matrix A to be factorized. To
1445*  obtain j-th column of A the routine luf_factorize calls the routine
1446*  col with the parameter j (1 <= j <= n). In response the routine col
1447*  should store row indices and numerical values of non-zero elements
1448*  of j-th column of A to locations ind[1,...,len] and val[1,...,len],
1449*  respectively, where len is the number of non-zeros in j-th column
1450*  returned on exit. Neither zero nor duplicate elements are allowed.
1451*
1452*  The parameter info is a transit pointer passed to the routine col.
1453*
1454*  RETURNS
1455*
1456*  0  LU-factorization has been successfully computed.
1457*
1458*  LUF_ESING
1459*     The specified matrix is singular within the working precision.
1460*     (On some elimination step the active submatrix is exactly zero,
1461*     so no pivot can be chosen.)
1462*
1463*  LUF_ECOND
1464*     The specified matrix is ill-conditioned.
1465*     (On some elimination step too intensive growth of elements of the
1466*     active submatix has been detected.)
1467*
1468*  If matrix A is well scaled, the return code LUF_ECOND may also mean
1469*  that the threshold pivoting tolerance piv_tol should be increased.
1470*
1471*  In case of non-zero return code the factorization becomes invalid.
1472*  It should not be used in other operations until the cause of failure
1473*  has been eliminated and the factorization has been recomputed again
1474*  with the routine luf_factorize.
1475*
1476*  REPAIRING SINGULAR MATRIX
1477*
1478*  If the routine luf_factorize returns non-zero code, it provides all
1479*  necessary information that can be used for "repairing" the matrix A,
1480*  where "repairing" means replacing linearly dependent columns of the
1481*  matrix A by appropriate columns of the unity matrix. This feature is
1482*  needed when this routine is used for factorizing the basis matrix
1483*  within the simplex method procedure.
1484*
1485*  On exit linearly dependent columns of the (partially transformed)
1486*  matrix U have numbers rank+1, rank+2, ..., n, where rank is estimated
1487*  rank of the matrix A stored by the routine to the member luf->rank.
1488*  The correspondence between columns of A and U is the same as between
1489*  columns of V and U. Thus, linearly dependent columns of the matrix A
1490*  have numbers qq_col[rank+1], qq_col[rank+2], ..., qq_col[n], where
1491*  qq_col is the column-like representation of the permutation matrix Q.
1492*  It is understood that each j-th linearly dependent column of the
1493*  matrix U should be replaced by the unity vector, where all elements
1494*  are zero except the unity diagonal element u[j,j]. On the other hand
1495*  j-th row of the matrix U corresponds to the row of the matrix V (and
1496*  therefore of the matrix A) with the number pp_row[j], where pp_row is
1497*  the row-like representation of the permutation matrix P. Thus, each
1498*  j-th linearly dependent column of the matrix U should be replaced by
1499*  column of the unity matrix with the number pp_row[j].
1500*
1501*  The code that repairs the matrix A may look like follows:
1502*
1503*     for (j = rank+1; j <= n; j++)
1504*     {  replace the column qq_col[j] of the matrix A by the column
1505*        pp_row[j] of the unity matrix;
1506*     }
1507*
1508*  where rank, pp_row, and qq_col are members of the structure LUF. */
1509
1510int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
1511      int ind[], double val[]), void *info)
1512{     int *pp_row, *pp_col, *qq_row, *qq_col;
1513      double max_gro = luf->max_gro;
1514      int i, j, k, p, q, t, ret;
1515      if (n < 1)
1516         xfault("luf_factorize: n = %d; invalid parameter\n", n);
1517      if (n > N_MAX)
1518         xfault("luf_factorize: n = %d; matrix too big\n", n);
1519      /* invalidate the factorization */
1520      luf->valid = 0;
1521      /* reallocate arrays, if necessary */
1522      reallocate(luf, n);
1523      pp_row = luf->pp_row;
1524      pp_col = luf->pp_col;
1525      qq_row = luf->qq_row;
1526      qq_col = luf->qq_col;
1527      /* estimate initial size of the SVA, if not specified */
1528      if (luf->sv_size == 0 && luf->new_sva == 0)
1529         luf->new_sva = 5 * (n + 10);
1530more: /* reallocate the sparse vector area, if required */
1531      if (luf->new_sva > 0)
1532      {  if (luf->sv_ind != NULL) xfree(luf->sv_ind);
1533         if (luf->sv_val != NULL) xfree(luf->sv_val);
1534         luf->sv_size = luf->new_sva;
1535         luf->sv_ind = xcalloc(1+luf->sv_size, sizeof(int));
1536         luf->sv_val = xcalloc(1+luf->sv_size, sizeof(double));
1537         luf->new_sva = 0;
1538      }
1539      /* initialize LU-factorization data structures */
1540      if (initialize(luf, col, info))
1541      {  /* overflow of the sparse vector area */
1542         luf->new_sva = luf->sv_size + luf->sv_size;
1543         xassert(luf->new_sva > luf->sv_size);
1544         goto more;
1545      }
1546      /* main elimination loop */
1547      for (k = 1; k <= n; k++)
1548      {  /* choose a pivot element v[p,q] */
1549         if (find_pivot(luf, &p, &q))
1550         {  /* no pivot can be chosen, because the active submatrix is
1551               exactly zero */
1552            luf->rank = k - 1;
1553            ret = LUF_ESING;
1554            goto done;
1555         }
1556         /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th
1557            rows and k-th and j'-th columns of the matrix U = P*V*Q to
1558            move the element u[i',j'] to the position u[k,k] */
1559         i = pp_col[p], j = qq_row[q];
1560         xassert(k <= i && i <= n && k <= j && j <= n);
1561         /* permute k-th and i-th rows of the matrix U */
1562         t = pp_row[k];
1563         pp_row[i] = t, pp_col[t] = i;
1564         pp_row[k] = p, pp_col[p] = k;
1565         /* permute k-th and j-th columns of the matrix U */
1566         t = qq_col[k];
1567         qq_col[j] = t, qq_row[t] = j;
1568         qq_col[k] = q, qq_row[q] = k;
1569         /* eliminate subdiagonal elements of k-th column of the matrix
1570            U = P*V*Q using the pivot element u[k,k] = v[p,q] */
1571         if (eliminate(luf, p, q))
1572         {  /* overflow of the sparse vector area */
1573            luf->new_sva = luf->sv_size + luf->sv_size;
1574            xassert(luf->new_sva > luf->sv_size);
1575            goto more;
1576         }
1577         /* check relative growth of elements of the matrix V */
1578         if (luf->big_v > max_gro * luf->max_a)
1579         {  /* the growth is too intensive, therefore most probably the
1580               matrix A is ill-conditioned */
1581            luf->rank = k - 1;
1582            ret = LUF_ECOND;
1583            goto done;
1584         }
1585      }
1586      /* now the matrix U = P*V*Q is upper triangular, the matrix V has
1587         been built in row-wise format, and the matrix F has been built
1588         in column-wise format */
1589      /* defragment the sparse vector area in order to merge all free
1590         locations in one continuous extent */
1591      luf_defrag_sva(luf);
1592      /* build the matrix V in column-wise format */
1593      if (build_v_cols(luf))
1594      {  /* overflow of the sparse vector area */
1595         luf->new_sva = luf->sv_size + luf->sv_size;
1596         xassert(luf->new_sva > luf->sv_size);
1597         goto more;
1598      }
1599      /* build the matrix F in row-wise format */
1600      if (build_f_rows(luf))
1601      {  /* overflow of the sparse vector area */
1602         luf->new_sva = luf->sv_size + luf->sv_size;
1603         xassert(luf->new_sva > luf->sv_size);
1604         goto more;
1605      }
1606      /* the LU-factorization has been successfully computed */
1607      luf->valid = 1;
1608      luf->rank = n;
1609      ret = 0;
1610      /* if there are few free locations in the sparse vector area, try
1611         increasing its size in the future */
1612      t = 3 * (n + luf->nnz_v) + 2 * luf->nnz_f;
1613      if (luf->sv_size < t)
1614      {  luf->new_sva = luf->sv_size;
1615         while (luf->new_sva < t)
1616         {  k = luf->new_sva;
1617            luf->new_sva = k + k;
1618            xassert(luf->new_sva > k);
1619         }
1620      }
1621done: /* return to the calling program */
1622      return ret;
1623}
1624
1625/***********************************************************************
1626*  NAME
1627*
1628*  luf_f_solve - solve system F*x = b or F'*x = b
1629*
1630*  SYNOPSIS
1631*
1632*  #include "glpluf.h"
1633*  void luf_f_solve(LUF *luf, int tr, double x[]);
1634*
1635*  DESCRIPTION
1636*
1637*  The routine luf_f_solve solves either the system F*x = b (if the
1638*  flag tr is zero) or the system F'*x = b (if the flag tr is non-zero),
1639*  where the matrix F is a component of LU-factorization specified by
1640*  the parameter luf, F' is a matrix transposed to F.
1641*
1642*  On entry the array x should contain elements of the right-hand side
1643*  vector b in locations x[1], ..., x[n], where n is the order of the
1644*  matrix F. On exit this array will contain elements of the solution
1645*  vector x in the same locations. */
1646
1647void luf_f_solve(LUF *luf, int tr, double x[])
1648{     int n = luf->n;
1649      int *fr_ptr = luf->fr_ptr;
1650      int *fr_len = luf->fr_len;
1651      int *fc_ptr = luf->fc_ptr;
1652      int *fc_len = luf->fc_len;
1653      int *pp_row = luf->pp_row;
1654      int *sv_ind = luf->sv_ind;
1655      double *sv_val = luf->sv_val;
1656      int i, j, k, beg, end, ptr;
1657      double xk;
1658      if (!luf->valid)
1659         xfault("luf_f_solve: LU-factorization is not valid\n");
1660      if (!tr)
1661      {  /* solve the system F*x = b */
1662         for (j = 1; j <= n; j++)
1663         {  k = pp_row[j];
1664            xk = x[k];
1665            if (xk != 0.0)
1666            {  beg = fc_ptr[k];
1667               end = beg + fc_len[k] - 1;
1668               for (ptr = beg; ptr <= end; ptr++)
1669                  x[sv_ind[ptr]] -= sv_val[ptr] * xk;
1670            }
1671         }
1672      }
1673      else
1674      {  /* solve the system F'*x = b */
1675         for (i = n; i >= 1; i--)
1676         {  k = pp_row[i];
1677            xk = x[k];
1678            if (xk != 0.0)
1679            {  beg = fr_ptr[k];
1680               end = beg + fr_len[k] - 1;
1681               for (ptr = beg; ptr <= end; ptr++)
1682                  x[sv_ind[ptr]] -= sv_val[ptr] * xk;
1683            }
1684         }
1685      }
1686      return;
1687}
1688
1689/***********************************************************************
1690*  NAME
1691*
1692*  luf_v_solve - solve system V*x = b or V'*x = b
1693*
1694*  SYNOPSIS
1695*
1696*  #include "glpluf.h"
1697*  void luf_v_solve(LUF *luf, int tr, double x[]);
1698*
1699*  DESCRIPTION
1700*
1701*  The routine luf_v_solve solves either the system V*x = b (if the
1702*  flag tr is zero) or the system V'*x = b (if the flag tr is non-zero),
1703*  where the matrix V is a component of LU-factorization specified by
1704*  the parameter luf, V' is a matrix transposed to V.
1705*
1706*  On entry the array x should contain elements of the right-hand side
1707*  vector b in locations x[1], ..., x[n], where n is the order of the
1708*  matrix V. On exit this array will contain elements of the solution
1709*  vector x in the same locations. */
1710
1711void luf_v_solve(LUF *luf, int tr, double x[])
1712{     int n = luf->n;
1713      int *vr_ptr = luf->vr_ptr;
1714      int *vr_len = luf->vr_len;
1715      double *vr_piv = luf->vr_piv;
1716      int *vc_ptr = luf->vc_ptr;
1717      int *vc_len = luf->vc_len;
1718      int *pp_row = luf->pp_row;
1719      int *qq_col = luf->qq_col;
1720      int *sv_ind = luf->sv_ind;
1721      double *sv_val = luf->sv_val;
1722      double *b = luf->work;
1723      int i, j, k, beg, end, ptr;
1724      double temp;
1725      if (!luf->valid)
1726         xfault("luf_v_solve: LU-factorization is not valid\n");
1727      for (k = 1; k <= n; k++) b[k] = x[k], x[k] = 0.0;
1728      if (!tr)
1729      {  /* solve the system V*x = b */
1730         for (k = n; k >= 1; k--)
1731         {  i = pp_row[k], j = qq_col[k];
1732            temp = b[i];
1733            if (temp != 0.0)
1734            {  x[j] = (temp /= vr_piv[i]);
1735               beg = vc_ptr[j];
1736               end = beg + vc_len[j] - 1;
1737               for (ptr = beg; ptr <= end; ptr++)
1738                  b[sv_ind[ptr]] -= sv_val[ptr] * temp;
1739            }
1740         }
1741      }
1742      else
1743      {  /* solve the system V'*x = b */
1744         for (k = 1; k <= n; k++)
1745         {  i = pp_row[k], j = qq_col[k];
1746            temp = b[j];
1747            if (temp != 0.0)
1748            {  x[i] = (temp /= vr_piv[i]);
1749               beg = vr_ptr[i];
1750               end = beg + vr_len[i] - 1;
1751               for (ptr = beg; ptr <= end; ptr++)
1752                  b[sv_ind[ptr]] -= sv_val[ptr] * temp;
1753            }
1754         }
1755      }
1756      return;
1757}
1758
1759/***********************************************************************
1760*  NAME
1761*
1762*  luf_a_solve - solve system A*x = b or A'*x = b
1763*
1764*  SYNOPSIS
1765*
1766*  #include "glpluf.h"
1767*  void luf_a_solve(LUF *luf, int tr, double x[]);
1768*
1769*  DESCRIPTION
1770*
1771*  The routine luf_a_solve solves either the system A*x = b (if the
1772*  flag tr is zero) or the system A'*x = b (if the flag tr is non-zero),
1773*  where the parameter luf specifies LU-factorization of the matrix A,
1774*  A' is a matrix transposed to A.
1775*
1776*  On entry the array x should contain elements of the right-hand side
1777*  vector b in locations x[1], ..., x[n], where n is the order of the
1778*  matrix A. On exit this array will contain elements of the solution
1779*  vector x in the same locations. */
1780
1781void luf_a_solve(LUF *luf, int tr, double x[])
1782{     if (!luf->valid)
1783         xfault("luf_a_solve: LU-factorization is not valid\n");
1784      if (!tr)
1785      {  /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
1786         luf_f_solve(luf, 0, x);
1787         luf_v_solve(luf, 0, x);
1788      }
1789      else
1790      {  /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
1791         luf_v_solve(luf, 1, x);
1792         luf_f_solve(luf, 1, x);
1793      }
1794      return;
1795}
1796
1797/***********************************************************************
1798*  NAME
1799*
1800*  luf_delete_it - delete LU-factorization
1801*
1802*  SYNOPSIS
1803*
1804*  #include "glpluf.h"
1805*  void luf_delete_it(LUF *luf);
1806*
1807*  DESCRIPTION
1808*
1809*  The routine luf_delete deletes LU-factorization specified by the
1810*  parameter luf and frees all the memory allocated to this program
1811*  object. */
1812
1813void luf_delete_it(LUF *luf)
1814{     if (luf->fr_ptr != NULL) xfree(luf->fr_ptr);
1815      if (luf->fr_len != NULL) xfree(luf->fr_len);
1816      if (luf->fc_ptr != NULL) xfree(luf->fc_ptr);
1817      if (luf->fc_len != NULL) xfree(luf->fc_len);
1818      if (luf->vr_ptr != NULL) xfree(luf->vr_ptr);
1819      if (luf->vr_len != NULL) xfree(luf->vr_len);
1820      if (luf->vr_cap != NULL) xfree(luf->vr_cap);
1821      if (luf->vr_piv != NULL) xfree(luf->vr_piv);
1822      if (luf->vc_ptr != NULL) xfree(luf->vc_ptr);
1823      if (luf->vc_len != NULL) xfree(luf->vc_len);
1824      if (luf->vc_cap != NULL) xfree(luf->vc_cap);
1825      if (luf->pp_row != NULL) xfree(luf->pp_row);
1826      if (luf->pp_col != NULL) xfree(luf->pp_col);
1827      if (luf->qq_row != NULL) xfree(luf->qq_row);
1828      if (luf->qq_col != NULL) xfree(luf->qq_col);
1829      if (luf->sv_ind != NULL) xfree(luf->sv_ind);
1830      if (luf->sv_val != NULL) xfree(luf->sv_val);
1831      if (luf->sv_prev != NULL) xfree(luf->sv_prev);
1832      if (luf->sv_next != NULL) xfree(luf->sv_next);
1833      if (luf->vr_max != NULL) xfree(luf->vr_max);
1834      if (luf->rs_head != NULL) xfree(luf->rs_head);
1835      if (luf->rs_prev != NULL) xfree(luf->rs_prev);
1836      if (luf->rs_next != NULL) xfree(luf->rs_next);
1837      if (luf->cs_head != NULL) xfree(luf->cs_head);
1838      if (luf->cs_prev != NULL) xfree(luf->cs_prev);
1839      if (luf->cs_next != NULL) xfree(luf->cs_next);
1840      if (luf->flag != NULL) xfree(luf->flag);
1841      if (luf->work != NULL) xfree(luf->work);
1842      xfree(luf);
1843      return;
1844}
1845
1846/* eof */
Note: See TracBrowser for help on using the repository browser.