glpini01.c

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

/* glpini01.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 "glpini.h"

/*----------------------------------------------------------------------
-- triang - find maximal triangular part of a rectangular matrix.
--
-- *Synopsis*
--
-- int triang(int m, int n,
--    void *info, int (*mat)(void *info, int k, int ndx[]),
--    int rn[], int cn[]);
--
-- *Description*
--
-- For a given rectangular (sparse) matrix A with m rows and n columns
-- the routine triang tries to find such permutation matrices P and Q
-- that the first rows and columns of the matrix B = P*A*Q form a lower
-- triangular submatrix of as greatest size as possible:
--
--                   1                       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
--    B = P*A*Q =    * * * * * . . 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 x x x x x x x x x x
--                m  x x x x x x x x x x x x x
--
-- where: '*' - elements of the lower triangular part, '.' - structural
-- zeros, 'x' - other (either non-zero or zero) elements.
--
-- The parameter info is a transit pointer passed to the formal routine
-- mat (see below).
--
-- The formal routine mat specifies the given matrix A in both row- and
-- column-wise formats. In order to obtain an i-th row of the matrix A
-- the routine triang calls the routine mat with the parameter k = +i,
-- 1 <= i <= m. In response the routine mat should store column indices
-- of (non-zero) elements of the i-th row to the locations ndx[1], ...,
-- ndx[len], where len is number of non-zeros in the i-th row returned
-- on exit. Analogously, in order to obtain a j-th column of the matrix
-- A, the routine mat is called with the parameter k = -j, 1 <= j <= n,
-- and should return pattern of the j-th column in the same way as for
-- row patterns. Note that the routine mat may be called more than once
-- for the same rows and columns.
--
-- On exit the routine computes two resultant arrays rn and cn, which
-- define the permutation matrices P and Q, respectively. The array rn
-- should have at least 1+m locations, where rn[i] = i' (1 <= i <= m)
-- means that i-th row of the original matrix A corresponds to i'-th row
-- of the matrix B = P*A*Q. Similarly, the array cn should have at least
-- 1+n locations, where cn[j] = j' (1 <= j <= n) means that j-th column
-- of the matrix A corresponds to j'-th column of the matrix B.
--
-- *Returns*
--
-- The routine triang returns the size of the lower tringular part of
-- the matrix B = P*A*Q (see the figure above).
--
-- *Complexity*
--
-- The time complexity of the routine triang is O(nnz), where nnz is
-- number of non-zeros in the given matrix A.
--
-- *Algorithm*
--
-- The routine triang starts from the matrix B = P*Q*A, where P and Q
-- are unity matrices, so initially B = A.
--
-- Before the next iteration B = (B1 | B2 | B3), where B1 is partially
-- built a lower triangular submatrix, B2 is the active submatrix, and
-- B3 is a submatrix that contains rejected columns. Thus, the current
-- matrix B looks like follows (initially k1 = 1 and k2 = n):
--
--       1         k1         k2         n
--    1  x . . . . . . . . . . . . . # # #
--       x x . . . . . . . . . . . . # # #
--       x x x . . . . . . . . . . # # # #
--       x x x x . . . . . . . . . # # # #
--       x x x x x . . . . . . . # # # # #
--    k1 x x x x x * * * * * * * # # # # #
--       x x x x x * * * * * * * # # # # #
--       x x x x x * * * * * * * # # # # #
--       x x x x x * * * * * * * # # # # #
--    m  x x x x x * * * * * * * # # # # #
--       <--B1---> <----B2-----> <---B3-->
--
-- On each iteartion the routine looks for a singleton row, i.e. some
-- row that has the only non-zero in the active submatrix B2. If such
-- row exists and the corresponding non-zero is b[i,j], where (by the
-- definition) k1 <= i <= m and k1 <= j <= k2, the routine permutes
-- k1-th and i-th rows and k1-th and j-th columns of the matrix B (in
-- order to place the element in the position b[k1,k1]), removes the
-- k1-th column from the active submatrix B2, and adds this column to
-- the submatrix B1. If no row singletons exist, but B2 is not empty
-- yet, the routine chooses a j-th column, which has maximal number of
-- non-zeros among other columns of B2, removes this column from B2 and
-- adds it to the submatrix B3 in the hope that new row singletons will
-- appear in the active submatrix. */

