glplux.c

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

/* glplux.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 "glplux.h"
#define xfault xerror
#define dmp_create_poolx(size) dmp_create_pool()

/*----------------------------------------------------------------------
// lux_create - create LU-factorization.
//
// SYNOPSIS
//
// #include "glplux.h"
// LUX *lux_create(int n);
//
// DESCRIPTION
//
// The routine lux_create creates LU-factorization data structure for
// a matrix of the order n. Initially the factorization corresponds to
// the unity matrix (F = V = P = Q = I, so A = I).
//
// RETURNS
//
// The routine returns a pointer to the created LU-factorization data
// structure, which represents the unity matrix of the order n. */

LUX *lux_create(int n)
{     LUX *lux;
      int k;
      if (n < 1)
         xfault("lux_create: n = %d; invalid parameter\n", n);
      lux = xmalloc(sizeof(LUX));
      lux->n = n;
      lux->pool = dmp_create_poolx(sizeof(LUXELM));
      lux->F_row = xcalloc(1+n, sizeof(LUXELM *));
      lux->F_col = xcalloc(1+n, sizeof(LUXELM *));
      lux->V_piv = xcalloc(1+n, sizeof(mpq_t));
      lux->V_row = xcalloc(1+n, sizeof(LUXELM *));
      lux->V_col = xcalloc(1+n, sizeof(LUXELM *));
      lux->P_row = xcalloc(1+n, sizeof(int));
      lux->P_col = xcalloc(1+n, sizeof(int));
      lux->Q_row = xcalloc(1+n, sizeof(int));
      lux->Q_col = xcalloc(1+n, sizeof(int));
      for (k = 1; k <= n; k++)
      {  lux->F_row[k] = lux->F_col[k] = NULL;
         mpq_init(lux->V_piv[k]);
         mpq_set_si(lux->V_piv[k], 1, 1);
         lux->V_row[k] = lux->V_col[k] = NULL;
         lux->P_row[k] = lux->P_col[k] = k;
         lux->Q_row[k] = lux->Q_col[k] = k;
      }
      lux->rank = n;
      return lux;
}

/*----------------------------------------------------------------------
// initialize - initialize LU-factorization data structures.
//
// This routine initializes data structures for subsequent computing
// the LU-factorization of a given matrix A, which is specified by the
// formal routine col. On exit V = A and F = P = Q = I, where I is the
// unity matrix. */

