glpmat.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

-- on exit.
--
-- *Algorithm*
--
-- The routine min_degree is based on some subroutines from the package
-- SPARSPAK (see comments in the module glpqmd). */

void min_degree(int n, int A_ptr[], int A_ind[], int P_per[])
{     int i, j, ne, t, pos, len;
      int *xadj, *adjncy, *deg, *marker, *rchset, *nbrhd, *qsize,
         *qlink, nofsub;
      /* determine number of non-zeros in complete pattern */
      ne = A_ptr[n+1] - 1;
      ne += ne;
      /* allocate working arrays */
      xadj = xcalloc(1+n+1, sizeof(int));
      adjncy = xcalloc(1+ne, sizeof(int));
      deg = xcalloc(1+n, sizeof(int));
      marker = xcalloc(1+n, sizeof(int));
      rchset = xcalloc(1+n, sizeof(int));
      nbrhd = xcalloc(1+n, sizeof(int));
      qsize = xcalloc(1+n, sizeof(int));
      qlink = xcalloc(1+n, sizeof(int));
      /* determine row lengths in complete pattern */
      for (i = 1; i <= n; i++) xadj[i] = 0;
      for (i = 1; i <= n; i++)
      {  for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
         {  j = A_ind[t];
            xassert(i < j && j <= n);
            xadj[i]++, xadj[j]++;
         }
      }
      /* set up row pointers for complete pattern */
      pos = 1;
      for (i = 1; i <= n; i++)
         len = xadj[i], pos += len, xadj[i] = pos;
      xadj[n+1] = pos;
      xassert(pos - 1 == ne);
      /* construct complete pattern */
      for (i = 1; i <= n; i++)
      {  for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
         {  j = A_ind[t];
            adjncy[--xadj[i]] = j, adjncy[--xadj[j]] = i;
         }
      }
      /* call the main minimimum degree ordering routine */
      genqmd(&n, xadj, adjncy, P_per, P_per + n, deg, marker, rchset,
         nbrhd, qsize, qlink, &nofsub);
      /* make sure that permutation matrix P is correct */
      for (i = 1; i <= n; i++)
      {  j = P_per[i];
         xassert(1 <= j && j <= n);
         xassert(P_per[n+j] == i);
      }
      /* free working arrays */
      xfree(xadj);
      xfree(adjncy);
      xfree(deg);
      xfree(marker);
      xfree(rchset);
      xfree(nbrhd);
      xfree(qsize);
      xfree(qlink);
      return;
}

/*----------------------------------------------------------------------
-- chol_symbolic - compute Cholesky factorization (symbolic phase).
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[]);
--
-- *Description*
--
-- The routine chol_symbolic implements the symbolic phase of Cholesky
-- factorization A = U'*U, where A is a given sparse symmetric positive
-- definite matrix, U is a resultant upper triangular factor, U' is a
-- matrix transposed to U.
--
-- The parameter n is the order of matrices A and U.
--
-- The pattern of the given matrix A is specified on entry in the arrays
-- A_ptr and A_ind in storage-by-rows format. Only the upper triangular
-- part without diagonal elements (which all are assumed to be non-zero)
-- should be specified as if A were upper triangular. The arrays A_ptr
-- and A_ind are not changed on exit.
--
-- The pattern of the matrix U without diagonal elements (which all are
-- assumed to be non-zero) is stored on exit from the routine in the
-- arrays U_ptr and U_ind in storage-by-rows format. The array U_ptr
-- should be allocated on entry, however, its content is ignored. The
-- array U_ind is allocated by the routine which returns a pointer to it
-- on exit.
--
-- *Returns*
--
-- The routine returns a pointer to the array U_ind.
--
-- *Method*
--
-- The routine chol_symbolic computes the pattern of the matrix U in a
-- row-wise manner. No pivoting is used.
--
-- It is known that to compute the pattern of row k of the matrix U we
-- need to merge the pattern of row k of the matrix A and the patterns
-- of each row i of U, where u[i,k] is non-zero (these rows are already
-- computed and placed above row k).
--
-- However, to reduce the number of rows to be merged the routine uses
-- an advanced algorithm proposed in:
--
-- D.J.Rose, R.E.Tarjan, and G.S.Lueker. Algorithmic aspects of vertex
-- elimination on graphs. SIAM J. Comput. 5, 1976, 266-83.
--
-- The authors of the cited paper show that we have the same result if
-- we merge row k of the matrix A and such rows of the matrix U (among
-- rows 1, ..., k-1) whose leftmost non-diagonal non-zero element is
-- placed in k-th column. This feature signficantly reduces the number
-- of rows to be merged, especially on the final steps, where rows of
-- the matrix U become quite dense.
--
-- To determine rows, which should be merged on k-th step, for a fixed
-- time the routine uses linked lists of row numbers of the matrix U.
-- Location head[k] contains the number of a first row, whose leftmost
-- non-diagonal non-zero element is placed in column k, and location
-- next[i] contains the number of a next row with the same property as
-- row i. */

