glpssx01.c

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

/* glpssx01.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 "glpssx.h"
#define xfault xerror

/*----------------------------------------------------------------------
// ssx_create - create simplex solver workspace.
//
// This routine creates the workspace used by simplex solver routines,
// and returns a pointer to it.
//
// Parameters m, n, and nnz specify, respectively, the number of rows,
// columns, and non-zero constraint coefficients.
//
// This routine only allocates the memory for the workspace components,
// so the workspace needs to be saturated by data. */

SSX *ssx_create(int m, int n, int nnz)
{     SSX *ssx;
      int i, j, k;
      if (m < 1)
         xfault("ssx_create: m = %d; invalid number of rows\n", m);
      if (n < 1)
         xfault("ssx_create: n = %d; invalid number of columns\n", n);
      if (nnz < 0)
         xfault("ssx_create: nnz = %d; invalid number of non-zero const"
            "raint coefficients\n", nnz);
      ssx = xmalloc(sizeof(SSX));
      ssx->m = m;
      ssx->n = n;
      ssx->type = xcalloc(1+m+n, sizeof(int));
      ssx->lb = xcalloc(1+m+n, sizeof(mpq_t));
      for (k = 1; k <= m+n; k++) mpq_init(ssx->lb[k]);
      ssx->ub = xcalloc(1+m+n, sizeof(mpq_t));
      for (k = 1; k <= m+n; k++) mpq_init(ssx->ub[k]);
      ssx->coef = xcalloc(1+m+n, sizeof(mpq_t));
      for (k = 0; k <= m+n; k++) mpq_init(ssx->coef[k]);
      ssx->A_ptr = xcalloc(1+n+1, sizeof(int));
      ssx->A_ptr[n+1] = nnz+1;
      ssx->A_ind = xcalloc(1+nnz, sizeof(int));
      ssx->A_val = xcalloc(1+nnz, sizeof(mpq_t));
      for (k = 1; k <= nnz; k++) mpq_init(ssx->A_val[k]);
      ssx->stat = xcalloc(1+m+n, sizeof(int));
      ssx->Q_row = xcalloc(1+m+n, sizeof(int));
      ssx->Q_col = xcalloc(1+m+n, sizeof(int));
      ssx->binv = bfx_create_binv();
      ssx->bbar = xcalloc(1+m, sizeof(mpq_t));
      for (i = 0; i <= m; i++) mpq_init(ssx->bbar[i]);
      ssx->pi = xcalloc(1+m, sizeof(mpq_t));
      for (i = 1; i <= m; i++) mpq_init(ssx->pi[i]);
      ssx->cbar = xcalloc(1+n, sizeof(mpq_t));
      for (j = 1; j <= n; j++) mpq_init(ssx->cbar[j]);
      ssx->rho = xcalloc(1+m, sizeof(mpq_t));
      for (i = 1; i <= m; i++) mpq_init(ssx->rho[i]);
      ssx->ap = xcalloc(1+n, sizeof(mpq_t));
      for (j = 1; j <= n; j++) mpq_init(ssx->ap[j]);
      ssx->aq = xcalloc(1+m, sizeof(mpq_t));
      for (i = 1; i <= m; i++) mpq_init(ssx->aq[i]);
      mpq_init(ssx->delta);
      return ssx;
}

/*----------------------------------------------------------------------
// ssx_factorize - factorize the current basis matrix.
//
// This routine computes factorization of the current basis matrix B
// and returns the singularity flag. If the matrix B is non-singular,
// the flag is zero, otherwise non-zero. */

static int basis_col(void *info, int j, int ind[], mpq_t val[])
{     /* this auxiliary routine provides row indices and numeric values
         of non-zero elements in j-th column of the matrix B */
      SSX *ssx = info;
      int m = ssx->m;
      int n = ssx->n;
      int *A_ptr = ssx->A_ptr;
      int *A_ind = ssx->A_ind;
      mpq_t *A_val = ssx->A_val;
      int *Q_col = ssx->Q_col;
      int k, len, ptr;
      xassert(1 <= j && j <= m);
      k = Q_col[j]; /* x[k] = xB[j] */
      xassert(1 <= k && k <= m+n);
      /* j-th column of the matrix B is k-th column of the augmented
         constraint matrix (I | -A) */
      if (k <= m)
      {  /* it is a column of the unity matrix I */
         len = 1, ind[1] = k, mpq_set_si(val[1], 1, 1);
      }
      else
      {  /* it is a column of the original constraint matrix -A */
         len = 0;
         for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
         {  len++;
            ind[len] = A_ind[ptr];
            mpq_neg(val[len], A_val[ptr]);
         }
      }
      return len;
}

int ssx_factorize(SSX *ssx)
{     int ret;
      ret = bfx_factorize(ssx->binv, ssx->m, basis_col, ssx);
      return ret;
}

