glpmat.c

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

/* glpmat.c */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glplib.h"
#include "glpmat.h"
#include "glpqmd.h"

/*----------------------------------------------------------------------
-- check_fvs - check sparse vector in full-vector storage format.
--
-- SYNOPSIS
--
-- #include "glpmat.h"
-- int check_fvs(int n, int nnz, int ind[], double vec[]);
--
-- DESCRIPTION
--
-- The routine check_fvs checks if a given vector of dimension n in
-- full-vector storage format has correct representation.
--
-- RETURNS
--
-- The routine returns one of the following codes:
--
-- 0 - the vector is correct;
-- 1 - the number of elements (n) is negative;
-- 2 - the number of non-zero elements (nnz) is negative;
-- 3 - some element index is out of range;
-- 4 - some element index is duplicate;
-- 5 - some non-zero element is out of pattern. */

int check_fvs(int n, int nnz, int ind[], double vec[])
{     int i, t, ret, *flag = NULL;
      /* check the number of elements */
      if (n < 0)
      {  ret = 1;
         goto done;
      }
      /* check the number of non-zero elements */
      if (nnz < 0)
      {  ret = 2;
         goto done;
      }
      /* check vector indices */
      flag = xcalloc(1+n, sizeof(int));
      for (i = 1; i <= n; i++) flag[i] = 0;
      for (t = 1; t <= nnz; t++)
      {  i = ind[t];
         if (!(1 <= i && i <= n))
         {  ret = 3;
            goto done;
         }
         if (flag[i])
         {  ret = 4;
            goto done;
         }
         flag[i] = 1;
      }
      /* check vector elements */
      for (i = 1; i <= n; i++)
      {  if (!flag[i] && vec[i] != 0.0)
         {  ret = 5;
            goto done;
         }
      }
      /* the vector is ok */
      ret = 0;
done: if (flag != NULL) xfree(flag);
      return ret;
}

/*----------------------------------------------------------------------
-- check_pattern - check pattern of sparse matrix.
--
-- SYNOPSIS
--
-- #include "glpmat.h"
-- int check_pattern(int m, int n, int A_ptr[], int A_ind[]);
--
-- DESCRIPTION
--
-- The routine check_pattern checks the pattern of a given mxn matrix
-- in storage-by-rows format.
--
-- RETURNS
--
-- The routine returns one of the following codes:
--
-- 0 - the pattern is correct;
-- 1 - the number of rows (m) is negative;
-- 2 - the number of columns (n) is negative;
-- 3 - A_ptr[1] is not 1;
-- 4 - some column index is out of range;
-- 5 - some column indices are duplicate. */

int check_pattern(int m, int n, int A_ptr[], int A_ind[])
{     int i, j, ptr, ret, *flag = NULL;
      /* check the number of rows */
      if (m < 0)
      {  ret = 1;
         goto done;
      }
      /* check the number of columns */
      if (n < 0)
      {  ret = 2;
         goto done;
      }
      /* check location A_ptr[1] */
      if (A_ptr[1] != 1)
      {  ret = 3;
         goto done;
      }
      /* check row patterns */
      flag = xcalloc(1+n, sizeof(int));
      for (j = 1; j <= n; j++) flag[j] = 0;
      for (i = 1; i <= m; i++)
      {  /* check pattern of row i */
         for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
         {  j = A_ind[ptr];
            /* check column index */
            if (!(1 <= j && j <= n))
            {  ret = 4;
               goto done;
            }
            /* check for duplication */
            if (flag[j])
            {  ret = 5;
               goto done;
            }
            flag[j] = 1;
         }
         /* clear flags */
         for (ptr = A_ptr[i]; ptr < A_ptr[i+1]; ptr++)
         {  j = A_ind[ptr];
            flag[j] = 0;
         }
      }
      /* the pattern is ok */
      ret = 0;
done: if (flag != NULL) xfree(flag);
      return ret;
}

/*----------------------------------------------------------------------
-- transpose - transpose sparse matrix.
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- void transpose(int m, int n, int A_ptr[], int A_ind[],
--    double A_val[], int AT_ptr[], int AT_ind[], double AT_val[]);
--
-- *Description*
--
-- For a given mxn sparse matrix A the routine transpose builds a nxm
-- sparse matrix A' which is a matrix transposed to A.
--
-- The arrays A_ptr, A_ind, and A_val specify a given mxn matrix A to
-- be transposed in storage-by-rows format. The parameter A_val can be
-- NULL, in which case numeric values are not copied. The arrays A_ptr,
-- A_ind, and A_val are not changed on exit.
--
-- On entry the arrays AT_ptr, AT_ind, and AT_val must be allocated,
-- but their content is ignored. On exit the routine stores a resultant
-- nxm matrix A' in these arrays in storage-by-rows format. Note that
-- if the parameter A_val is NULL, the array AT_val is not used.
--
-- The routine transpose has a side effect that elements in rows of the
-- resultant matrix A' follow in ascending their column indices. */

