glpini01.c

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

         for (t = t; t >= 1; t--)
         {  i = ndx[t];
            xassert(1 <= i && i <= m);
            /* the non-zero a[i,j] has left the active submatrix */
            len = rs_len[i];
            xassert(len >= 1);
            /* remove the i-th row from the linked list of rows with
               active length len */
            if (rs_prev[i] == 0)
               rs_head[len] = rs_next[i];
            else
               rs_next[rs_prev[i]] = rs_next[i];
            if (rs_next[i] == 0)
               /* nop */;
            else
               rs_prev[rs_next[i]] = rs_prev[i];
            /* decrease the active length of the i-th row */
            rs_len[i] = --len;
            /* return the i-th row to the corresponding 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;
         }
      }
      /* other rows of the matrix A, which are still active, correspond
         to rows k1, ..., m of the matrix B (in arbitrary order) */
      for (i = 1; i <= m; i++) if (rn[i] == 0) rn[i] = k1++;
      /* but for columns this is not needed, because now the submatrix
         B2 has no columns */
      for (j = 1; j <= n; j++) xassert(cn[j] != 0);
      /* perform some optional checks */
      /* make sure that rn is a permutation of {1, ..., m} and cn is a
         permutation of {1, ..., n} */
      rn_inv = rs_len; /* used as working array */
      for (ii = 1; ii <= m; ii++) rn_inv[ii] = 0;
      for (i = 1; i <= m; i++)
      {  ii = rn[i];
         xassert(1 <= ii && ii <= m);
         xassert(rn_inv[ii] == 0);
         rn_inv[ii] = i;
      }
      cn_inv = rs_head; /* used as working array */
      for (jj = 1; jj <= n; jj++) cn_inv[jj] = 0;
      for (j = 1; j <= n; j++)
      {  jj = cn[j];
         xassert(1 <= jj && jj <= n);
         xassert(cn_inv[jj] == 0);
         cn_inv[jj] = j;
      }
      /* make sure that the matrix B = P*A*Q really has the form, which
         was declared */
      for (ii = 1; ii <= size; ii++)
      {  int diag = 0;
         i = rn_inv[ii];
         t = mat(info, +i, ndx);
         xassert(0 <= t && t <= n);
         for (t = t; t >= 1; t--)
         {  j = ndx[t];
            xassert(1 <= j && j <= n);
            jj = cn[j];
            if (jj <= size) xassert(jj <= ii);
            if (jj == ii)
            {  xassert(!diag);
               diag = 1;
            }
         }
         xassert(diag);
      }
      /* free working arrays */
      xfree(ndx);
      xfree(rs_len);
      xfree(rs_head);
      xfree(rs_prev);
      xfree(rs_next);
      xfree(cs_prev);
      xfree(cs_next);
      /* return to the calling program */
      return size;
}

/*----------------------------------------------------------------------
-- adv_basis - construct advanced initial LP basis.
--
-- *Synopsis*
--
-- #include "glpini.h"
-- void adv_basis(glp_prob *lp);
--
-- *Description*
--
-- The routine adv_basis constructs an advanced initial basis for an LP
-- problem object, which the parameter lp points to.
--
-- In order to build the initial basis the routine does the following:
--
-- 1) includes in the basis all non-fixed auxiliary variables;
--
-- 2) includes in the basis as many as possible non-fixed structural
--    variables preserving triangular form of the basis matrix;
--
-- 3) includes in the basis appropriate (fixed) auxiliary variables
--    in order to complete the basis.
--
-- As a result the initial basis has minimum of fixed variables and the
-- corresponding basis matrix is triangular. */

static int mat(void *info, int k, int ndx[])
{     /* this auxiliary routine returns the pattern of a given row or
         a given column of the augmented constraint matrix A~ = (I|-A),
         in which columns of fixed variables are implicitly cleared */
      LPX *lp = info;
      int m = lpx_get_num_rows(lp);
      int n = lpx_get_num_cols(lp);
      int typx, i, j, lll, len = 0;
      if (k > 0)
      {  /* the pattern of the i-th row is required */
         i = +k;
         xassert(1 <= i && i <= m);
#if 0 /* 22/XII-2003 */
         /* if the auxiliary variable x[i] is non-fixed, include its
            element (placed in the i-th column) in the pattern */
         lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
         if (typx != LPX_FX) ndx[++len] = i;
         /* include in the pattern elements placed in columns, which
            correspond to non-fixed structural varables */
         i_beg = aa_ptr[i];
         i_end = i_beg + aa_len[i] - 1;
         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
         {  j = m + sv_ndx[i_ptr];
            lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
            if (typx != LPX_FX) ndx[++len] = j;
         }
#else
         lll = lpx_get_mat_row(lp, i, ndx, NULL);
         for (k = 1; k <= lll; k++)
         {  lpx_get_col_bnds(lp, ndx[k], &typx, NULL, NULL);
            if (typx != LPX_FX) ndx[++len] = m + ndx[k];
         }
         lpx_get_row_bnds(lp, i, &typx, NULL, NULL);
         if (typx != LPX_FX) ndx[++len] = i;
#endif
      }
      else
      {  /* the pattern of the j-th column is required */
         j = -k;
         xassert(1 <= j && j <= m+n);
         /* if the (auxiliary or structural) variable x[j] is fixed,
            the pattern of its column is empty */
         if (j <= m)
            lpx_get_row_bnds(lp, j, &typx, NULL, NULL);
         else
            lpx_get_col_bnds(lp, j-m, &typx, NULL, NULL);
         if (typx != LPX_FX)
         {  if (j <= m)
            {  /* x[j] is non-fixed auxiliary variable */
               ndx[++len] = j;
            }
            else
            {  /* x[j] is non-fixed structural variables */
#if 0 /* 22/XII-2003 */
               j_beg = aa_ptr[j];
               j_end = j_beg + aa_len[j] - 1;
               for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
                  ndx[++len] = sv_ndx[j_ptr];
#else
               len = lpx_get_mat_col(lp, j-m, ndx, NULL);
#endif
            }
         }
      }
      /* return the length of the row/column pattern */
      return len;
}