/*----------------------------------------------------------------------
// ssx_get_xNj - determine value of non-basic variable.
//
// This routine determines the value of non-basic variable xN[j] in the
// current basic solution defined as follows:
//
//    0,             if xN[j] is free variable
//    lN[j],         if xN[j] is on its lower bound
//    uN[j],         if xN[j] is on its upper bound
//    lN[j] = uN[j], if xN[j] is fixed variable
//
// where lN[j] and uN[j] are lower and upper bounds of xN[j]. */

void ssx_get_xNj(SSX *ssx, int j, mpq_t x)
{     int m = ssx->m;
      int n = ssx->n;
      mpq_t *lb = ssx->lb;
      mpq_t *ub = ssx->ub;
      int *stat = ssx->stat;
      int *Q_col = ssx->Q_col;
      int k;
      xassert(1 <= j && j <= n);
      k = Q_col[m+j]; /* x[k] = xN[j] */
      xassert(1 <= k && k <= m+n);
      switch (stat[k])
      {  case SSX_NL:
            /* xN[j] is on its lower bound */
            mpq_set(x, lb[k]); break;
         case SSX_NU:
            /* xN[j] is on its upper bound */
            mpq_set(x, ub[k]); break;
         case SSX_NF:
            /* xN[j] is free variable */
            mpq_set_si(x, 0, 1); break;
         case SSX_NS:
            /* xN[j] is fixed variable */
            mpq_set(x, lb[k]); break;
         default:
            xassert(stat != stat);
      }
      return;
}

/*----------------------------------------------------------------------
// ssx_eval_bbar - compute values of basic variables.
//
// This routine computes values of basic variables xB in the current
// basic solution as follows:
//
//    beta = - inv(B) * N * xN,
//
// where B is the basis matrix, N is the matrix of non-basic columns,
// xN is a vector of current values of non-basic variables. */