void transpose(int m, int n, int A_ptr[], int A_ind[], double A_val[],
      int AT_ptr[], int AT_ind[], double AT_val[])
{     int i, j, t, beg, end, pos, len;
      /* determine row lengths of resultant matrix */
      for (j = 1; j <= n; j++) AT_ptr[j] = 0;
      for (i = 1; i <= m; i++)
      {  beg = A_ptr[i], end = A_ptr[i+1];
         for (t = beg; t < end; t++) AT_ptr[A_ind[t]]++;
      }
      /* set up row pointers of resultant matrix */
      pos = 1;
      for (j = 1; j <= n; j++)
         len = AT_ptr[j], pos += len, AT_ptr[j] = pos;
      AT_ptr[n+1] = pos;
      /* build resultant matrix */
      for (i = m; i >= 1; i--)
      {  beg = A_ptr[i], end = A_ptr[i+1];
         for (t = beg; t < end; t++)
         {  pos = --AT_ptr[A_ind[t]];
            AT_ind[pos] = i;
            if (A_val != NULL) AT_val[pos] = A_val[t];
         }
      }
      return;
}

/*----------------------------------------------------------------------
-- adat_symbolic - compute S = P*A*D*A'*P' (symbolic phase).
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- int *adat_symbolic(int m, int n, int P_per[], int A_ptr[],
--    int A_ind[], int S_ptr[]);
--
-- *Description*
--
-- The routine adat_symbolic implements the symbolic phase to compute
-- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
-- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
-- transposed to A, P' is an inverse of P.
--
-- The parameter m is the number of rows in A and the order of P.
--
-- The parameter n is the number of columns in A and the order of D.
--
-- The array P_per specifies permutation matrix P. It is not changed on
-- exit.
--
-- The arrays A_ptr and A_ind specify the pattern of matrix A. They are
-- not changed on exit.
--
-- On exit the routine stores the pattern of upper triangular part of
-- matrix S without diagonal elements in the arrays S_ptr and S_ind in
-- storage-by-rows format. The array S_ptr should be allocated on entry,
-- however, its content is ignored. The array S_ind is allocated by the
-- routine itself which returns a pointer to it.
--
-- *Returns*
--
-- The routine returns a pointer to the array S_ind. */

int *adat_symbolic(int m, int n, int P_per[], int A_ptr[], int A_ind[],
      int S_ptr[])
{     int i, j, t, ii, jj, tt, k, size, len;
      int *S_ind, *AT_ptr, *AT_ind, *ind, *map, *temp;
      /* build the pattern of A', which is a matrix transposed to A, to
         efficiently access A in column-wise manner */
      AT_ptr = xcalloc(1+n+1, sizeof(int));
      AT_ind = xcalloc(A_ptr[m+1], sizeof(int));
      transpose(m, n, A_ptr, A_ind, NULL, AT_ptr, AT_ind, NULL);
      /* allocate the array S_ind */
      size = A_ptr[m+1] - 1;
      if (size < m) size = m;
      S_ind = xcalloc(1+size, sizeof(int));
      /* allocate and initialize working arrays */
      ind = xcalloc(1+m, sizeof(int));
      map = xcalloc(1+m, sizeof(int));
      for (jj = 1; jj <= m; jj++) map[jj] = 0;
      /* compute pattern of S; note that symbolically S = B*B', where
         B = P*A, B' is matrix transposed to B */
      S_ptr[1] = 1;
      for (ii = 1; ii <= m; ii++)
      {  /* compute pattern of ii-th row of S */
         len = 0;
         i = P_per[ii]; /* i-th row of A = ii-th row of B */
         for (t = A_ptr[i]; t < A_ptr[i+1]; t++)
         {  k = A_ind[t];
            /* walk through k-th column of A */
            for (tt = AT_ptr[k]; tt < AT_ptr[k+1]; tt++)
            {  j = AT_ind[tt];
               jj = P_per[m+j]; /* j-th row of A = jj-th row of B */
               /* a[i,k] != 0 and a[j,k] != 0 ergo s[ii,jj] != 0 */
               if (ii < jj && !map[jj]) ind[++len] = jj, map[jj] = 1;
            }
         }
         /* now (ind) is pattern of ii-th row of S */
         S_ptr[ii+1] = S_ptr[ii] + len;
         /* at least (S_ptr[ii+1] - 1) locations should be available in
            the array S_ind */
         if (S_ptr[ii+1] - 1 > size)
         {  temp = S_ind;
            size += size;
            S_ind = xcalloc(1+size, sizeof(int));
            memcpy(&S_ind[1], &temp[1], (S_ptr[ii] - 1) * sizeof(int));
            xfree(temp);
         }
         xassert(S_ptr[ii+1] - 1 <= size);
         /* (ii-th row of S) := (ind) */
         memcpy(&S_ind[S_ptr[ii]], &ind[1], len * sizeof(int));
         /* clear the row pattern map */
         for (t = 1; t <= len; t++) map[ind[t]] = 0;
      }
      /* free working arrays */
      xfree(AT_ptr);
      xfree(AT_ind);
      xfree(ind);
      xfree(map);
      /* reallocate the array S_ind to free unused locations */
      temp = S_ind;
      size = S_ptr[m+1] - 1;
      S_ind = xcalloc(1+size, sizeof(int));
      memcpy(&S_ind[1], &temp[1], size * sizeof(int));
      xfree(temp);
      return S_ind;
}