void adv_basis(glp_prob *lp)
{     int m = lpx_get_num_rows(lp);
      int n = lpx_get_num_cols(lp);
      int i, j, jj, k, size;
      int *rn, *cn, *rn_inv, *cn_inv;
      int typx, *tagx = xcalloc(1+m+n, sizeof(int));
      double lb, ub;
      xprintf("Crashing...\n");
      if (m == 0)
         xerror("glp_adv_basis: problem has no rows\n");
      if (n == 0)
         xerror("glp_adv_basis: problem has no columns\n");
      /* use the routine triang (see above) to find maximal triangular
         part of the augmented constraint matrix A~ = (I|-A); in order
         to prevent columns of fixed variables to be included in the
         triangular part, such columns are implictly removed from the
         matrix A~ by the routine adv_mat */
      rn = xcalloc(1+m, sizeof(int));
      cn = xcalloc(1+m+n, sizeof(int));
      size = triang(m, m+n, lp, mat, rn, cn);
      if (lpx_get_int_parm(lp, LPX_K_MSGLEV) >= 3)
         xprintf("Size of triangular part = %d\n", size);
      /* the first size rows and columns of the matrix P*A~*Q (where
         P and Q are permutation matrices defined by the arrays rn and
         cn) form a lower triangular matrix; build the arrays (rn_inv
         and cn_inv), which define the matrices inv(P) and inv(Q) */
      rn_inv = xcalloc(1+m, sizeof(int));
      cn_inv = xcalloc(1+m+n, sizeof(int));
      for (i = 1; i <= m; i++) rn_inv[rn[i]] = i;
      for (j = 1; j <= m+n; j++) cn_inv[cn[j]] = j;
      /* include the columns of the matrix A~, which correspond to the
         first size columns of the matrix P*A~*Q, in the basis */
      for (k = 1; k <= m+n; k++) tagx[k] = -1;
      for (jj = 1; jj <= size; jj++)
      {  j = cn_inv[jj];
         /* the j-th column of A~ is the jj-th column of P*A~*Q */
         tagx[j] = LPX_BS;
      }
      /* if size < m, we need to add appropriate columns of auxiliary
         variables to the basis */
      for (jj = size + 1; jj <= m; jj++)
      {  /* the jj-th column of P*A~*Q should be replaced by the column
            of the auxiliary variable, for which the only unity element
            is placed in the position [jj,jj] */
         i = rn_inv[jj];
         /* the jj-th row of P*A~*Q is the i-th row of A~, but in the
            i-th row of A~ the unity element belongs to the i-th column
            of A~; therefore the disired column corresponds to the i-th
            auxiliary variable (note that this column doesn't belong to
            the triangular part found by the routine triang) */
         xassert(1 <= i && i <= m);
         xassert(cn[i] > size);
         tagx[i] = LPX_BS;
      }
      /* free working arrays */
      xfree(rn);
      xfree(cn);
      xfree(rn_inv);
      xfree(cn_inv);
      /* build tags of non-basic variables */
      for (k = 1; k <= m+n; k++)
      {  if (tagx[k] != LPX_BS)
         {  if (k <= m)
               lpx_get_row_bnds(lp, k, &typx, &lb, &ub);
            else
               lpx_get_col_bnds(lp, k-m, &typx, &lb, &ub);
            switch (typx)
            {  case LPX_FR:
                  tagx[k] = LPX_NF; break;
               case LPX_LO:
                  tagx[k] = LPX_NL; break;
               case LPX_UP:
                  tagx[k] = LPX_NU; break;
               case LPX_DB:
                  tagx[k] =
                     (fabs(lb) <= fabs(ub) ? LPX_NL : LPX_NU);
                  break;
               case LPX_FX:
                  tagx[k] = LPX_NS; break;
               default:
                  xassert(typx != typx);
            }
         }
      }
      for (k = 1; k <= m+n; k++)
      {  if (k <= m)
            lpx_set_row_stat(lp, k, tagx[k]);
         else
            lpx_set_col_stat(lp, k-m, tagx[k]);
      }
      xfree(tagx);
      return;
}

/* eof */