int *chol_symbolic(int n, int A_ptr[], int A_ind[], int U_ptr[])
{     int i, j, k, t, len, size, beg, end, min_j, *U_ind, *head, *next,
         *ind, *map, *temp;
      /* initially we assume that on computing the pattern of U fill-in
         will double the number of non-zeros in A */
      size = A_ptr[n+1] - 1;
      if (size < n) size = n;
      size += size;
      U_ind = xcalloc(1+size, sizeof(int));
      /* allocate and initialize working arrays */
      head = xcalloc(1+n, sizeof(int));
      for (i = 1; i <= n; i++) head[i] = 0;
      next = xcalloc(1+n, sizeof(int));
      ind = xcalloc(1+n, sizeof(int));
      map = xcalloc(1+n, sizeof(int));
      for (j = 1; j <= n; j++) map[j] = 0;
      /* compute the pattern of matrix U */
      U_ptr[1] = 1;
      for (k = 1; k <= n; k++)
      {  /* compute the pattern of k-th row of U, which is the union of
            k-th row of A and those rows of U (among 1, ..., k-1) whose
            leftmost non-diagonal non-zero is placed in k-th column */
         /* (ind) := (k-th row of A) */
         len = A_ptr[k+1] - A_ptr[k];
         memcpy(&ind[1], &A_ind[A_ptr[k]], len * sizeof(int));
         for (t = 1; t <= len; t++)
         {  j = ind[t];
            xassert(k < j && j <= n);
            map[j] = 1;
         }
         /* walk through rows of U whose leftmost non-diagonal non-zero
            is placed in k-th column */
         for (i = head[k]; i != 0; i = next[i])
         {  /* (ind) := (ind) union (i-th row of U) */
            beg = U_ptr[i], end = U_ptr[i+1];
            for (t = beg; t < end; t++)
            {  j = U_ind[t];
               if (j > k && !map[j]) ind[++len] = j, map[j] = 1;
            }
         }
         /* now (ind) is the pattern of k-th row of U */
         U_ptr[k+1] = U_ptr[k] + len;
         /* at least (U_ptr[k+1] - 1) locations should be available in
            the array U_ind */
         if (U_ptr[k+1] - 1 > size)
         {  temp = U_ind;
            size += size;
            U_ind = xcalloc(1+size, sizeof(int));
            memcpy(&U_ind[1], &temp[1], (U_ptr[k] - 1) * sizeof(int));
            xfree(temp);
         }
         xassert(U_ptr[k+1] - 1 <= size);
         /* (k-th row of U) := (ind) */
         memcpy(&U_ind[U_ptr[k]], &ind[1], len * sizeof(int));
         /* determine column index of leftmost non-diagonal non-zero in
            k-th row of U and clear the row pattern map */
         min_j = n + 1;
         for (t = 1; t <= len; t++)
         {  j = ind[t], map[j] = 0;
            if (min_j > j) min_j = j;
         }
         /* include k-th row into corresponding linked list */
         if (min_j <= n) next[k] = head[min_j], head[min_j] = k;
      }
      /* free working arrays */
      xfree(head);
      xfree(next);
      xfree(ind);
      xfree(map);
      /* reallocate the array U_ind to free unused locations */
      temp = U_ind;
      size = U_ptr[n+1] - 1;
      U_ind = xcalloc(1+size, sizeof(int));
      memcpy(&U_ind[1], &temp[1], size * sizeof(int));
      xfree(temp);
      return U_ind;
}

