COIN-OR::LEMON - Graph Library

source: lemon-project-template-glpk/deps/glpk/src/glpmat.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: 32.2 KB
Line 
1/* glpmat.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 "glpenv.h"
26#include "glpmat.h"
27#include "glpqmd.h"
28#include "amd/amd.h"
29#include "colamd/colamd.h"
30
31/*----------------------------------------------------------------------
32-- check_fvs - check sparse vector in full-vector storage format.
33--
34-- SYNOPSIS
35--
36-- #include "glpmat.h"
37-- int check_fvs(int n, int nnz, int ind[], double vec[]);
38--
39-- DESCRIPTION
40--
41-- The routine check_fvs checks if a given vector of dimension n in
42-- full-vector storage format has correct representation.
43--
44-- RETURNS
45--
46-- The routine returns one of the following codes:
47--
48-- 0 - the vector is correct;
49-- 1 - the number of elements (n) is negative;
50-- 2 - the number of non-zero elements (nnz) is negative;
51-- 3 - some element index is out of range;
52-- 4 - some element index is duplicate;
53-- 5 - some non-zero element is out of pattern. */
54
55int check_fvs(int n, int nnz, int ind[], double vec[])
56{     int i, t, ret, *flag = NULL;
57      /* check the number of elements */
58      if (n < 0)
59      {  ret = 1;
60         goto done;
61      }
62      /* check the number of non-zero elements */
63      if (nnz < 0)
64      {  ret = 2;
65         goto done;
66      }
67      /* check vector indices */
68      flag = xcalloc(1+n, sizeof(int));
69      for (i = 1; i <= n; i++) flag[i] = 0;
70      for (t = 1; t <= nnz; t++)
71      {  i = ind[t];
72         if (!(1 <= i && i <= n))
73         {  ret = 3;
74            goto done;
75         }
76         if (flag[i])
77         {  ret = 4;
78            goto done;
79         }
80         flag[i] = 1;
81      }
82      /* check vector elements */
83      for (i = 1; i <= n; i++)
84      {  if (!flag[i] && vec[i] != 0.0)
85         {  ret = 5;
86            goto done;
87         }
88      }
89      /* the vector is ok */
90      ret = 0;
91done: if (flag != NULL) xfree(flag);
92      return ret;
93}
94
95/*----------------------------------------------------------------------
96-- check_pattern - check pattern of sparse matrix.
97--
98-- SYNOPSIS
99--
100-- #include "glpmat.h"
101-- int check_pattern(int m, int n, int A_ptr[], int A_ind[]);
102--
103-- DESCRIPTION
104--
105-- The routine check_pattern checks the pattern of a given mxn matrix
106-- in storage-by-rows format.
107--
108-- RETURNS
109--
110-- The routine returns one of the following codes:
111--
112-- 0 - the pattern is correct;
113-- 1 - the number of rows (m) is negative;
114-- 2 - the number of columns (n) is negative;
115-- 3 - A_ptr[1] is not 1;
116-- 4 - some column index is out of range;
117-- 5 - some column indices are duplicate. */
118
119int check_pattern(int m, int n, int A_ptr[], int A_ind[])
120{     int i, j, ptr, ret, *flag = NULL;
121      /* check the number of rows */
122      if (m < 0)
123      {  ret = 1;
124         goto done;
125      }
126      /* check the number of columns */
127      if (n < 0)
128      {  ret = 2;
129         goto done;
130      }
131      /* check location A_ptr[1] */
132      if (A_ptr[1] != 1)
133      {  ret = 3;
134         goto done;
135      }
136      /* check row patterns */
137      flag = xcalloc(1+n, sizeof(int));
138      for (j = 1; j <= n; j++) flag[j] = 0;
139      for (i = 1; i <= m; i++)
140      {  /* check pattern of row i */
141         for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
142         {  j = A_ind[ptr];
143            /* check column index */
144            if (!(1 <= j && j <= n))
145            {  ret = 4;
146               goto done;
147            }
148            /* check for duplication */
149            if (flag[j])
150            {  ret = 5;
151               goto done;
152            }
153            flag[j] = 1;
154         }
155         /* clear flags */
156         for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
157         {  j = A_ind[ptr];
158            flag[j] = 0;
159         }
160      }
161      /* the pattern is ok */
162      ret = 0;
163done: if (flag != NULL) xfree(flag);
164      return ret;
165}
166
167/*----------------------------------------------------------------------
168-- transpose - transpose sparse matrix.
169--
170-- *Synopsis*
171--
172-- #include "glpmat.h"
173-- void transpose(int m, int n, int A_ptr[], int A_ind[],
174--    double A_val[], int AT_ptr[], int AT_ind[], double AT_val[]);
175--
176-- *Description*
177--
178-- For a given mxn sparse matrix A the routine transpose builds a nxm
179-- sparse matrix A' which is a matrix transposed to A.
180--
181-- The arrays A_ptr, A_ind, and A_val specify a given mxn matrix A to
182-- be transposed in storage-by-rows format. The parameter A_val can be
183-- NULL, in which case numeric values are not copied. The arrays A_ptr,
184-- A_ind, and A_val are not changed on exit.
185--
186-- On entry the arrays AT_ptr, AT_ind, and AT_val must be allocated,
187-- but their content is ignored. On exit the routine stores a resultant
188-- nxm matrix A' in these arrays in storage-by-rows format. Note that
189-- if the parameter A_val is NULL, the array AT_val is not used.
190--
191-- The routine transpose has a side effect that elements in rows of the
192-- resultant matrix A' follow in ascending their column indices. */
193
194void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[],
195      int AT_ptr[], int AT_ind[], double AT_val[])
196{     int i, j, t, beg, end, pos, len;
197      /* determine row lengths of resultant matrix */
198      for (j = 1; j <= n; j++) AT_ptr[j] = 0;
199      for (i = 1; i <= m; i++)
200      {  beg = A_ptr[i], end = A_ptr[i+1];
201         for (t = beg; t < end; t++) AT_ptr[A_ind[t]]++;
202      }
203      /* set up row pointers of resultant matrix */
204      pos = 1;
205      for (j = 1; j <= n; j++)
206         len = AT_ptr[j], pos += len, AT_ptr[j] = pos;
207      AT_ptr[n+1] = pos;
208      /* build resultant matrix */
209      for (i = m; i >= 1; i--)
210      {  beg = A_ptr[i], end = A_ptr[i+1];
211         for (t = beg; t < end; t++)
212         {  pos = --AT_ptr[A_ind[t]];
213            AT_ind[pos] = i;
214            if (A_val != NULL) AT_val[pos] = A_val[t];
215         }
216      }
217      return;
218}
219
220/*----------------------------------------------------------------------
221-- adat_symbolic - compute S = P*A*D*A'*P' (symbolic phase).
222--
223-- *Synopsis*
224--
225-- #include "glpmat.h"
226-- int *adat_symbolic(int m, int n, int P_per[], int A_ptr[],
227--    int A_ind[], int S_ptr[]);
228--
229-- *Description*
230--
231-- The routine adat_symbolic implements the symbolic phase to compute
232-- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
233-- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
234-- transposed to A, P' is an inverse of P.
235--
236-- The parameter m is the number of rows in A and the order of P.
237--
238-- The parameter n is the number of columns in A and the order of D.
239--
240-- The array P_per specifies permutation matrix P. It is not changed on
241-- exit.
242--
243-- The arrays A_ptr and A_ind specify the pattern of matrix A. They are
244-- not changed on exit.
245--
246-- On exit the routine stores the pattern of upper triangular part of
247-- matrix S without diagonal elements in the arrays S_ptr and S_ind in
248-- storage-by-rows format. The array S_ptr should be allocated on entry,
249-- however, its content is ignored. The array S_ind is allocated by the
250-- routine itself which returns a pointer to it.
251--
252-- *Returns*
253--
254-- The routine returns a pointer to the array S_ind. */
255
256int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[],
257      int S_ptr[])
258{     int i, j, t, ii, jj, tt, k, size, len;
259      int *S_ind, *AT_ptr, *AT_ind, *ind, *map, *temp;
260      /* build the pattern of A', which is a matrix transposed to A, to
261         efficiently access A in column-wise manner */
262      AT_ptr = xcalloc(1+n+1, sizeof(int));
263      AT_ind = xcalloc(A_ptr[m+1], sizeof(int));
264      transpose(m, n, A_ptr, A_ind, NULL, AT_ptr, AT_ind, NULL);
265      /* allocate the array S_ind */
266      size = A_ptr[m+1] - 1;
267      if (size < m) size = m;
268      S_ind = xcalloc(1+size, sizeof(int));
269      /* allocate and initialize working arrays */
270      ind = xcalloc(1+m, sizeof(int));
271      map = xcalloc(1+m, sizeof(int));
272      for (jj = 1; jj <= m; jj++) map[jj] = 0;
273      /* compute pattern of S; note that symbolically S = B*B', where
274         B = P*A, B' is matrix transposed to B */
275      S_ptr[1] = 1;
276      for (ii = 1; ii <= m; ii++)
277      {  /* compute pattern of ii-th row of S */
278         len = 0;
279         i = P_per[ii]; /* i-th row of A = ii-th row of B */
280         for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
281         {  k = A_ind[t];
282            /* walk through k-th column of A */
283            for (tt = AT_ptr[k]; tt < AT_ptr[k+1]; tt++)
284            {  j = AT_ind[tt];
285               jj = P_per[m+j]; /* j-th row of A = jj-th row of B */
286               /* a[i,k] != 0 and a[j,k] != 0 ergo s[ii,jj] != 0 */
287               if (ii < jj && !map[jj]) ind[++len] = jj, map[jj] = 1;
288            }
289         }
290         /* now (ind) is pattern of ii-th row of S */
291         S_ptr[ii+1] = S_ptr[ii] + len;
292         /* at least (S_ptr[ii+1] - 1) locations should be available in
293            the array S_ind */
294         if (S_ptr[ii+1] - 1 > size)
295         {  temp = S_ind;
296            size += size;
297            S_ind = xcalloc(1+size, sizeof(int));
298            memcpy(&S_ind[1], &temp[1], (S_ptr[ii] - 1) * sizeof(int));
299            xfree(temp);
300         }
301         xassert(S_ptr[ii+1] - 1 <= size);
302         /* (ii-th row of S) := (ind) */
303         memcpy(&S_ind[S_ptr[ii]], &ind[1], len * sizeof(int));
304         /* clear the row pattern map */
305         for (t = 1; t <= len; t++) map[ind[t]] = 0;
306      }
307      /* free working arrays */
308      xfree(AT_ptr);
309      xfree(AT_ind);
310      xfree(ind);
311      xfree(map);
312      /* reallocate the array S_ind to free unused locations */
313      temp = S_ind;
314      size = S_ptr[m+1] - 1;
315      S_ind = xcalloc(1+size, sizeof(int));
316      memcpy(&S_ind[1], &temp[1], size * sizeof(int));
317      xfree(temp);
318      return S_ind;
319}
320
321/*----------------------------------------------------------------------
322-- adat_numeric - compute S = P*A*D*A'*P' (numeric phase).
323--
324-- *Synopsis*
325--
326-- #include "glpmat.h"
327-- void adat_numeric(int m, int n, int P_per[],
328--    int A_ptr[], int A_ind[], double A_val[], double D_diag[],
329--    int S_ptr[], int S_ind[], double S_val[], double S_diag[]);
330--
331-- *Description*
332--
333-- The routine adat_numeric implements the numeric phase to compute
334-- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
335-- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
336-- transposed to A, P' is an inverse of P.
337--
338-- The parameter m is the number of rows in A and the order of P.
339--
340-- The parameter n is the number of columns in A and the order of D.
341--
342-- The matrix P is specified in the array P_per, which is not changed
343-- on exit.
344--
345-- The matrix A is specified in the arrays A_ptr, A_ind, and A_val in
346-- storage-by-rows format. These arrays are not changed on exit.
347--
348-- Diagonal elements of the matrix D are specified in the array D_diag,
349-- where D_diag[0] is not used, D_diag[i] = d[i,i] for i = 1, ..., n.
350-- The array D_diag is not changed on exit.
351--
352-- The pattern of the upper triangular part of the matrix S without
353-- diagonal elements (previously computed by the routine adat_symbolic)
354-- is specified in the arrays S_ptr and S_ind, which are not changed on
355-- exit. Numeric values of non-diagonal elements of S are stored in
356-- corresponding locations of the array S_val, and values of diagonal
357-- elements of S are stored in locations S_diag[1], ..., S_diag[n]. */
358
359void adat_numeric(int m, int n, int P_per[],
360      int A_ptr[], int A_ind[], double A_val[], double D_diag[],
361      int S_ptr[], int S_ind[], double S_val[], double S_diag[])
362{     int i, j, t, ii, jj, tt, beg, end, beg1, end1, k;
363      double sum, *work;
364      work = xcalloc(1+n, sizeof(double));
365      for (j = 1; j <= n; j++) work[j] = 0.0;
366      /* compute S = B*D*B', where B = P*A, B' is a matrix transposed
367         to B */
368      for (ii = 1; ii <= m; ii++)
369      {  i = P_per[ii]; /* i-th row of A = ii-th row of B */
370         /* (work) := (i-th row of A) */
371         beg = A_ptr[i], end = A_ptr[i+1];
372         for (t = beg; t < end; t++)
373            work[A_ind[t]] = A_val[t];
374         /* compute ii-th row of S */
375         beg = S_ptr[ii], end = S_ptr[ii+1];
376         for (t = beg; t < end; t++)
377         {  jj = S_ind[t];
378            j = P_per[jj]; /* j-th row of A = jj-th row of B */
379            /* s[ii,jj] := sum a[i,k] * d[k,k] * a[j,k] */
380            sum = 0.0;
381            beg1 = A_ptr[j], end1 = A_ptr[j+1];
382            for (tt = beg1; tt < end1; tt++)
383            {  k = A_ind[tt];
384               sum += work[k] * D_diag[k] * A_val[tt];
385            }
386            S_val[t] = sum;
387         }
388         /* s[ii,ii] := sum a[i,k] * d[k,k] * a[i,k] */
389         sum = 0.0;
390         beg = A_ptr[i], end = A_ptr[i+1];
391         for (t = beg; t < end; t++)
392         {  k = A_ind[t];
393            sum += A_val[t] * D_diag[k] * A_val[t];
394            work[k] = 0.0;
395         }
396         S_diag[ii] = sum;
397      }
398      xfree(work);
399      return;
400}
401
402/*----------------------------------------------------------------------
403-- min_degree - minimum degree ordering.
404--
405-- *Synopsis*
406--
407-- #include "glpmat.h"
408-- void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]);
409--
410-- *Description*
411--
412-- The routine min_degree uses the minimum degree ordering algorithm
413-- to find a permutation matrix P for a given sparse symmetric positive
414-- matrix A which minimizes the number of non-zeros in upper triangular
415-- factor U for Cholesky factorization P*A*P' = U'*U.
416--
417-- The parameter n is the order of matrices A and P.
418--
419-- The pattern of the given matrix A is specified on entry in the arrays
420-- A_ptr and A_ind in storage-by-rows format. Only the upper triangular
421-- part without diagonal elements (which all are assumed to be non-zero)
422-- should be specified as if A were upper triangular. The arrays A_ptr
423-- and A_ind are not changed on exit.
424--
425-- The permutation matrix P is stored by the routine in the array P_per
426-- on exit.
427--
428-- *Algorithm*
429--
430-- The routine min_degree is based on some subroutines from the package
431-- SPARSPAK (see comments in the module glpqmd). */
432
433void min_degree(int n, int A_ptr[], int A_ind[], int P_per[])
434{     int i, j, ne, t, pos, len;
435      int *xadj, *adjncy, *deg, *marker, *rchset, *nbrhd, *qsize,
436         *qlink, nofsub;
437      /* determine number of non-zeros in complete pattern */
438      ne = A_ptr[n+1] - 1;
439      ne += ne;
440      /* allocate working arrays */
441      xadj = xcalloc(1+n+1, sizeof(int));
442      adjncy = xcalloc(1+ne, sizeof(int));
443      deg = xcalloc(1+n, sizeof(int));
444      marker = xcalloc(1+n, sizeof(int));
445      rchset = xcalloc(1+n, sizeof(int));
446      nbrhd = xcalloc(1+n, sizeof(int));
447      qsize = xcalloc(1+n, sizeof(int));
448      qlink = xcalloc(1+n, sizeof(int));
449      /* determine row lengths in complete pattern */
450      for (i = 1; i <= n; i++) xadj[i] = 0;
451      for (i = 1; i <= n; i++)
452      {  for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
453         {  j = A_ind[t];
454            xassert(i < j && j <= n);
455            xadj[i]++, xadj[j]++;
456         }
457      }
458      /* set up row pointers for complete pattern */
459      pos = 1;
460      for (i = 1; i <= n; i++)
461         len = xadj[i], pos += len, xadj[i] = pos;
462      xadj[n+1] = pos;
463      xassert(pos - 1 == ne);
464      /* construct complete pattern */
465      for (i = 1; i <= n; i++)
466      {  for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
467         {  j = A_ind[t];
468            adjncy[--xadj[i]] = j, adjncy[--xadj[j]] = i;
469         }
470      }
471      /* call the main minimimum degree ordering routine */
472      genqmd(&n, xadj, adjncy, P_per, P_per + n, deg, marker, rchset,
473         nbrhd, qsize, qlink, &nofsub);
474      /* make sure that permutation matrix P is correct */
475      for (i = 1; i <= n; i++)
476      {  j = P_per[i];
477         xassert(1 <= j && j <= n);
478         xassert(P_per[n+j] == i);
479      }
480      /* free working arrays */
481      xfree(xadj);
482      xfree(adjncy);
483      xfree(deg);
484      xfree(marker);
485      xfree(rchset);
486      xfree(nbrhd);
487      xfree(qsize);
488      xfree(qlink);
489      return;
490}
491
492/**********************************************************************/
493
494void amd_order1(int n, int A_ptr[], int A_ind[], int P_per[])
495{     /* approximate minimum degree ordering (AMD) */
496      int k, ret;
497      double Control[AMD_CONTROL], Info[AMD_INFO];
498      /* get the default parameters */
499      amd_defaults(Control);
500#if 0
501      /* and print them */
502      amd_control(Control);
503#endif
504      /* make all indices 0-based */
505      for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--;
506      for (k = 1; k <= n+1; k++) A_ptr[k]--;
507      /* call the ordering routine */
508      ret = amd_order(n, &A_ptr[1], &A_ind[1], &P_per[1], Control, Info)
509         ;
510#if 0
511      amd_info(Info);
512#endif
513      xassert(ret == AMD_OK || ret == AMD_OK_BUT_JUMBLED);
514      /* retsore 1-based indices */
515      for (k = 1; k <= n+1; k++) A_ptr[k]++;
516      for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++;
517      /* patch up permutation matrix */
518      memset(&P_per[n+1], 0, n * sizeof(int));
519      for (k = 1; k <= n; k++)
520      {  P_per[k]++;
521         xassert(1 <= P_per[k] && P_per[k] <= n);
522         xassert(P_per[n+P_per[k]] == 0);
523         P_per[n+P_per[k]] = k;
524      }
525      return;
526}
527
528/**********************************************************************/
529
530static void *allocate(size_t n, size_t size)
531{     void *ptr;
532      ptr = xcalloc(n, size);
533      memset(ptr, 0, n * size);
534      return ptr;
535}
536
537static void release(void *ptr)
538{     xfree(ptr);
539      return;
540}
541
542void symamd_ord(int n, int A_ptr[], int A_ind[], int P_per[])
543{     /* approximate minimum degree ordering (SYMAMD) */
544      int k, ok;
545      int stats[COLAMD_STATS];
546      /* make all indices 0-based */
547      for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]--;
548      for (k = 1; k <= n+1; k++) A_ptr[k]--;
549      /* call the ordering routine */
550      ok = symamd(n, &A_ind[1], &A_ptr[1], &P_per[1], NULL, stats,
551         allocate, release);
552#if 0
553      symamd_report(stats);
554#endif
555      xassert(ok);
556      /* restore 1-based indices */
557      for (k = 1; k <= n+1; k++) A_ptr[k]++;
558      for (k = 1; k < A_ptr[n+1]; k++) A_ind[k]++;
559      /* patch up permutation matrix */
560      memset(&P_per[n+1], 0, n * sizeof(int));
561      for (k = 1; k <= n; k++)
562      {  P_per[k]++;
563         xassert(1 <= P_per[k] && P_per[k] <= n);
564         xassert(P_per[n+P_per[k]] == 0);
565         P_per[n+P_per[k]] = k;
566      }
567      return;
568}
569
570/*----------------------------------------------------------------------
571-- chol_symbolic - compute Cholesky factorization (symbolic phase).
572--
573-- *Synopsis*
574--
575-- #include "glpmat.h"
576-- int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]);
577--
578-- *Description*
579--
580-- The routine chol_symbolic implements the symbolic phase of Cholesky
581-- factorization A = U'*U, where A is a given sparse symmetric positive
582-- definite matrix, U is a resultant upper triangular factor, U' is a
583-- matrix transposed to U.
584--
585-- The parameter n is the order of matrices A and U.
586--
587-- The pattern of the given matrix A is specified on entry in the arrays
588-- A_ptr and A_ind in storage-by-rows format. Only the upper triangular
589-- part without diagonal elements (which all are assumed to be non-zero)
590-- should be specified as if A were upper triangular. The arrays A_ptr
591-- and A_ind are not changed on exit.
592--
593-- The pattern of the matrix U without diagonal elements (which all are
594-- assumed to be non-zero) is stored on exit from the routine in the
595-- arrays U_ptr and U_ind in storage-by-rows format. The array U_ptr
596-- should be allocated on entry, however, its content is ignored. The
597-- array U_ind is allocated by the routine which returns a pointer to it
598-- on exit.
599--
600-- *Returns*
601--
602-- The routine returns a pointer to the array U_ind.
603--
604-- *Method*
605--
606-- The routine chol_symbolic computes the pattern of the matrix U in a
607-- row-wise manner. No pivoting is used.
608--
609-- It is known that to compute the pattern of row k of the matrix U we
610-- need to merge the pattern of row k of the matrix A and the patterns
611-- of each row i of U, where u[i,k] is non-zero (these rows are already
612-- computed and placed above row k).
613--
614-- However, to reduce the number of rows to be merged the routine uses
615-- an advanced algorithm proposed in:
616--
617-- D.J.Rose, R.E.Tarjan, and G.S.Lueker. Algorithmic aspects of vertex
618-- elimination on graphs. SIAM J. Comput. 5, 1976, 266-83.
619--
620-- The authors of the cited paper show that we have the same result if
621-- we merge row k of the matrix A and such rows of the matrix U (among
622-- rows 1, ..., k-1) whose leftmost non-diagonal non-zero element is
623-- placed in k-th column. This feature signficantly reduces the number
624-- of rows to be merged, especially on the final steps, where rows of
625-- the matrix U become quite dense.
626--
627-- To determine rows, which should be merged on k-th step, for a fixed
628-- time the routine uses linked lists of row numbers of the matrix U.
629-- Location head[k] contains the number of a first row, whose leftmost
630-- non-diagonal non-zero element is placed in column k, and location
631-- next[i] contains the number of a next row with the same property as
632-- row i. */
633
634int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[])
635{     int i, j, k, t, len, size, beg, end, min_j, *U_ind, *head, *next,
636         *ind, *map, *temp;
637      /* initially we assume that on computing the pattern of U fill-in
638         will double the number of non-zeros in A */
639      size = A_ptr[n+1] - 1;
640      if (size < n) size = n;
641      size += size;
642      U_ind = xcalloc(1+size, sizeof(int));
643      /* allocate and initialize working arrays */
644      head = xcalloc(1+n, sizeof(int));
645      for (i = 1; i <= n; i++) head[i] = 0;
646      next = xcalloc(1+n, sizeof(int));
647      ind = xcalloc(1+n, sizeof(int));
648      map = xcalloc(1+n, sizeof(int));
649      for (j = 1; j <= n; j++) map[j] = 0;
650      /* compute the pattern of matrix U */
651      U_ptr[1] = 1;
652      for (k = 1; k <= n; k++)
653      {  /* compute the pattern of k-th row of U, which is the union of
654            k-th row of A and those rows of U (among 1, ..., k-1) whose
655            leftmost non-diagonal non-zero is placed in k-th column */
656         /* (ind) := (k-th row of A) */
657         len = A_ptr[k+1] - A_ptr[k];
658         memcpy(&ind[1], &A_ind[A_ptr[k]], len * sizeof(int));
659         for (t = 1; t <= len; t++)
660         {  j = ind[t];
661            xassert(k < j && j <= n);
662            map[j] = 1;
663         }
664         /* walk through rows of U whose leftmost non-diagonal non-zero
665            is placed in k-th column */
666         for (i = head[k]; i != 0; i = next[i])
667         {  /* (ind) := (ind) union (i-th row of U) */
668            beg = U_ptr[i], end = U_ptr[i+1];
669            for (t = beg; t < end; t++)
670            {  j = U_ind[t];
671               if (j > k && !map[j]) ind[++len] = j, map[j] = 1;
672            }
673         }
674         /* now (ind) is the pattern of k-th row of U */
675         U_ptr[k+1] = U_ptr[k] + len;
676         /* at least (U_ptr[k+1] - 1) locations should be available in
677            the array U_ind */
678         if (U_ptr[k+1] - 1 > size)
679         {  temp = U_ind;
680            size += size;
681            U_ind = xcalloc(1+size, sizeof(int));
682            memcpy(&U_ind[1], &temp[1], (U_ptr[k] - 1) * sizeof(int));
683            xfree(temp);
684         }
685         xassert(U_ptr[k+1] - 1 <= size);
686         /* (k-th row of U) := (ind) */
687         memcpy(&U_ind[U_ptr[k]], &ind[1], len * sizeof(int));
688         /* determine column index of leftmost non-diagonal non-zero in
689            k-th row of U and clear the row pattern map */
690         min_j = n + 1;
691         for (t = 1; t <= len; t++)
692         {  j = ind[t], map[j] = 0;
693            if (min_j > j) min_j = j;
694         }
695         /* include k-th row into corresponding linked list */
696         if (min_j <= n) next[k] = head[min_j], head[min_j] = k;
697      }
698      /* free working arrays */
699      xfree(head);
700      xfree(next);
701      xfree(ind);
702      xfree(map);
703      /* reallocate the array U_ind to free unused locations */
704      temp = U_ind;
705      size = U_ptr[n+1] - 1;
706      U_ind = xcalloc(1+size, sizeof(int));
707      memcpy(&U_ind[1], &temp[1], size * sizeof(int));
708      xfree(temp);
709      return U_ind;
710}
711
712/*----------------------------------------------------------------------
713-- chol_numeric - compute Cholesky factorization (numeric phase).
714--
715-- *Synopsis*
716--
717-- #include "glpmat.h"
718-- int chol_numeric(int n,
719--    int A_ptr[], int A_ind[], double A_val[], double A_diag[],
720--    int U_ptr[], int U_ind[], double U_val[], double U_diag[]);
721--
722-- *Description*
723--
724-- The routine chol_symbolic implements the numeric phase of Cholesky
725-- factorization A = U'*U, where A is a given sparse symmetric positive
726-- definite matrix, U is a resultant upper triangular factor, U' is a
727-- matrix transposed to U.
728--
729-- The parameter n is the order of matrices A and U.
730--
731-- Upper triangular part of the matrix A without diagonal elements is
732-- specified in the arrays A_ptr, A_ind, and A_val in storage-by-rows
733-- format. Diagonal elements of A are specified in the array A_diag,
734-- where A_diag[0] is not used, A_diag[i] = a[i,i] for i = 1, ..., n.
735-- The arrays A_ptr, A_ind, A_val, and A_diag are not changed on exit.
736--
737-- The pattern of the matrix U without diagonal elements (previously
738-- computed with the routine chol_symbolic) is specified in the arrays
739-- U_ptr and U_ind, which are not changed on exit. Numeric values of
740-- non-diagonal elements of U are stored in corresponding locations of
741-- the array U_val, and values of diagonal elements of U are stored in
742-- locations U_diag[1], ..., U_diag[n].
743--
744-- *Returns*
745--
746-- The routine returns the number of non-positive diagonal elements of
747-- the matrix U which have been replaced by a huge positive number (see
748-- the method description below). Zero return code means the matrix A
749-- has been successfully factorized.
750--
751-- *Method*
752--
753-- The routine chol_numeric computes the matrix U in a row-wise manner
754-- using standard gaussian elimination technique. No pivoting is used.
755--
756-- Initially the routine sets U = A, and before k-th elimination step
757-- the matrix U is the following:
758--
759--       1       k         n
760--    1  x x x x x x x x x x
761--       . x x x x x x x x x
762--       . . x x x x x x x x
763--       . . . x x x x x x x
764--    k  . . . . * * * * * *
765--       . . . . * * * * * *
766--       . . . . * * * * * *
767--       . . . . * * * * * *
768--       . . . . * * * * * *
769--    n  . . . . * * * * * *
770--
771-- where 'x' are elements of already computed rows, '*' are elements of
772-- the active submatrix. (Note that the lower triangular part of the
773-- active submatrix being symmetric is not stored and diagonal elements
774-- are stored separately in the array U_diag.)
775--
776-- The matrix A is assumed to be positive definite. However, if it is
777-- close to semi-definite, on some elimination step a pivot u[k,k] may
778-- happen to be non-positive due to round-off errors. In this case the
779-- routine uses a technique proposed in:
780--
781-- S.J.Wright. The Cholesky factorization in interior-point and barrier
782-- methods. Preprint MCS-P600-0596, Mathematics and Computer Science
783-- Division, Argonne National Laboratory, Argonne, Ill., May 1996.
784--
785-- The routine just replaces non-positive u[k,k] by a huge positive
786-- number. This involves non-diagonal elements in k-th row of U to be
787-- close to zero that, in turn, involves k-th component of a solution
788-- vector to be close to zero. Note, however, that this technique works
789-- only if the system A*x = b is consistent. */
790
791int chol_numeric(int n,
792      int A_ptr[], int A_ind[], double A_val[], double A_diag[],
793      int U_ptr[], int U_ind[], double U_val[], double U_diag[])
794{     int i, j, k, t, t1, beg, end, beg1, end1, count = 0;
795      double ukk, uki, *work;
796      work = xcalloc(1+n, sizeof(double));
797      for (j = 1; j <= n; j++) work[j] = 0.0;
798      /* U := (upper triangle of A) */
799      /* note that the upper traingle of A is a subset of U */
800      for (i = 1; i <= n; i++)
801      {  beg = A_ptr[i], end = A_ptr[i+1];
802         for (t = beg; t < end; t++)
803            j = A_ind[t], work[j] = A_val[t];
804         beg = U_ptr[i], end = U_ptr[i+1];
805         for (t = beg; t < end; t++)
806            j = U_ind[t], U_val[t] = work[j], work[j] = 0.0;
807         U_diag[i] = A_diag[i];
808      }
809      /* main elimination loop */
810      for (k = 1; k <= n; k++)
811      {  /* transform k-th row of U */
812         ukk = U_diag[k];
813         if (ukk > 0.0)
814            U_diag[k] = ukk = sqrt(ukk);
815         else
816            U_diag[k] = ukk = DBL_MAX, count++;
817         /* (work) := (transformed k-th row) */
818         beg = U_ptr[k], end = U_ptr[k+1];
819         for (t = beg; t < end; t++)
820            work[U_ind[t]] = (U_val[t] /= ukk);
821         /* transform other rows of U */
822         for (t = beg; t < end; t++)
823         {  i = U_ind[t];
824            xassert(i > k);
825            /* (i-th row) := (i-th row) - u[k,i] * (k-th row) */
826            uki = work[i];
827            beg1 = U_ptr[i], end1 = U_ptr[i+1];
828            for (t1 = beg1; t1 < end1; t1++)
829               U_val[t1] -= uki * work[U_ind[t1]];
830            U_diag[i] -= uki * uki;
831         }
832         /* (work) := 0 */
833         for (t = beg; t < end; t++)
834            work[U_ind[t]] = 0.0;
835      }
836      xfree(work);
837      return count;
838}
839
840/*----------------------------------------------------------------------
841-- u_solve - solve upper triangular system U*x = b.
842--
843-- *Synopsis*
844--
845-- #include "glpmat.h"
846-- void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
847--    double U_diag[], double x[]);
848--
849-- *Description*
850--
851-- The routine u_solve solves an linear system U*x = b, where U is an
852-- upper triangular matrix.
853--
854-- The parameter n is the order of matrix U.
855--
856-- The matrix U without diagonal elements is specified in the arrays
857-- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
858-- of U are specified in the array U_diag, where U_diag[0] is not used,
859-- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
860-- changed on exit.
861--
862-- The right-hand side vector b is specified on entry in the array x,
863-- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
864-- the routine stores computed components of the vector of unknowns x
865-- in the array x in the same manner. */
866
867void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
868      double U_diag[], double x[])
869{     int i, t, beg, end;
870      double temp;
871      for (i = n; i >= 1; i--)
872      {  temp = x[i];
873         beg = U_ptr[i], end = U_ptr[i+1];
874         for (t = beg; t < end; t++)
875            temp -= U_val[t] * x[U_ind[t]];
876         xassert(U_diag[i] != 0.0);
877         x[i] = temp / U_diag[i];
878      }
879      return;
880}
881
882/*----------------------------------------------------------------------
883-- ut_solve - solve lower triangular system U'*x = b.
884--
885-- *Synopsis*
886--
887-- #include "glpmat.h"
888-- void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
889--    double U_diag[], double x[]);
890--
891-- *Description*
892--
893-- The routine ut_solve solves an linear system U'*x = b, where U is a
894-- matrix transposed to an upper triangular matrix.
895--
896-- The parameter n is the order of matrix U.
897--
898-- The matrix U without diagonal elements is specified in the arrays
899-- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
900-- of U are specified in the array U_diag, where U_diag[0] is not used,
901-- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
902-- changed on exit.
903--
904-- The right-hand side vector b is specified on entry in the array x,
905-- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
906-- the routine stores computed components of the vector of unknowns x
907-- in the array x in the same manner. */
908
909void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
910      double U_diag[], double x[])
911{     int i, t, beg, end;
912      double temp;
913      for (i = 1; i <= n; i++)
914      {  xassert(U_diag[i] != 0.0);
915         temp = (x[i] /= U_diag[i]);
916         if (temp == 0.0) continue;
917         beg = U_ptr[i], end = U_ptr[i+1];
918         for (t = beg; t < end; t++)
919            x[U_ind[t]] -= U_val[t] * temp;
920      }
921      return;
922}
923
924/* eof */
Note: See TracBrowser for help on using the repository browser.