static void initialize(LUX *lux, int (*col)(void *info, int j,
      int ind[], mpq_t val[]), void *info, LUXWKA *wka)
{     int n = lux->n;
      DMP *pool = lux->pool;
      LUXELM **F_row = lux->F_row;
      LUXELM **F_col = lux->F_col;
      mpq_t *V_piv = lux->V_piv;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *P_row = lux->P_row;
      int *P_col = lux->P_col;
      int *Q_row = lux->Q_row;
      int *Q_col = lux->Q_col;
      int *R_len = wka->R_len;
      int *R_head = wka->R_head;
      int *R_prev = wka->R_prev;
      int *R_next = wka->R_next;
      int *C_len = wka->C_len;
      int *C_head = wka->C_head;
      int *C_prev = wka->C_prev;
      int *C_next = wka->C_next;
      LUXELM *fij, *vij;
      int i, j, k, len, *ind;
      mpq_t *val;
      /* F := I */
      for (i = 1; i <= n; i++)
      {  while (F_row[i] != NULL)
         {  fij = F_row[i], F_row[i] = fij->r_next;
            mpq_clear(fij->val);
            dmp_free_atom(pool, fij, sizeof(LUXELM));
         }
      }
      for (j = 1; j <= n; j++) F_col[j] = NULL;
      /* V := 0 */
      for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1);
      for (i = 1; i <= n; i++)
      {  while (V_row[i] != NULL)
         {  vij = V_row[i], V_row[i] = vij->r_next;
            mpq_clear(vij->val);
            dmp_free_atom(pool, vij, sizeof(LUXELM));
         }
      }
      for (j = 1; j <= n; j++) V_col[j] = NULL;
      /* V := A */
      ind = xcalloc(1+n, sizeof(int));
      val = xcalloc(1+n, sizeof(mpq_t));
      for (k = 1; k <= n; k++) mpq_init(val[k]);
      for (j = 1; j <= n; j++)
      {  /* obtain j-th column of matrix A */
         len = col(info, j, ind, val);
         if (!(0 <= len && len <= n))
            xfault("lux_decomp: j = %d: len = %d; invalid column length"
               "\n", j, len);
         /* copy elements of j-th column to matrix V */
         for (k = 1; k <= len; k++)
         {  /* get row index of a[i,j] */
            i = ind[k];
            if (!(1 <= i && i <= n))
               xfault("lux_decomp: j = %d: i = %d; row index out of ran"
                  "ge\n", j, i);
            /* check for duplicate indices */
            if (V_row[i] != NULL && V_row[i]->j == j)
               xfault("lux_decomp: j = %d: i = %d; duplicate row indice"
                  "s not allowed\n", j, i);
            /* check for zero value */
            if (mpq_sgn(val[k]) == 0)
               xfault("lux_decomp: j = %d: i = %d; zero elements not al"
                  "lowed\n", j, i);
            /* add new element v[i,j] = a[i,j] to V */
            vij = dmp_get_atom(pool, sizeof(LUXELM));
            vij->i = i, vij->j = j;
            mpq_init(vij->val);
            mpq_set(vij->val, val[k]);
            vij->r_prev = NULL;
            vij->r_next = V_row[i];
            vij->c_prev = NULL;
            vij->c_next = V_col[j];
            if (vij->r_next != NULL) vij->r_next->r_prev = vij;
            if (vij->c_next != NULL) vij->c_next->c_prev = vij;
            V_row[i] = V_col[j] = vij;
         }
      }
      xfree(ind);
      for (k = 1; k <= n; k++) mpq_clear(val[k]);
      xfree(val);
      /* P := Q := I */
      for (k = 1; k <= n; k++)
         P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k;
      /* the rank of A and V is not determined yet */
      lux->rank = -1;
      /* initially the entire matrix V is active */
      /* determine its row lengths */
      for (i = 1; i <= n; i++)
      {  len = 0;
         for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++;
         R_len[i] = len;
      }
      /* build linked lists of active rows */
      for (len = 0; len <= n; len++) R_head[len] = 0;
      for (i = 1; i <= n; i++)
      {  len = R_len[i];
         R_prev[i] = 0;
         R_next[i] = R_head[len];
         if (R_next[i] != 0) R_prev[R_next[i]] = i;
         R_head[len] = i;
      }
      /* determine its column lengths */
      for (j = 1; j <= n; j++)
      {  len = 0;
         for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++;
         C_len[j] = len;
      }
      /* build linked lists of active columns */
      for (len = 0; len <= n; len++) C_head[len] = 0;
      for (j = 1; j <= n; j++)
      {  len = C_len[j];
         C_prev[j] = 0;
         C_next[j] = C_head[len];
         if (C_next[j] != 0) C_prev[C_next[j]] = j;
         C_head[len] = j;
      }
      return;
}