/*----------------------------------------------------------------------
-- chol_numeric - compute Cholesky factorization (numeric phase).
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- int chol_numeric(int n,
--    int A_ptr[], int A_ind[], double A_val[], double A_diag[],
--    int U_ptr[], int U_ind[], double U_val[], double U_diag[]);
--
-- *Description*
--
-- The routine chol_symbolic implements the numeric phase of Cholesky
-- factorization A = U'*U, where A is a given sparse symmetric positive
-- definite matrix, U is a resultant upper triangular factor, U' is a
-- matrix transposed to U.
--
-- The parameter n is the order of matrices A and U.
--
-- Upper triangular part of the matrix A without diagonal elements is
-- specified in the arrays A_ptr, A_ind, and A_val in storage-by-rows
-- format. Diagonal elements of A are specified in the array A_diag,
-- where A_diag[0] is not used, A_diag[i] = a[i,i] for i = 1, ..., n.
-- The arrays A_ptr, A_ind, A_val, and A_diag are not changed on exit.
--
-- The pattern of the matrix U without diagonal elements (previously
-- computed with the routine chol_symbolic) is specified in the arrays
-- U_ptr and U_ind, which are not changed on exit. Numeric values of
-- non-diagonal elements of U are stored in corresponding locations of
-- the array U_val, and values of diagonal elements of U are stored in
-- locations U_diag[1], ..., U_diag[n].
--
-- *Returns*
--
-- The routine returns the number of non-positive diagonal elements of
-- the matrix U which have been replaced by a huge positive number (see
-- the method description below). Zero return code means the matrix A
-- has been successfully factorized.
--
-- *Method*
--
-- The routine chol_numeric computes the matrix U in a row-wise manner
-- using standard gaussian elimination technique. No pivoting is used.
--
-- Initially the routine sets U = A, and before k-th elimination step
-- the matrix U is the following:
--
--       1       k         n
--    1  x x x x x x x x x x
--       . x x x x x x x x x
--       . . x x x x x x x x
--       . . . x x x x x x x
--    k  . . . . * * * * * *
--       . . . . * * * * * *
--       . . . . * * * * * *
--       . . . . * * * * * *
--       . . . . * * * * * *
--    n  . . . . * * * * * *
--
-- where 'x' are elements of already computed rows, '*' are elements of
-- the active submatrix. (Note that the lower triangular part of the
-- active submatrix being symmetric is not stored and diagonal elements
-- are stored separately in the array U_diag.)
--
-- The matrix A is assumed to be positive definite. However, if it is
-- close to semi-definite, on some elimination step a pivot u[k,k] may
-- happen to be non-positive due to round-off errors. In this case the
-- routine uses a technique proposed in:
--
-- S.J.Wright. The Cholesky factorization in interior-point and barrier
-- methods. Preprint MCS-P600-0596, Mathematics and Computer Science
-- Division, Argonne National Laboratory, Argonne, Ill., May 1996.
--
-- The routine just replaces non-positive u[k,k] by a huge positive
-- number. This involves non-diagonal elements in k-th row of U to be
-- close to zero that, in turn, involves k-th component of a solution
-- vector to be close to zero. Note, however, that this technique works
-- only if the system A*x = b is consistent. */