static int triang(int m, int n,
      void *info, int (*mat)(void *info, int k, int ndx[]),
      int rn[], int cn[])
{     int *ndx; /* int ndx[1+max(m,n)]; */
      /* this array is used for querying row and column patterns of the
         given matrix A (the third parameter to the routine mat) */
      int *rs_len; /* int rs_len[1+m]; */
      /* rs_len[0] is not used;
         rs_len[i], 1 <= i <= m, is number of non-zeros in the i-th row
         of the matrix A, which (non-zeros) belong to the current active
         submatrix */
      int *rs_head; /* int rs_head[1+n]; */
      /* rs_head[len], 0 <= len <= n, is the number i of the first row
         of the matrix A, for which rs_len[i] = len */
      int *rs_prev; /* int rs_prev[1+m]; */
      /* rs_prev[0] is not used;
         rs_prev[i], 1 <= i <= m, is a number i' of the previous row of
         the matrix A, for which rs_len[i] = rs_len[i'] (zero marks the
         end of this linked list) */
      int *rs_next; /* int rs_next[1+m]; */
      /* rs_next[0] is not used;
         rs_next[i], 1 <= i <= m, is a number i' of the next row of the
         matrix A, for which rs_len[i] = rs_len[i'] (zero marks the end
         this linked list) */
      int cs_head;
      /* is a number j of the first column of the matrix A, which has
         maximal number of non-zeros among other columns */
      int *cs_prev; /* cs_prev[1+n]; */
      /* cs_prev[0] is not used;
         cs_prev[j], 1 <= j <= n, is a number of the previous column of
         the matrix A with the same or greater number of non-zeros than
         in the j-th column (zero marks the end of this linked list) */
      int *cs_next; /* cs_next[1+n]; */
      /* cs_next[0] is not used;
         cs_next[j], 1 <= j <= n, is a number of the next column of
         the matrix A with the same or lesser number of non-zeros than
         in the j-th column (zero marks the end of this linked list) */
      int i, j, ii, jj, k1, k2, len, t, size = 0;
      int *head, *rn_inv, *cn_inv;
      if (!(m > 0 && n > 0))
         xerror("triang: m = %d; n = %d; invalid dimension\n", m, n);
      /* allocate working arrays */
      ndx = xcalloc(1+(m >= n ? m : n), sizeof(int));
      rs_len = xcalloc(1+m, sizeof(int));
      rs_head = xcalloc(1+n, sizeof(int));
      rs_prev = xcalloc(1+m, sizeof(int));
      rs_next = xcalloc(1+m, sizeof(int));
      cs_prev = xcalloc(1+n, sizeof(int));
      cs_next = xcalloc(1+n, sizeof(int));
      /* build linked lists of columns of the matrix A with the same
         number of non-zeros */
      head = rs_len; /* currently rs_len is used as working array */
      for (len = 0; len <= m; len ++) head[len] = 0;
      for (j = 1; j <= n; j++)
      {  /* obtain length of the j-th column */
         len = mat(info, -j, ndx);
         xassert(0 <= len && len <= m);
         /* include the j-th column in the corresponding linked list */
         cs_prev[j] = head[len];
         head[len] = j;
      }
      /* merge all linked lists of columns in one linked list, where
         columns are ordered by descending of their lengths */
      cs_head = 0;
      for (len = 0; len <= m; len++)
      {  for (j = head[len]; j != 0; j = cs_prev[j])
         {  cs_next[j] = cs_head;
            cs_head = j;
         }
      }
      jj = 0;
      for (j = cs_head; j != 0; j = cs_next[j])
      {  cs_prev[j] = jj;
         jj = j;
      }
      /* build initial doubly linked lists of rows of the matrix A with
         the same number of non-zeros */
      for (len = 0; len <= n; len++) rs_head[len] = 0;
      for (i = 1; i <= m; i++)
      {  /* obtain length of the i-th row */
         rs_len[i] = len = mat(info, +i, ndx);
         xassert(0 <= len && len <= n);
         /* include the i-th row in the correspondng linked list */
         rs_prev[i] = 0;
         rs_next[i] = rs_head[len];
         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
         rs_head[len] = i;
      }
      /* initially all rows and columns of the matrix A are active */
      for (i = 1; i <= m; i++) rn[i] = 0;
      for (j = 1; j <= n; j++) cn[j] = 0;
      /* set initial bounds of the active submatrix */
      k1 = 1, k2 = n;
      /* main loop starts here */
      while (k1 <= k2)
      {  i = rs_head[1];
         if (i != 0)
         {  /* the i-th row of the matrix A is a row singleton, since
               it has the only non-zero in the active submatrix */
            xassert(rs_len[i] == 1);
            /* determine the number j of an active column of the matrix
               A, in which this non-zero is placed */
            j = 0;
            t = mat(info, +i, ndx);
            xassert(0 <= t && t <= n);
            for (t = t; t >= 1; t--)
            {  jj = ndx[t];
               xassert(1 <= jj && jj <= n);
               if (cn[jj] == 0)
               {  xassert(j == 0);
                  j = jj;
               }
            }
            xassert(j != 0);
            /* the singleton is a[i,j]; move a[i,j] to the position
               b[k1,k1] of the matrix B */
            rn[i] = cn[j] = k1;
            /* shift the left bound of the active submatrix */
            k1++;
            /* increase the size of the lower triangular part */
            size++;
         }
         else
         {  /* the current active submatrix has no row singletons */
            /* remove an active column with maximal number of non-zeros
               from the active submatrix */
            j = cs_head;
            xassert(j != 0);
            cn[j] = k2;
            /* shift the right bound of the active submatrix */
            k2--;
         }
         /* the j-th column of the matrix A has been removed from the
            active submatrix */
         /* remove the j-th column from the linked list */
         if (cs_prev[j] == 0)
            cs_head = cs_next[j];
         else
            cs_next[cs_prev[j]] = cs_next[j];
         if (cs_next[j] == 0)
            /* nop */;
         else
            cs_prev[cs_next[j]] = cs_prev[j];
         /* go through non-zeros of the j-th columns and update active
            lengths of the corresponding rows */
         t = mat(info, -j, ndx);
         xassert(0 <= t && t <= m);