void ssx_eval_bbar(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      mpq_t *coef = ssx->coef;
      int *A_ptr = ssx->A_ptr;
      int *A_ind = ssx->A_ind;
      mpq_t *A_val = ssx->A_val;
      int *Q_col = ssx->Q_col;
      mpq_t *bbar = ssx->bbar;
      int i, j, k, ptr;
      mpq_t x, temp;
      mpq_init(x);
      mpq_init(temp);
      /* bbar := 0 */
      for (i = 1; i <= m; i++)
         mpq_set_si(bbar[i], 0, 1);
      /* bbar := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n] */
      for (j = 1; j <= n; j++)
      {  ssx_get_xNj(ssx, j, x);
         if (mpq_sgn(x) == 0) continue;
         k = Q_col[m+j]; /* x[k] = xN[j] */
         if (k <= m)
         {  /* N[j] is a column of the unity matrix I */
            mpq_sub(bbar[k], bbar[k], x);
         }
         else
         {  /* N[j] is a column of the original constraint matrix -A */
            for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
            {  mpq_mul(temp, A_val[ptr], x);
               mpq_add(bbar[A_ind[ptr]], bbar[A_ind[ptr]], temp);
            }
         }
      }
      /* bbar := inv(B) * bbar */
      bfx_ftran(ssx->binv, bbar, 0);
#if 1
      /* compute value of the objective function */
      /* bbar[0] := c[0] */
      mpq_set(bbar[0], coef[0]);
      /* bbar[0] := bbar[0] + sum{i in B} cB[i] * xB[i] */
      for (i = 1; i <= m; i++)
      {  k = Q_col[i]; /* x[k] = xB[i] */
         if (mpq_sgn(coef[k]) == 0) continue;
         mpq_mul(temp, coef[k], bbar[i]);
         mpq_add(bbar[0], bbar[0], temp);
      }
      /* bbar[0] := bbar[0] + sum{j in N} cN[j] * xN[j] */
      for (j = 1; j <= n; j++)
      {  k = Q_col[m+j]; /* x[k] = xN[j] */
         if (mpq_sgn(coef[k]) == 0) continue;
         ssx_get_xNj(ssx, j, x);
         mpq_mul(temp, coef[k], x);
         mpq_add(bbar[0], bbar[0], temp);
      }
#endif
      mpq_clear(x);
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
// ssx_eval_pi - compute values of simplex multipliers.
//
// This routine computes values of simplex multipliers (shadow prices)
// pi in the current basic solution as follows:
//
//    pi = inv(B') * cB,
//
// where B' is a matrix transposed to the basis matrix B, cB is a vector
// of objective coefficients at basic variables xB. */

void ssx_eval_pi(SSX *ssx)
{     int m = ssx->m;
      mpq_t *coef = ssx->coef;
      int *Q_col = ssx->Q_col;
      mpq_t *pi = ssx->pi;
      int i;
      /* pi := cB */
      for (i = 1; i <= m; i++) mpq_set(pi[i], coef[Q_col[i]]);
      /* pi := inv(B') * cB */
      bfx_btran(ssx->binv, pi);
      return;
}

/*----------------------------------------------------------------------
// ssx_eval_dj - compute reduced cost of non-basic variable.
//
// This routine computes reduced cost d[j] of non-basic variable xN[j]
// in the current basic solution as follows:
//
//    d[j] = cN[j] - N[j] * pi,
//
// where cN[j] is an objective coefficient at xN[j], N[j] is a column
// of the augmented constraint matrix (I | -A) corresponding to xN[j],
// pi is the vector of simplex multipliers (shadow prices). */

void ssx_eval_dj(SSX *ssx, int j, mpq_t dj)
{     int m = ssx->m;
      int n = ssx->n;
      mpq_t *coef = ssx->coef;
      int *A_ptr = ssx->A_ptr;
      int *A_ind = ssx->A_ind;
      mpq_t *A_val = ssx->A_val;
      int *Q_col = ssx->Q_col;
      mpq_t *pi = ssx->pi;
      int k, ptr, end;
      mpq_t temp;
      mpq_init(temp);
      xassert(1 <= j && j <= n);
      k = Q_col[m+j]; /* x[k] = xN[j] */
      xassert(1 <= k && k <= m+n);
      /* j-th column of the matrix N is k-th column of the augmented
         constraint matrix (I | -A) */
      if (k <= m)
      {  /* it is a column of the unity matrix I */
         mpq_sub(dj, coef[k], pi[k]);
      }
      else
      {  /* it is a column of the original constraint matrix -A */
         mpq_set(dj, coef[k]);
         for (ptr = A_ptr[k-m], end = A_ptr[k-m+1]; ptr < end; ptr++)
         {  mpq_mul(temp, A_val[ptr], pi[A_ind[ptr]]);
            mpq_add(dj, dj, temp);
         }
      }
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
// ssx_eval_cbar - compute reduced costs of all non-basic variables.
//
// This routine computes the vector of reduced costs pi in the current
// basic solution for all non-basic variables, including fixed ones. */

void ssx_eval_cbar(SSX *ssx)
{     int n = ssx->n;
      mpq_t *cbar = ssx->cbar;
      int j;
      for (j = 1; j <= n; j++)
         ssx_eval_dj(ssx, j, cbar[j]);
      return;
}

/*----------------------------------------------------------------------
// ssx_eval_rho - compute p-th row of the inverse.
//
// This routine computes p-th row of the matrix inv(B), where B is the
// current basis matrix.
//
// p-th row of the inverse is computed using the following formula:
//
//    rho = inv(B') * e[p],
//
// where B' is a matrix transposed to B, e[p] is a unity vector, which
// contains one in p-th position. */

void ssx_eval_rho(SSX *ssx)
{     int m = ssx->m;
      int p = ssx->p;
      mpq_t *rho = ssx->rho;
      int i;
      xassert(1 <= p && p <= m);
      /* rho := 0 */
      for (i = 1; i <= m; i++) mpq_set_si(rho[i], 0, 1);
      /* rho := e[p] */
      mpq_set_si(rho[p], 1, 1);
      /* rho := inv(B') * rho */
      bfx_btran(ssx->binv, rho);
      return;
}

/*----------------------------------------------------------------------
// ssx_eval_row - compute pivot row of the simplex table.
//
// This routine computes p-th (pivot) row of the current simplex table
// A~ = - inv(B) * N using the following formula:
//
//    A~[p] = - N' * inv(B') * e[p] = - N' * rho[p],
//
// where N' is a matrix transposed to the matrix N, rho[p] is p-th row
// of the inverse inv(B). */

void ssx_eval_row(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int *A_ptr = ssx->A_ptr;
      int *A_ind = ssx->A_ind;
      mpq_t *A_val = ssx->A_val;
      int *Q_col = ssx->Q_col;
      mpq_t *rho = ssx->rho;
      mpq_t *ap = ssx->ap;
      int j, k, ptr;
      mpq_t temp;
      mpq_init(temp);
      for (j = 1; j <= n; j++)
      {  /* ap[j] := - N'[j] * rho (inner product) */
         k = Q_col[m+j]; /* x[k] = xN[j] */
         if (k <= m)
            mpq_neg(ap[j], rho[k]);
         else
         {  mpq_set_si(ap[j], 0, 1);
            for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
            {  mpq_mul(temp, A_val[ptr], rho[A_ind[ptr]]);
               mpq_add(ap[j], ap[j], temp);
            }
         }
      }
      mpq_clear(temp);
      return;
}

/*----------------------------------------------------------------------
// ssx_eval_col - compute pivot column of the simplex table.
//
// This routine computes q-th (pivot) column of the current simplex
// table A~ = - inv(B) * N using the following formula:
//
//    A~[q] = - inv(B) * N[q],
//
// where N[q] is q-th column of the matrix N corresponding to chosen
// non-basic variable xN[q]. */

void ssx_eval_col(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int *A_ptr = ssx->A_ptr;
      int *A_ind = ssx->A_ind;
      mpq_t *A_val = ssx->A_val;
      int *Q_col = ssx->Q_col;
      int q = ssx->q;
      mpq_t *aq = ssx->aq;
      int i, k, ptr;
      xassert(1 <= q && q <= n);
      /* aq := 0 */
      for (i = 1; i <= m; i++) mpq_set_si(aq[i], 0, 1);
      /* aq := N[q] */