int chol_numeric(int n,
      int A_ptr[], int A_ind[], double A_val[], double A_diag[],
      int U_ptr[], int U_ind[], double U_val[], double U_diag[])
{     int i, j, k, t, t1, beg, end, beg1, end1, count = 0;
      double ukk, uki, *work;
      work = xcalloc(1+n, sizeof(double));
      for (j = 1; j <= n; j++) work[j] = 0.0;
      /* U := (upper triangle of A) */
      /* note that the upper traingle of A is a subset of U */
      for (i = 1; i <= n; i++)
      {  beg = A_ptr[i], end = A_ptr[i+1];
         for (t = beg; t < end; t++)
            j = A_ind[t], work[j] = A_val[t];
         beg = U_ptr[i], end = U_ptr[i+1];
         for (t = beg; t < end; t++)
            j = U_ind[t], U_val[t] = work[j], work[j] = 0.0;
         U_diag[i] = A_diag[i];
      }
      /* main elimination loop */
      for (k = 1; k <= n; k++)
      {  /* transform k-th row of U */
         ukk = U_diag[k];
         if (ukk > 0.0)
            U_diag[k] = ukk = sqrt(ukk);
         else
            U_diag[k] = ukk = DBL_MAX, count++;
         /* (work) := (transformed k-th row) */
         beg = U_ptr[k], end = U_ptr[k+1];
         for (t = beg; t < end; t++)
            work[U_ind[t]] = (U_val[t] /= ukk);
         /* transform other rows of U */
         for (t = beg; t < end; t++)
         {  i = U_ind[t];
            xassert(i > k);
            /* (i-th row) := (i-th row) - u[k,i] * (k-th row) */
            uki = work[i];
            beg1 = U_ptr[i], end1 = U_ptr[i+1];
            for (t1 = beg1; t1 < end1; t1++)
               U_val[t1] -= uki * work[U_ind[t1]];
            U_diag[i] -= uki * uki;
         }
         /* (work) := 0 */
         for (t = beg; t < end; t++)
            work[U_ind[t]] = 0.0;
      }
      xfree(work);
      return count;
}

/*----------------------------------------------------------------------
-- u_solve - solve upper triangular system U*x = b.
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
--    double U_diag[], double x[]);
--
-- *Description*
--
-- The routine u_solve solves an linear system U*x = b, where U is an
-- upper triangular matrix.
--
-- The parameter n is the order of matrix U.
--
-- The matrix U without diagonal elements is specified in the arrays
-- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
-- of U are specified in the array U_diag, where U_diag[0] is not used,
-- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
-- changed on exit.
--
-- The right-hand side vector b is specified on entry in the array x,
-- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
-- the routine stores computed components of the vector of unknowns x
-- in the array x in the same manner. */

void u_solve(int n, int U_ptr[], int U_ind[], double U_val[],
      double U_diag[], double x[])
{     int i, t, beg, end;
      double temp;
      for (i = n; i >= 1; i--)
      {  temp = x[i];
         beg = U_ptr[i], end = U_ptr[i+1];
         for (t = beg; t < end; t++)
            temp -= U_val[t] * x[U_ind[t]];
         xassert(U_diag[i] != 0.0);
         x[i] = temp / U_diag[i];
      }
      return;
}

/*----------------------------------------------------------------------
-- ut_solve - solve lower triangular system U'*x = b.
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
--    double U_diag[], double x[]);
--
-- *Description*
--
-- The routine ut_solve solves an linear system U'*x = b, where U is a
-- matrix transposed to an upper triangular matrix.
--
-- The parameter n is the order of matrix U.
--
-- The matrix U without diagonal elements is specified in the arrays
-- U_ptr, U_ind, and U_val in storage-by-rows format. Diagonal elements
-- of U are specified in the array U_diag, where U_diag[0] is not used,
-- U_diag[i] = u[i,i] for i = 1, ..., n. All these four arrays are not
-- changed on exit.
--
-- The right-hand side vector b is specified on entry in the array x,
-- where x[0] is not used, and x[i] = b[i] for i = 1, ..., n. On exit
-- the routine stores computed components of the vector of unknowns x
-- in the array x in the same manner. */

void ut_solve(int n, int U_ptr[], int U_ind[], double U_val[],
      double U_diag[], double x[])
{     int i, t, beg, end;
      double temp;
      for (i = 1; i <= n; i++)
      {  xassert(U_diag[i] != 0.0);
         temp = (x[i] /= U_diag[i]);
         if (temp == 0.0) continue;
         beg = U_ptr[i], end = U_ptr[i+1];
         for (t = beg; t < end; t++)
            x[U_ind[t]] -= U_val[t] * temp;
      }
      return;
}

/* eof */