alpar@1: /* glpssx01.c */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpenv.h"
alpar@1: #include "glpssx.h"
alpar@1: #define xfault xerror
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_create - create simplex solver workspace.
alpar@1: //
alpar@1: // This routine creates the workspace used by simplex solver routines,
alpar@1: // and returns a pointer to it.
alpar@1: //
alpar@1: // Parameters m, n, and nnz specify, respectively, the number of rows,
alpar@1: // columns, and non-zero constraint coefficients.
alpar@1: //
alpar@1: // This routine only allocates the memory for the workspace components,
alpar@1: // so the workspace needs to be saturated by data. */
alpar@1: 
alpar@1: SSX *ssx_create(int m, int n, int nnz)
alpar@1: {     SSX *ssx;
alpar@1:       int i, j, k;
alpar@1:       if (m < 1)
alpar@1:          xfault("ssx_create: m = %d; invalid number of rows\n", m);
alpar@1:       if (n < 1)
alpar@1:          xfault("ssx_create: n = %d; invalid number of columns\n", n);
alpar@1:       if (nnz < 0)
alpar@1:          xfault("ssx_create: nnz = %d; invalid number of non-zero const"
alpar@1:             "raint coefficients\n", nnz);
alpar@1:       ssx = xmalloc(sizeof(SSX));
alpar@1:       ssx->m = m;
alpar@1:       ssx->n = n;
alpar@1:       ssx->type = xcalloc(1+m+n, sizeof(int));
alpar@1:       ssx->lb = xcalloc(1+m+n, sizeof(mpq_t));
alpar@1:       for (k = 1; k <= m+n; k++) mpq_init(ssx->lb[k]);
alpar@1:       ssx->ub = xcalloc(1+m+n, sizeof(mpq_t));
alpar@1:       for (k = 1; k <= m+n; k++) mpq_init(ssx->ub[k]);
alpar@1:       ssx->coef = xcalloc(1+m+n, sizeof(mpq_t));
alpar@1:       for (k = 0; k <= m+n; k++) mpq_init(ssx->coef[k]);
alpar@1:       ssx->A_ptr = xcalloc(1+n+1, sizeof(int));
alpar@1:       ssx->A_ptr[n+1] = nnz+1;
alpar@1:       ssx->A_ind = xcalloc(1+nnz, sizeof(int));
alpar@1:       ssx->A_val = xcalloc(1+nnz, sizeof(mpq_t));
alpar@1:       for (k = 1; k <= nnz; k++) mpq_init(ssx->A_val[k]);
alpar@1:       ssx->stat = xcalloc(1+m+n, sizeof(int));
alpar@1:       ssx->Q_row = xcalloc(1+m+n, sizeof(int));
alpar@1:       ssx->Q_col = xcalloc(1+m+n, sizeof(int));
alpar@1:       ssx->binv = bfx_create_binv();
alpar@1:       ssx->bbar = xcalloc(1+m, sizeof(mpq_t));
alpar@1:       for (i = 0; i <= m; i++) mpq_init(ssx->bbar[i]);
alpar@1:       ssx->pi = xcalloc(1+m, sizeof(mpq_t));
alpar@1:       for (i = 1; i <= m; i++) mpq_init(ssx->pi[i]);
alpar@1:       ssx->cbar = xcalloc(1+n, sizeof(mpq_t));
alpar@1:       for (j = 1; j <= n; j++) mpq_init(ssx->cbar[j]);
alpar@1:       ssx->rho = xcalloc(1+m, sizeof(mpq_t));
alpar@1:       for (i = 1; i <= m; i++) mpq_init(ssx->rho[i]);
alpar@1:       ssx->ap = xcalloc(1+n, sizeof(mpq_t));
alpar@1:       for (j = 1; j <= n; j++) mpq_init(ssx->ap[j]);
alpar@1:       ssx->aq = xcalloc(1+m, sizeof(mpq_t));
alpar@1:       for (i = 1; i <= m; i++) mpq_init(ssx->aq[i]);
alpar@1:       mpq_init(ssx->delta);
alpar@1:       return ssx;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_factorize - factorize the current basis matrix.
alpar@1: //
alpar@1: // This routine computes factorization of the current basis matrix B
alpar@1: // and returns the singularity flag. If the matrix B is non-singular,
alpar@1: // the flag is zero, otherwise non-zero. */
alpar@1: 
alpar@1: static int basis_col(void *info, int j, int ind[], mpq_t val[])
alpar@1: {     /* this auxiliary routine provides row indices and numeric values
alpar@1:          of non-zero elements in j-th column of the matrix B */
alpar@1:       SSX *ssx = info;
alpar@1:       int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       int *A_ptr = ssx->A_ptr;
alpar@1:       int *A_ind = ssx->A_ind;
alpar@1:       mpq_t *A_val = ssx->A_val;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       int k, len, ptr;
alpar@1:       xassert(1 <= j && j <= m);
alpar@1:       k = Q_col[j]; /* x[k] = xB[j] */
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1:       /* j-th column of the matrix B is k-th column of the augmented
alpar@1:          constraint matrix (I | -A) */
alpar@1:       if (k <= m)
alpar@1:       {  /* it is a column of the unity matrix I */
alpar@1:          len = 1, ind[1] = k, mpq_set_si(val[1], 1, 1);
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* it is a column of the original constraint matrix -A */
alpar@1:          len = 0;
alpar@1:          for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
alpar@1:          {  len++;
alpar@1:             ind[len] = A_ind[ptr];
alpar@1:             mpq_neg(val[len], A_val[ptr]);
alpar@1:          }
alpar@1:       }
alpar@1:       return len;
alpar@1: }
alpar@1: 
alpar@1: int ssx_factorize(SSX *ssx)
alpar@1: {     int ret;
alpar@1:       ret = bfx_factorize(ssx->binv, ssx->m, basis_col, ssx);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_get_xNj - determine value of non-basic variable.
alpar@1: //
alpar@1: // This routine determines the value of non-basic variable xN[j] in the
alpar@1: // current basic solution defined as follows:
alpar@1: //
alpar@1: //    0,             if xN[j] is free variable
alpar@1: //    lN[j],         if xN[j] is on its lower bound
alpar@1: //    uN[j],         if xN[j] is on its upper bound
alpar@1: //    lN[j] = uN[j], if xN[j] is fixed variable
alpar@1: //
alpar@1: // where lN[j] and uN[j] are lower and upper bounds of xN[j]. */
alpar@1: 
alpar@1: void ssx_get_xNj(SSX *ssx, int j, mpq_t x)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       mpq_t *lb = ssx->lb;
alpar@1:       mpq_t *ub = ssx->ub;
alpar@1:       int *stat = ssx->stat;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       int k;
alpar@1:       xassert(1 <= j && j <= n);
alpar@1:       k = Q_col[m+j]; /* x[k] = xN[j] */
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1:       switch (stat[k])
alpar@1:       {  case SSX_NL:
alpar@1:             /* xN[j] is on its lower bound */
alpar@1:             mpq_set(x, lb[k]); break;
alpar@1:          case SSX_NU:
alpar@1:             /* xN[j] is on its upper bound */
alpar@1:             mpq_set(x, ub[k]); break;
alpar@1:          case SSX_NF:
alpar@1:             /* xN[j] is free variable */
alpar@1:             mpq_set_si(x, 0, 1); break;
alpar@1:          case SSX_NS:
alpar@1:             /* xN[j] is fixed variable */
alpar@1:             mpq_set(x, lb[k]); break;
alpar@1:          default:
alpar@1:             xassert(stat != stat);
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_eval_bbar - compute values of basic variables.
alpar@1: //
alpar@1: // This routine computes values of basic variables xB in the current
alpar@1: // basic solution as follows:
alpar@1: //
alpar@1: //    beta = - inv(B) * N * xN,
alpar@1: //
alpar@1: // where B is the basis matrix, N is the matrix of non-basic columns,
alpar@1: // xN is a vector of current values of non-basic variables. */
alpar@1: 
alpar@1: void ssx_eval_bbar(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       mpq_t *coef = ssx->coef;
alpar@1:       int *A_ptr = ssx->A_ptr;
alpar@1:       int *A_ind = ssx->A_ind;
alpar@1:       mpq_t *A_val = ssx->A_val;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       mpq_t *bbar = ssx->bbar;
alpar@1:       int i, j, k, ptr;
alpar@1:       mpq_t x, temp;
alpar@1:       mpq_init(x);
alpar@1:       mpq_init(temp);
alpar@1:       /* bbar := 0 */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:          mpq_set_si(bbar[i], 0, 1);
alpar@1:       /* bbar := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n] */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  ssx_get_xNj(ssx, j, x);
alpar@1:          if (mpq_sgn(x) == 0) continue;
alpar@1:          k = Q_col[m+j]; /* x[k] = xN[j] */
alpar@1:          if (k <= m)
alpar@1:          {  /* N[j] is a column of the unity matrix I */
alpar@1:             mpq_sub(bbar[k], bbar[k], x);
alpar@1:          }
alpar@1:          else
alpar@1:          {  /* N[j] is a column of the original constraint matrix -A */
alpar@1:             for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
alpar@1:             {  mpq_mul(temp, A_val[ptr], x);
alpar@1:                mpq_add(bbar[A_ind[ptr]], bbar[A_ind[ptr]], temp);
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* bbar := inv(B) * bbar */
alpar@1:       bfx_ftran(ssx->binv, bbar, 0);
alpar@1: #if 1
alpar@1:       /* compute value of the objective function */
alpar@1:       /* bbar[0] := c[0] */
alpar@1:       mpq_set(bbar[0], coef[0]);
alpar@1:       /* bbar[0] := bbar[0] + sum{i in B} cB[i] * xB[i] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  k = Q_col[i]; /* x[k] = xB[i] */
alpar@1:          if (mpq_sgn(coef[k]) == 0) continue;
alpar@1:          mpq_mul(temp, coef[k], bbar[i]);
alpar@1:          mpq_add(bbar[0], bbar[0], temp);
alpar@1:       }
alpar@1:       /* bbar[0] := bbar[0] + sum{j in N} cN[j] * xN[j] */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  k = Q_col[m+j]; /* x[k] = xN[j] */
alpar@1:          if (mpq_sgn(coef[k]) == 0) continue;
alpar@1:          ssx_get_xNj(ssx, j, x);
alpar@1:          mpq_mul(temp, coef[k], x);
alpar@1:          mpq_add(bbar[0], bbar[0], temp);
alpar@1:       }
alpar@1: #endif
alpar@1:       mpq_clear(x);
alpar@1:       mpq_clear(temp);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_eval_pi - compute values of simplex multipliers.
alpar@1: //
alpar@1: // This routine computes values of simplex multipliers (shadow prices)
alpar@1: // pi in the current basic solution as follows:
alpar@1: //
alpar@1: //    pi = inv(B') * cB,
alpar@1: //
alpar@1: // where B' is a matrix transposed to the basis matrix B, cB is a vector
alpar@1: // of objective coefficients at basic variables xB. */
alpar@1: 
alpar@1: void ssx_eval_pi(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       mpq_t *coef = ssx->coef;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       mpq_t *pi = ssx->pi;
alpar@1:       int i;
alpar@1:       /* pi := cB */
alpar@1:       for (i = 1; i <= m; i++) mpq_set(pi[i], coef[Q_col[i]]);
alpar@1:       /* pi := inv(B') * cB */
alpar@1:       bfx_btran(ssx->binv, pi);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_eval_dj - compute reduced cost of non-basic variable.
alpar@1: //
alpar@1: // This routine computes reduced cost d[j] of non-basic variable xN[j]
alpar@1: // in the current basic solution as follows:
alpar@1: //
alpar@1: //    d[j] = cN[j] - N[j] * pi,
alpar@1: //
alpar@1: // where cN[j] is an objective coefficient at xN[j], N[j] is a column
alpar@1: // of the augmented constraint matrix (I | -A) corresponding to xN[j],
alpar@1: // pi is the vector of simplex multipliers (shadow prices). */
alpar@1: 
alpar@1: void ssx_eval_dj(SSX *ssx, int j, mpq_t dj)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       mpq_t *coef = ssx->coef;
alpar@1:       int *A_ptr = ssx->A_ptr;
alpar@1:       int *A_ind = ssx->A_ind;
alpar@1:       mpq_t *A_val = ssx->A_val;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       mpq_t *pi = ssx->pi;
alpar@1:       int k, ptr, end;
alpar@1:       mpq_t temp;
alpar@1:       mpq_init(temp);
alpar@1:       xassert(1 <= j && j <= n);
alpar@1:       k = Q_col[m+j]; /* x[k] = xN[j] */
alpar@1:       xassert(1 <= k && k <= m+n);
alpar@1:       /* j-th column of the matrix N is k-th column of the augmented
alpar@1:          constraint matrix (I | -A) */
alpar@1:       if (k <= m)
alpar@1:       {  /* it is a column of the unity matrix I */
alpar@1:          mpq_sub(dj, coef[k], pi[k]);
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* it is a column of the original constraint matrix -A */
alpar@1:          mpq_set(dj, coef[k]);
alpar@1:          for (ptr = A_ptr[k-m], end = A_ptr[k-m+1]; ptr < end; ptr++)
alpar@1:          {  mpq_mul(temp, A_val[ptr], pi[A_ind[ptr]]);
alpar@1:             mpq_add(dj, dj, temp);
alpar@1:          }
alpar@1:       }
alpar@1:       mpq_clear(temp);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_eval_cbar - compute reduced costs of all non-basic variables.
alpar@1: //
alpar@1: // This routine computes the vector of reduced costs pi in the current
alpar@1: // basic solution for all non-basic variables, including fixed ones. */
alpar@1: 
alpar@1: void ssx_eval_cbar(SSX *ssx)
alpar@1: {     int n = ssx->n;
alpar@1:       mpq_t *cbar = ssx->cbar;
alpar@1:       int j;
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:          ssx_eval_dj(ssx, j, cbar[j]);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_eval_rho - compute p-th row of the inverse.
alpar@1: //
alpar@1: // This routine computes p-th row of the matrix inv(B), where B is the
alpar@1: // current basis matrix.
alpar@1: //
alpar@1: // p-th row of the inverse is computed using the following formula:
alpar@1: //
alpar@1: //    rho = inv(B') * e[p],
alpar@1: //
alpar@1: // where B' is a matrix transposed to B, e[p] is a unity vector, which
alpar@1: // contains one in p-th position. */
alpar@1: 
alpar@1: void ssx_eval_rho(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int p = ssx->p;
alpar@1:       mpq_t *rho = ssx->rho;
alpar@1:       int i;
alpar@1:       xassert(1 <= p && p <= m);
alpar@1:       /* rho := 0 */
alpar@1:       for (i = 1; i <= m; i++) mpq_set_si(rho[i], 0, 1);
alpar@1:       /* rho := e[p] */
alpar@1:       mpq_set_si(rho[p], 1, 1);
alpar@1:       /* rho := inv(B') * rho */
alpar@1:       bfx_btran(ssx->binv, rho);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_eval_row - compute pivot row of the simplex table.
alpar@1: //
alpar@1: // This routine computes p-th (pivot) row of the current simplex table
alpar@1: // A~ = - inv(B) * N using the following formula:
alpar@1: //
alpar@1: //    A~[p] = - N' * inv(B') * e[p] = - N' * rho[p],
alpar@1: //
alpar@1: // where N' is a matrix transposed to the matrix N, rho[p] is p-th row
alpar@1: // of the inverse inv(B). */
alpar@1: 
alpar@1: void ssx_eval_row(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       int *A_ptr = ssx->A_ptr;
alpar@1:       int *A_ind = ssx->A_ind;
alpar@1:       mpq_t *A_val = ssx->A_val;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       mpq_t *rho = ssx->rho;
alpar@1:       mpq_t *ap = ssx->ap;
alpar@1:       int j, k, ptr;
alpar@1:       mpq_t temp;
alpar@1:       mpq_init(temp);
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  /* ap[j] := - N'[j] * rho (inner product) */
alpar@1:          k = Q_col[m+j]; /* x[k] = xN[j] */
alpar@1:          if (k <= m)
alpar@1:             mpq_neg(ap[j], rho[k]);
alpar@1:          else
alpar@1:          {  mpq_set_si(ap[j], 0, 1);
alpar@1:             for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
alpar@1:             {  mpq_mul(temp, A_val[ptr], rho[A_ind[ptr]]);
alpar@1:                mpq_add(ap[j], ap[j], temp);
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       mpq_clear(temp);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_eval_col - compute pivot column of the simplex table.
alpar@1: //
alpar@1: // This routine computes q-th (pivot) column of the current simplex
alpar@1: // table A~ = - inv(B) * N using the following formula:
alpar@1: //
alpar@1: //    A~[q] = - inv(B) * N[q],
alpar@1: //
alpar@1: // where N[q] is q-th column of the matrix N corresponding to chosen
alpar@1: // non-basic variable xN[q]. */
alpar@1: 
alpar@1: void ssx_eval_col(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       int *A_ptr = ssx->A_ptr;
alpar@1:       int *A_ind = ssx->A_ind;
alpar@1:       mpq_t *A_val = ssx->A_val;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       int q = ssx->q;
alpar@1:       mpq_t *aq = ssx->aq;
alpar@1:       int i, k, ptr;
alpar@1:       xassert(1 <= q && q <= n);
alpar@1:       /* aq := 0 */
alpar@1:       for (i = 1; i <= m; i++) mpq_set_si(aq[i], 0, 1);
alpar@1:       /* aq := N[q] */
alpar@1:       k = Q_col[m+q]; /* x[k] = xN[q] */
alpar@1:       if (k <= m)
alpar@1:       {  /* N[q] is a column of the unity matrix I */
alpar@1:          mpq_set_si(aq[k], 1, 1);
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* N[q] is a column of the original constraint matrix -A */
alpar@1:          for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
alpar@1:             mpq_neg(aq[A_ind[ptr]], A_val[ptr]);
alpar@1:       }
alpar@1:       /* aq := inv(B) * aq */
alpar@1:       bfx_ftran(ssx->binv, aq, 1);
alpar@1:       /* aq := - aq */
alpar@1:       for (i = 1; i <= m; i++) mpq_neg(aq[i], aq[i]);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_chuzc - choose pivot column.
alpar@1: //
alpar@1: // This routine chooses non-basic variable xN[q] whose reduced cost
alpar@1: // indicates possible improving of the objective function to enter it
alpar@1: // in the basis.
alpar@1: //
alpar@1: // Currently the standard (textbook) pricing is used, i.e. that
alpar@1: // non-basic variable is preferred which has greatest reduced cost (in
alpar@1: // magnitude).
alpar@1: //
alpar@1: // If xN[q] has been chosen, the routine stores its number q and also
alpar@1: // sets the flag q_dir that indicates direction in which xN[q] has to
alpar@1: // change (+1 means increasing, -1 means decreasing).
alpar@1: //
alpar@1: // If the choice cannot be made, because the current basic solution is
alpar@1: // dual feasible, the routine sets the number q to 0. */
alpar@1: 
alpar@1: void ssx_chuzc(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       int dir = (ssx->dir == SSX_MIN ? +1 : -1);
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       int *stat = ssx->stat;
alpar@1:       mpq_t *cbar = ssx->cbar;
alpar@1:       int j, k, s, q, q_dir;
alpar@1:       double best, temp;
alpar@1:       /* nothing is chosen so far */
alpar@1:       q = 0, q_dir = 0, best = 0.0;
alpar@1:       /* look through the list of non-basic variables */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  k = Q_col[m+j]; /* x[k] = xN[j] */
alpar@1:          s = dir * mpq_sgn(cbar[j]);
alpar@1:          if ((stat[k] == SSX_NF || stat[k] == SSX_NL) && s < 0 ||
alpar@1:              (stat[k] == SSX_NF || stat[k] == SSX_NU) && s > 0)
alpar@1:          {  /* reduced cost of xN[j] indicates possible improving of
alpar@1:                the objective function */
alpar@1:             temp = fabs(mpq_get_d(cbar[j]));
alpar@1:             xassert(temp != 0.0);
alpar@1:             if (q == 0 || best < temp)
alpar@1:                q = j, q_dir = - s, best = temp;
alpar@1:          }
alpar@1:       }
alpar@1:       ssx->q = q, ssx->q_dir = q_dir;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_chuzr - choose pivot row.
alpar@1: //
alpar@1: // This routine looks through elements of q-th column of the simplex
alpar@1: // table and chooses basic variable xB[p] which should leave the basis.
alpar@1: //
alpar@1: // The choice is based on the standard (textbook) ratio test.
alpar@1: //
alpar@1: // If xB[p] has been chosen, the routine stores its number p and also
alpar@1: // sets its non-basic status p_stat which should be assigned to xB[p]
alpar@1: // when it has left the basis and become xN[q].
alpar@1: //
alpar@1: // Special case p < 0 means that xN[q] is double-bounded variable and
alpar@1: // it reaches its opposite bound before any basic variable does that,
alpar@1: // so the current basis remains unchanged.
alpar@1: //
alpar@1: // If the choice cannot be made, because xN[q] can infinitely change in
alpar@1: // the feasible direction, the routine sets the number p to 0. */
alpar@1: 
alpar@1: void ssx_chuzr(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       int *type = ssx->type;
alpar@1:       mpq_t *lb = ssx->lb;
alpar@1:       mpq_t *ub = ssx->ub;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       mpq_t *bbar = ssx->bbar;
alpar@1:       int q = ssx->q;
alpar@1:       mpq_t *aq = ssx->aq;
alpar@1:       int q_dir = ssx->q_dir;
alpar@1:       int i, k, s, t, p, p_stat;
alpar@1:       mpq_t teta, temp;
alpar@1:       mpq_init(teta);
alpar@1:       mpq_init(temp);
alpar@1:       xassert(1 <= q && q <= n);
alpar@1:       xassert(q_dir == +1 || q_dir == -1);
alpar@1:       /* nothing is chosen so far */
alpar@1:       p = 0, p_stat = 0;
alpar@1:       /* look through the list of basic variables */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  s = q_dir * mpq_sgn(aq[i]);
alpar@1:          if (s < 0)
alpar@1:          {  /* xB[i] decreases */
alpar@1:             k = Q_col[i]; /* x[k] = xB[i] */
alpar@1:             t = type[k];
alpar@1:             if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
alpar@1:             {  /* xB[i] has finite lower bound */
alpar@1:                mpq_sub(temp, bbar[i], lb[k]);
alpar@1:                mpq_div(temp, temp, aq[i]);
alpar@1:                mpq_abs(temp, temp);
alpar@1:                if (p == 0 || mpq_cmp(teta, temp) > 0)
alpar@1:                {  p = i;
alpar@1:                   p_stat = (t == SSX_FX ? SSX_NS : SSX_NL);
alpar@1:                   mpq_set(teta, temp);
alpar@1:                }
alpar@1:             }
alpar@1:          }
alpar@1:          else if (s > 0)
alpar@1:          {  /* xB[i] increases */
alpar@1:             k = Q_col[i]; /* x[k] = xB[i] */
alpar@1:             t = type[k];
alpar@1:             if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
alpar@1:             {  /* xB[i] has finite upper bound */
alpar@1:                mpq_sub(temp, bbar[i], ub[k]);
alpar@1:                mpq_div(temp, temp, aq[i]);
alpar@1:                mpq_abs(temp, temp);
alpar@1:                if (p == 0 || mpq_cmp(teta, temp) > 0)
alpar@1:                {  p = i;
alpar@1:                   p_stat = (t == SSX_FX ? SSX_NS : SSX_NU);
alpar@1:                   mpq_set(teta, temp);
alpar@1:                }
alpar@1:             }
alpar@1:          }
alpar@1:          /* if something has been chosen and the ratio test indicates
alpar@1:             exact degeneracy, the search can be finished */
alpar@1:          if (p != 0 && mpq_sgn(teta) == 0) break;
alpar@1:       }
alpar@1:       /* if xN[q] is double-bounded, check if it can reach its opposite
alpar@1:          bound before any basic variable */
alpar@1:       k = Q_col[m+q]; /* x[k] = xN[q] */
alpar@1:       if (type[k] == SSX_DB)
alpar@1:       {  mpq_sub(temp, ub[k], lb[k]);
alpar@1:          if (p == 0 || mpq_cmp(teta, temp) > 0)
alpar@1:          {  p = -1;
alpar@1:             p_stat = -1;
alpar@1:             mpq_set(teta, temp);
alpar@1:          }
alpar@1:       }
alpar@1:       ssx->p = p;
alpar@1:       ssx->p_stat = p_stat;
alpar@1:       /* if xB[p] has been chosen, determine its actual change in the
alpar@1:          adjacent basis (it has the same sign as q_dir) */
alpar@1:       if (p != 0)
alpar@1:       {  xassert(mpq_sgn(teta) >= 0);
alpar@1:          if (q_dir > 0)
alpar@1:             mpq_set(ssx->delta, teta);
alpar@1:          else
alpar@1:             mpq_neg(ssx->delta, teta);
alpar@1:       }
alpar@1:       mpq_clear(teta);
alpar@1:       mpq_clear(temp);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_update_bbar - update values of basic variables.
alpar@1: //
alpar@1: // This routine recomputes the current values of basic variables for
alpar@1: // the adjacent basis.
alpar@1: //
alpar@1: // The simplex table for the current basis is the following:
alpar@1: //
alpar@1: //    xB[i] = sum{j in 1..n} alfa[i,j] * xN[q],  i = 1,...,m
alpar@1: //
alpar@1: // therefore
alpar@1: //
alpar@1: //    delta xB[i] = alfa[i,q] * delta xN[q],  i = 1,...,m
alpar@1: //
alpar@1: // where delta xN[q] = xN.new[q] - xN[q] is the change of xN[q] in the
alpar@1: // adjacent basis, and delta xB[i] = xB.new[i] - xB[i] is the change of
alpar@1: // xB[i]. This gives formulae for recomputing values of xB[i]:
alpar@1: //
alpar@1: //    xB.new[p] = xN[q] + delta xN[q]
alpar@1: //
alpar@1: // (because xN[q] becomes xB[p] in the adjacent basis), and
alpar@1: //
alpar@1: //    xB.new[i] = xB[i] + alfa[i,q] * delta xN[q],  i != p
alpar@1: //
alpar@1: // for other basic variables. */
alpar@1: 
alpar@1: void ssx_update_bbar(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       mpq_t *bbar = ssx->bbar;
alpar@1:       mpq_t *cbar = ssx->cbar;
alpar@1:       int p = ssx->p;
alpar@1:       int q = ssx->q;
alpar@1:       mpq_t *aq = ssx->aq;
alpar@1:       int i;
alpar@1:       mpq_t temp;
alpar@1:       mpq_init(temp);
alpar@1:       xassert(1 <= q && q <= n);
alpar@1:       if (p < 0)
alpar@1:       {  /* xN[q] is double-bounded and goes to its opposite bound */
alpar@1:          /* nop */;
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* xN[q] becomes xB[p] in the adjacent basis */
alpar@1:          /* xB.new[p] = xN[q] + delta xN[q] */
alpar@1:          xassert(1 <= p && p <= m);
alpar@1:          ssx_get_xNj(ssx, q, temp);
alpar@1:          mpq_add(bbar[p], temp, ssx->delta);
alpar@1:       }
alpar@1:       /* update values of other basic variables depending on xN[q] */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  if (i == p) continue;
alpar@1:          /* xB.new[i] = xB[i] + alfa[i,q] * delta xN[q] */
alpar@1:          if (mpq_sgn(aq[i]) == 0) continue;
alpar@1:          mpq_mul(temp, aq[i], ssx->delta);
alpar@1:          mpq_add(bbar[i], bbar[i], temp);
alpar@1:       }
alpar@1: #if 1
alpar@1:       /* update value of the objective function */
alpar@1:       /* z.new = z + d[q] * delta xN[q] */
alpar@1:       mpq_mul(temp, cbar[q], ssx->delta);
alpar@1:       mpq_add(bbar[0], bbar[0], temp);
alpar@1: #endif
alpar@1:       mpq_clear(temp);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: -- ssx_update_pi - update simplex multipliers.
alpar@1: --
alpar@1: -- This routine recomputes the vector of simplex multipliers for the
alpar@1: -- adjacent basis. */
alpar@1: 
alpar@1: void ssx_update_pi(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       mpq_t *pi = ssx->pi;
alpar@1:       mpq_t *cbar = ssx->cbar;
alpar@1:       int p = ssx->p;
alpar@1:       int q = ssx->q;
alpar@1:       mpq_t *aq = ssx->aq;
alpar@1:       mpq_t *rho = ssx->rho;
alpar@1:       int i;
alpar@1:       mpq_t new_dq, temp;
alpar@1:       mpq_init(new_dq);
alpar@1:       mpq_init(temp);
alpar@1:       xassert(1 <= p && p <= m);
alpar@1:       xassert(1 <= q && q <= n);
alpar@1:       /* compute d[q] in the adjacent basis */
alpar@1:       mpq_div(new_dq, cbar[q], aq[p]);
alpar@1:       /* update the vector of simplex multipliers */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  if (mpq_sgn(rho[i]) == 0) continue;
alpar@1:          mpq_mul(temp, new_dq, rho[i]);
alpar@1:          mpq_sub(pi[i], pi[i], temp);
alpar@1:       }
alpar@1:       mpq_clear(new_dq);
alpar@1:       mpq_clear(temp);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_update_cbar - update reduced costs of non-basic variables.
alpar@1: //
alpar@1: // This routine recomputes the vector of reduced costs of non-basic
alpar@1: // variables for the adjacent basis. */
alpar@1: 
alpar@1: void ssx_update_cbar(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       mpq_t *cbar = ssx->cbar;
alpar@1:       int p = ssx->p;
alpar@1:       int q = ssx->q;
alpar@1:       mpq_t *ap = ssx->ap;
alpar@1:       int j;
alpar@1:       mpq_t temp;
alpar@1:       mpq_init(temp);
alpar@1:       xassert(1 <= p && p <= m);
alpar@1:       xassert(1 <= q && q <= n);
alpar@1:       /* compute d[q] in the adjacent basis */
alpar@1:       /* d.new[q] = d[q] / alfa[p,q] */
alpar@1:       mpq_div(cbar[q], cbar[q], ap[q]);
alpar@1:       /* update reduced costs of other non-basic variables */
alpar@1:       for (j = 1; j <= n; j++)
alpar@1:       {  if (j == q) continue;
alpar@1:          /* d.new[j] = d[j] - (alfa[p,j] / alfa[p,q]) * d[q] */
alpar@1:          if (mpq_sgn(ap[j]) == 0) continue;
alpar@1:          mpq_mul(temp, ap[j], cbar[q]);
alpar@1:          mpq_sub(cbar[j], cbar[j], temp);
alpar@1:       }
alpar@1:       mpq_clear(temp);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_change_basis - change current basis to adjacent one.
alpar@1: //
alpar@1: // This routine changes the current basis to the adjacent one swapping
alpar@1: // basic variable xB[p] and non-basic variable xN[q]. */
alpar@1: 
alpar@1: void ssx_change_basis(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       int *type = ssx->type;
alpar@1:       int *stat = ssx->stat;
alpar@1:       int *Q_row = ssx->Q_row;
alpar@1:       int *Q_col = ssx->Q_col;
alpar@1:       int p = ssx->p;
alpar@1:       int q = ssx->q;
alpar@1:       int p_stat = ssx->p_stat;
alpar@1:       int k, kp, kq;
alpar@1:       if (p < 0)
alpar@1:       {  /* special case: xN[q] goes to its opposite bound */
alpar@1:          xassert(1 <= q && q <= n);
alpar@1:          k = Q_col[m+q]; /* x[k] = xN[q] */
alpar@1:          xassert(type[k] == SSX_DB);
alpar@1:          switch (stat[k])
alpar@1:          {  case SSX_NL:
alpar@1:                stat[k] = SSX_NU;
alpar@1:                break;
alpar@1:             case SSX_NU:
alpar@1:                stat[k] = SSX_NL;
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(stat != stat);
alpar@1:          }
alpar@1:       }
alpar@1:       else
alpar@1:       {  /* xB[p] leaves the basis, xN[q] enters the basis */
alpar@1:          xassert(1 <= p && p <= m);
alpar@1:          xassert(1 <= q && q <= n);
alpar@1:          kp = Q_col[p];   /* x[kp] = xB[p] */
alpar@1:          kq = Q_col[m+q]; /* x[kq] = xN[q] */
alpar@1:          /* check non-basic status of xB[p] which becomes xN[q] */
alpar@1:          switch (type[kp])
alpar@1:          {  case SSX_FR:
alpar@1:                xassert(p_stat == SSX_NF);
alpar@1:                break;
alpar@1:             case SSX_LO:
alpar@1:                xassert(p_stat == SSX_NL);
alpar@1:                break;
alpar@1:             case SSX_UP:
alpar@1:                xassert(p_stat == SSX_NU);
alpar@1:                break;
alpar@1:             case SSX_DB:
alpar@1:                xassert(p_stat == SSX_NL || p_stat == SSX_NU);
alpar@1:                break;
alpar@1:             case SSX_FX:
alpar@1:                xassert(p_stat == SSX_NS);
alpar@1:                break;
alpar@1:             default:
alpar@1:                xassert(type != type);
alpar@1:          }
alpar@1:          /* swap xB[p] and xN[q] */
alpar@1:          stat[kp] = (char)p_stat, stat[kq] = SSX_BS;
alpar@1:          Q_row[kp] = m+q, Q_row[kq] = p;
alpar@1:          Q_col[p] = kq, Q_col[m+q] = kp;
alpar@1:          /* update factorization of the basis matrix */
alpar@1:          if (bfx_update(ssx->binv, p))
alpar@1:          {  if (ssx_factorize(ssx))
alpar@1:                xassert(("Internal error: basis matrix is singular", 0));
alpar@1:          }
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_delete - delete simplex solver workspace.
alpar@1: //
alpar@1: // This routine deletes the simplex solver workspace freeing all the
alpar@1: // memory allocated to this object. */
alpar@1: 
alpar@1: void ssx_delete(SSX *ssx)
alpar@1: {     int m = ssx->m;
alpar@1:       int n = ssx->n;
alpar@1:       int nnz = ssx->A_ptr[n+1]-1;
alpar@1:       int i, j, k;
alpar@1:       xfree(ssx->type);
alpar@1:       for (k = 1; k <= m+n; k++) mpq_clear(ssx->lb[k]);
alpar@1:       xfree(ssx->lb);
alpar@1:       for (k = 1; k <= m+n; k++) mpq_clear(ssx->ub[k]);
alpar@1:       xfree(ssx->ub);
alpar@1:       for (k = 0; k <= m+n; k++) mpq_clear(ssx->coef[k]);
alpar@1:       xfree(ssx->coef);
alpar@1:       xfree(ssx->A_ptr);
alpar@1:       xfree(ssx->A_ind);
alpar@1:       for (k = 1; k <= nnz; k++) mpq_clear(ssx->A_val[k]);
alpar@1:       xfree(ssx->A_val);
alpar@1:       xfree(ssx->stat);
alpar@1:       xfree(ssx->Q_row);
alpar@1:       xfree(ssx->Q_col);
alpar@1:       bfx_delete_binv(ssx->binv);
alpar@1:       for (i = 0; i <= m; i++) mpq_clear(ssx->bbar[i]);
alpar@1:       xfree(ssx->bbar);
alpar@1:       for (i = 1; i <= m; i++) mpq_clear(ssx->pi[i]);
alpar@1:       xfree(ssx->pi);
alpar@1:       for (j = 1; j <= n; j++) mpq_clear(ssx->cbar[j]);
alpar@1:       xfree(ssx->cbar);
alpar@1:       for (i = 1; i <= m; i++) mpq_clear(ssx->rho[i]);
alpar@1:       xfree(ssx->rho);
alpar@1:       for (j = 1; j <= n; j++) mpq_clear(ssx->ap[j]);
alpar@1:       xfree(ssx->ap);
alpar@1:       for (i = 1; i <= m; i++) mpq_clear(ssx->aq[i]);
alpar@1:       xfree(ssx->aq);
alpar@1:       mpq_clear(ssx->delta);
alpar@1:       xfree(ssx);
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /* eof */