/*----------------------------------------------------------------------
-- adat_numeric - compute S = P*A*D*A'*P' (numeric phase).
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- void adat_numeric(int m, int n, int P_per[],
--    int A_ptr[], int A_ind[], double A_val[], double D_diag[],
--    int S_ptr[], int S_ind[], double S_val[], double S_diag[]);
--
-- *Description*
--
-- The routine adat_numeric implements the numeric phase to compute
-- symmetric matrix S = P*A*D*A'*P', where P is a permutation matrix,
-- A is a given sparse matrix, D is a diagonal matrix, A' is a matrix
-- transposed to A, P' is an inverse of P.
--
-- The parameter m is the number of rows in A and the order of P.
--
-- The parameter n is the number of columns in A and the order of D.
--
-- The matrix P is specified in the array P_per, which is not changed
-- on exit.
--
-- The matrix A is specified in the arrays A_ptr, A_ind, and A_val in
-- storage-by-rows format. These arrays are not changed on exit.
--
-- Diagonal elements of the matrix D are specified in the array D_diag,
-- where D_diag[0] is not used, D_diag[i] = d[i,i] for i = 1, ..., n.
-- The array D_diag is not changed on exit.
--
-- The pattern of the upper triangular part of the matrix S without
-- diagonal elements (previously computed by the routine adat_symbolic)
-- is specified in the arrays S_ptr and S_ind, which are not changed on
-- exit. Numeric values of non-diagonal elements of S are stored in
-- corresponding locations of the array S_val, and values of diagonal
-- elements of S are stored in locations S_diag[1], ..., S_diag[n]. */

void adat_numeric(int m, int n, int P_per[],
      int A_ptr[], int A_ind[], double A_val[], double D_diag[],
      int S_ptr[], int S_ind[], double S_val[], double S_diag[])
{     int i, j, t, ii, jj, tt, beg, end, beg1, end1, k;
      double sum, *work;
      work = xcalloc(1+n, sizeof(double));
      for (j = 1; j <= n; j++) work[j] = 0.0;
      /* compute S = B*D*B', where B = P*A, B' is a matrix transposed
         to B */
      for (ii = 1; ii <= m; ii++)
      {  i = P_per[ii]; /* i-th row of A = ii-th row of B */
         /* (work) := (i-th row of A) */
         beg = A_ptr[i], end = A_ptr[i+1];
         for (t = beg; t < end; t++)
            work[A_ind[t]] = A_val[t];
         /* compute ii-th row of S */
         beg = S_ptr[ii], end = S_ptr[ii+1];
         for (t = beg; t < end; t++)
         {  jj = S_ind[t];
            j = P_per[jj]; /* j-th row of A = jj-th row of B */
            /* s[ii,jj] := sum a[i,k] * d[k,k] * a[j,k] */
            sum = 0.0;
            beg1 = A_ptr[j], end1 = A_ptr[j+1];
            for (tt = beg1; tt < end1; tt++)
            {  k = A_ind[tt];
               sum += work[k] * D_diag[k] * A_val[tt];
            }
            S_val[t] = sum;
         }
         /* s[ii,ii] := sum a[i,k] * d[k,k] * a[i,k] */
         sum = 0.0;
         beg = A_ptr[i], end = A_ptr[i+1];
         for (t = beg; t < end; t++)
         {  k = A_ind[t];
            sum += A_val[t] * D_diag[k] * A_val[t];
            work[k] = 0.0;
         }
         S_diag[ii] = sum;
      }
      xfree(work);
      return;
}

/*----------------------------------------------------------------------
-- min_degree - minimum degree ordering.
--
-- *Synopsis*
--
-- #include "glpmat.h"
-- void min_degree(int n, int A_ptr[], int A_ind[], int P_per[]);
--
-- *Description*
--
-- The routine min_degree uses the minimum degree ordering algorithm
-- to find a permutation matrix P for a given sparse symmetric positive
-- matrix A which minimizes the number of non-zeros in upper triangular
-- factor U for Cholesky factorization P*A*P' = U'*U.
--
-- The parameter n is the order of matrices A and P.
--
-- 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 permutation matrix P is stored by the routine in the array P_per