/*----------------------------------------------------------------------
// find_pivot - choose a pivot element.
//
// This routine chooses a pivot element v[p,q] in the active submatrix
// of matrix U = P*V*Q.
//
// It is assumed that on entry the matrix U has the following partially
// triangularized form:
//
//       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 rows and columns k, k+1, ..., n belong to the active submatrix
// (elements of the active submatrix are marked by '*').
//
// Since the matrix U = P*V*Q is not stored, the routine works with the
// matrix V. It is assumed that the row-wise representation corresponds
// to the matrix V, but the column-wise representation corresponds to
// the active submatrix of the matrix V, i.e. elements of the matrix V,
// which does not belong to the active submatrix, are missing from the
// column linked lists. It is also assumed that each active row of the
// matrix V is in the set R[len], where len is number of non-zeros in
// the row, and each active column of the matrix V is in the set C[len],
// where len is number of non-zeros in the column (in the latter case
// only elements of the active submatrix are counted; such elements are
// marked by '*' on the figure above).
//
// Due to exact arithmetic any non-zero element of the active submatrix
// can be chosen as a pivot. However, to keep sparsity of the matrix V
// the routine uses Markowitz strategy, trying to choose such element
// v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1),
// where nr[p] and nc[q] are the number of non-zero elements, resp., in
// p-th row and in q-th column of the active submatrix.
//
// In order to reduce the search, i.e. not to walk through all elements
// of the active submatrix, the routine exploits a technique proposed by
// I.Duff. This technique is based on using the sets R[len] and C[len]
// of active rows and columns.
//
// On exit the routine returns a pointer to a pivot v[p,q] chosen, or
// NULL, if the active submatrix is empty. */

static LUXELM *find_pivot(LUX *lux, LUXWKA *wka)
{     int n = lux->n;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *R_len = wka->R_len;
      int *R_head = wka->R_head;
      int *R_next = wka->R_next;
      int *C_len = wka->C_len;
      int *C_head = wka->C_head;
      int *C_next = wka->C_next;
      LUXELM *piv, *some, *vij;
      int i, j, len, min_len, ncand, piv_lim = 5;
      double best, cost;
      /* nothing is chosen so far */
      piv = NULL, best = DBL_MAX, ncand = 0;
      /* if in the active submatrix there is a column that has the only
         non-zero (column singleton), choose it as a pivot */
      j = C_head[1];
      if (j != 0)
      {  xassert(C_len[j] == 1);
         piv = V_col[j];
         xassert(piv != NULL && piv->c_next == NULL);
         goto done;
      }
      /* if in the active submatrix there is a row that has the only
         non-zero (row singleton), choose it as a pivot */
      i = R_head[1];
      if (i != 0)
      {  xassert(R_len[i] == 1);
         piv = V_row[i];
         xassert(piv != NULL && piv->r_next == NULL);
         goto done;
      }
      /* there are no singletons in the active submatrix; walk through
         other non-empty rows and columns */
      for (len = 2; len <= n; len++)
      {  /* consider active columns having len non-zeros */
         for (j = C_head[len]; j != 0; j = C_next[j])
         {  /* j-th column has len non-zeros */
            /* find an element in the row of minimal length */
            some = NULL, min_len = INT_MAX;
            for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
            {  if (min_len > R_len[vij->i])
                  some = vij, min_len = R_len[vij->i];
               /* if Markowitz cost of this element is not greater than
                  (len-1)**2, it can be chosen right now; this heuristic
                  reduces the search and works well in many cases */
               if (min_len <= len)
               {  piv = some;
                  goto done;
               }
            }
            /* j-th column has been scanned */
            /* the minimal element found is a next pivot candidate */
            xassert(some != NULL);
            ncand++;
            /* compute its Markowitz cost */
            cost = (double)(min_len - 1) * (double)(len - 1);
            /* choose between the current candidate and this element */
            if (cost < best) piv = some, best = cost;
            /* if piv_lim candidates have been considered, there is a
               doubt that a much better candidate exists; therefore it
               is the time to terminate the search */
            if (ncand == piv_lim) goto done;
         }
         /* now consider active rows having len non-zeros */
         for (i = R_head[len]; i != 0; i = R_next[i])
         {  /* i-th row has len non-zeros */
            /* find an element in the column of minimal length */
            some = NULL, min_len = INT_MAX;
            for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
            {  if (min_len > C_len[vij->j])
                  some = vij, min_len = C_len[vij->j];
               /* if Markowitz cost of this element is not greater than
                  (len-1)**2, it can be chosen right now; this heuristic
                  reduces the search and works well in many cases */
               if (min_len <= len)
               {  piv = some;
                  goto done;
               }
            }
            /* i-th row has been scanned */
            /* the minimal element found is a next pivot candidate */
            xassert(some != NULL);
            ncand++;