alpar@1: /* glpssx02.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: 
alpar@1: static void show_progress(SSX *ssx, int phase)
alpar@1: {     /* this auxiliary routine displays information about progress of
alpar@1:          the search */
alpar@1:       int i, def = 0;
alpar@1:       for (i = 1; i <= ssx->m; i++)
alpar@1:          if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
alpar@1:       xprintf("%s%6d:   %s = %22.15g   (%d)\n", phase == 1 ? " " : "*",
alpar@1:          ssx->it_cnt, phase == 1 ? "infsum" : "objval",
alpar@1:          mpq_get_d(ssx->bbar[0]), def);
alpar@1: #if 0
alpar@1:       ssx->tm_lag = utime();
alpar@1: #else
alpar@1:       ssx->tm_lag = xtime();
alpar@1: #endif
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_phase_I - find primal feasible solution.
alpar@1: //
alpar@1: // This routine implements phase I of the primal simplex method.
alpar@1: //
alpar@1: // On exit the routine returns one of the following codes:
alpar@1: //
alpar@1: // 0 - feasible solution found;
alpar@1: // 1 - problem has no feasible solution;
alpar@1: // 2 - iterations limit exceeded;
alpar@1: // 3 - time limit exceeded.
alpar@1: ----------------------------------------------------------------------*/
alpar@1: 
alpar@1: int ssx_phase_I(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:       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:       mpq_t *pi = ssx->pi;
alpar@1:       mpq_t *cbar = ssx->cbar;
alpar@1:       int *orig_type, orig_dir;
alpar@1:       mpq_t *orig_lb, *orig_ub, *orig_coef;
alpar@1:       int i, k, ret;
alpar@1:       /* save components of the original LP problem, which are changed
alpar@1:          by the routine */
alpar@1:       orig_type = xcalloc(1+m+n, sizeof(int));
alpar@1:       orig_lb = xcalloc(1+m+n, sizeof(mpq_t));
alpar@1:       orig_ub = xcalloc(1+m+n, sizeof(mpq_t));
alpar@1:       orig_coef = xcalloc(1+m+n, sizeof(mpq_t));
alpar@1:       for (k = 1; k <= m+n; k++)
alpar@1:       {  orig_type[k] = type[k];
alpar@1:          mpq_init(orig_lb[k]);
alpar@1:          mpq_set(orig_lb[k], lb[k]);
alpar@1:          mpq_init(orig_ub[k]);
alpar@1:          mpq_set(orig_ub[k], ub[k]);
alpar@1:       }
alpar@1:       orig_dir = ssx->dir;
alpar@1:       for (k = 0; k <= m+n; k++)
alpar@1:       {  mpq_init(orig_coef[k]);
alpar@1:          mpq_set(orig_coef[k], coef[k]);
alpar@1:       }
alpar@1:       /* build an artificial basic solution, which is primal feasible,
alpar@1:          and also build an auxiliary objective function to minimize the
alpar@1:          sum of infeasibilities for the original problem */
alpar@1:       ssx->dir = SSX_MIN;
alpar@1:       for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1);
alpar@1:       mpq_set_si(bbar[0], 0, 1);
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  int t;
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:          {  /* in the original problem x[k] has lower bound */
alpar@1:             if (mpq_cmp(bbar[i], lb[k]) < 0)
alpar@1:             {  /* which is violated */
alpar@1:                type[k] = SSX_UP;
alpar@1:                mpq_set(ub[k], lb[k]);
alpar@1:                mpq_set_si(lb[k], 0, 1);
alpar@1:                mpq_set_si(coef[k], -1, 1);
alpar@1:                mpq_add(bbar[0], bbar[0], ub[k]);
alpar@1:                mpq_sub(bbar[0], bbar[0], bbar[i]);
alpar@1:             }
alpar@1:          }
alpar@1:          if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
alpar@1:          {  /* in the original problem x[k] has upper bound */
alpar@1:             if (mpq_cmp(bbar[i], ub[k]) > 0)
alpar@1:             {  /* which is violated */
alpar@1:                type[k] = SSX_LO;
alpar@1:                mpq_set(lb[k], ub[k]);
alpar@1:                mpq_set_si(ub[k], 0, 1);
alpar@1:                mpq_set_si(coef[k], +1, 1);
alpar@1:                mpq_add(bbar[0], bbar[0], bbar[i]);
alpar@1:                mpq_sub(bbar[0], bbar[0], lb[k]);
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       /* now the initial basic solution should be primal feasible due
alpar@1:          to changes of bounds of some basic variables, which turned to
alpar@1:          implicit artifical variables */
alpar@1:       /* compute simplex multipliers and reduced costs */
alpar@1:       ssx_eval_pi(ssx);
alpar@1:       ssx_eval_cbar(ssx);
alpar@1:       /* display initial progress of the search */
alpar@1:       show_progress(ssx, 1);
alpar@1:       /* main loop starts here */
alpar@1:       for (;;)
alpar@1:       {  /* display current progress of the search */
alpar@1: #if 0
alpar@1:          if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
alpar@1: #else
alpar@1:          if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
alpar@1: #endif
alpar@1:             show_progress(ssx, 1);
alpar@1:          /* we do not need to wait until all artificial variables have
alpar@1:             left the basis */
alpar@1:          if (mpq_sgn(bbar[0]) == 0)
alpar@1:          {  /* the sum of infeasibilities is zero, therefore the current
alpar@1:                solution is primal feasible for the original problem */
alpar@1:             ret = 0;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* check if the iterations limit has been exhausted */
alpar@1:          if (ssx->it_lim == 0)
alpar@1:          {  ret = 2;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* check if the time limit has been exhausted */
alpar@1: #if 0
alpar@1:          if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
alpar@1: #else
alpar@1:          if (ssx->tm_lim >= 0.0 &&
alpar@1:              ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
alpar@1: #endif
alpar@1:          {  ret = 3;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* choose non-basic variable xN[q] */
alpar@1:          ssx_chuzc(ssx);
alpar@1:          /* if xN[q] cannot be chosen, the sum of infeasibilities is
alpar@1:             minimal but non-zero; therefore the original problem has no
alpar@1:             primal feasible solution */
alpar@1:          if (ssx->q == 0)
alpar@1:          {  ret = 1;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* compute q-th column of the simplex table */
alpar@1:          ssx_eval_col(ssx);
alpar@1:          /* choose basic variable xB[p] */
alpar@1:          ssx_chuzr(ssx);
alpar@1:          /* the sum of infeasibilities cannot be negative, therefore
alpar@1:             the auxiliary lp problem cannot have unbounded solution */
alpar@1:          xassert(ssx->p != 0);
alpar@1:          /* update values of basic variables */
alpar@1:          ssx_update_bbar(ssx);
alpar@1:          if (ssx->p > 0)
alpar@1:          {  /* compute p-th row of the inverse inv(B) */
alpar@1:             ssx_eval_rho(ssx);
alpar@1:             /* compute p-th row of the simplex table */
alpar@1:             ssx_eval_row(ssx);
alpar@1:             xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
alpar@1:             /* update simplex multipliers */
alpar@1:             ssx_update_pi(ssx);
alpar@1:             /* update reduced costs of non-basic variables */
alpar@1:             ssx_update_cbar(ssx);
alpar@1:          }
alpar@1:          /* xB[p] is leaving the basis; if it is implicit artificial
alpar@1:             variable, the corresponding residual vanishes; therefore
alpar@1:             bounds of this variable should be restored to the original
alpar@1:             values */
alpar@1:          if (ssx->p > 0)
alpar@1:          {  k = Q_col[ssx->p]; /* x[k] = xB[p] */
alpar@1:             if (type[k] != orig_type[k])
alpar@1:             {  /* x[k] is implicit artificial variable */
alpar@1:                type[k] = orig_type[k];
alpar@1:                mpq_set(lb[k], orig_lb[k]);
alpar@1:                mpq_set(ub[k], orig_ub[k]);
alpar@1:                xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU);
alpar@1:                ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL);
alpar@1:                if (type[k] == SSX_FX) ssx->p_stat = SSX_NS;
alpar@1:                /* nullify the objective coefficient at x[k] */
alpar@1:                mpq_set_si(coef[k], 0, 1);
alpar@1:                /* since coef[k] has been changed, we need to compute
alpar@1:                   new reduced cost of x[k], which it will have in the
alpar@1:                   adjacent basis */
alpar@1:                /* the formula d[j] = cN[j] - pi' * N[j] is used (note
alpar@1:                   that the vector pi is not changed, because it depends
alpar@1:                   on objective coefficients at basic variables, but in
alpar@1:                   the adjacent basis, for which the vector pi has been
alpar@1:                   just recomputed, x[k] is non-basic) */
alpar@1:                if (k <= m)
alpar@1:                {  /* x[k] is auxiliary variable */
alpar@1:                   mpq_neg(cbar[ssx->q], pi[k]);
alpar@1:                }
alpar@1:                else
alpar@1:                {  /* x[k] is structural variable */
alpar@1:                   int ptr;
alpar@1:                   mpq_t temp;
alpar@1:                   mpq_init(temp);
alpar@1:                   mpq_set_si(cbar[ssx->q], 0, 1);
alpar@1:                   for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
alpar@1:                   {  mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]);
alpar@1:                      mpq_add(cbar[ssx->q], cbar[ssx->q], temp);
alpar@1:                   }
alpar@1:                   mpq_clear(temp);
alpar@1:                }
alpar@1:             }
alpar@1:          }
alpar@1:          /* jump to the adjacent vertex of the polyhedron */
alpar@1:          ssx_change_basis(ssx);
alpar@1:          /* one simplex iteration has been performed */
alpar@1:          if (ssx->it_lim > 0) ssx->it_lim--;
alpar@1:          ssx->it_cnt++;
alpar@1:       }
alpar@1:       /* display final progress of the search */
alpar@1:       show_progress(ssx, 1);
alpar@1:       /* restore components of the original problem, which were changed
alpar@1:          by the routine */
alpar@1:       for (k = 1; k <= m+n; k++)
alpar@1:       {  type[k] = orig_type[k];
alpar@1:          mpq_set(lb[k], orig_lb[k]);
alpar@1:          mpq_clear(orig_lb[k]);
alpar@1:          mpq_set(ub[k], orig_ub[k]);
alpar@1:          mpq_clear(orig_ub[k]);
alpar@1:       }
alpar@1:       ssx->dir = orig_dir;
alpar@1:       for (k = 0; k <= m+n; k++)
alpar@1:       {  mpq_set(coef[k], orig_coef[k]);
alpar@1:          mpq_clear(orig_coef[k]);
alpar@1:       }
alpar@1:       xfree(orig_type);
alpar@1:       xfree(orig_lb);
alpar@1:       xfree(orig_ub);
alpar@1:       xfree(orig_coef);
alpar@1:       /* return to the calling program */
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_phase_II - find optimal solution.
alpar@1: //
alpar@1: // This routine implements phase II of the primal simplex method.
alpar@1: //
alpar@1: // On exit the routine returns one of the following codes:
alpar@1: //
alpar@1: // 0 - optimal solution found;
alpar@1: // 1 - problem has unbounded solution;
alpar@1: // 2 - iterations limit exceeded;
alpar@1: // 3 - time limit exceeded.
alpar@1: ----------------------------------------------------------------------*/
alpar@1: 
alpar@1: int ssx_phase_II(SSX *ssx)
alpar@1: {     int ret;
alpar@1:       /* display initial progress of the search */
alpar@1:       show_progress(ssx, 2);
alpar@1:       /* main loop starts here */
alpar@1:       for (;;)
alpar@1:       {  /* display current progress of the search */
alpar@1: #if 0
alpar@1:          if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
alpar@1: #else
alpar@1:          if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
alpar@1: #endif
alpar@1:             show_progress(ssx, 2);
alpar@1:          /* check if the iterations limit has been exhausted */
alpar@1:          if (ssx->it_lim == 0)
alpar@1:          {  ret = 2;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* check if the time limit has been exhausted */
alpar@1: #if 0
alpar@1:          if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
alpar@1: #else
alpar@1:          if (ssx->tm_lim >= 0.0 &&
alpar@1:              ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
alpar@1: #endif
alpar@1:          {  ret = 3;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* choose non-basic variable xN[q] */
alpar@1:          ssx_chuzc(ssx);
alpar@1:          /* if xN[q] cannot be chosen, the current basic solution is
alpar@1:             dual feasible and therefore optimal */
alpar@1:          if (ssx->q == 0)
alpar@1:          {  ret = 0;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* compute q-th column of the simplex table */
alpar@1:          ssx_eval_col(ssx);
alpar@1:          /* choose basic variable xB[p] */
alpar@1:          ssx_chuzr(ssx);
alpar@1:          /* if xB[p] cannot be chosen, the problem has no dual feasible
alpar@1:             solution (i.e. unbounded) */
alpar@1:          if (ssx->p == 0)
alpar@1:          {  ret = 1;
alpar@1:             break;
alpar@1:          }
alpar@1:          /* update values of basic variables */
alpar@1:          ssx_update_bbar(ssx);
alpar@1:          if (ssx->p > 0)
alpar@1:          {  /* compute p-th row of the inverse inv(B) */
alpar@1:             ssx_eval_rho(ssx);
alpar@1:             /* compute p-th row of the simplex table */
alpar@1:             ssx_eval_row(ssx);
alpar@1:             xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
alpar@1: #if 0
alpar@1:             /* update simplex multipliers */
alpar@1:             ssx_update_pi(ssx);
alpar@1: #endif
alpar@1:             /* update reduced costs of non-basic variables */
alpar@1:             ssx_update_cbar(ssx);
alpar@1:          }
alpar@1:          /* jump to the adjacent vertex of the polyhedron */
alpar@1:          ssx_change_basis(ssx);
alpar@1:          /* one simplex iteration has been performed */
alpar@1:          if (ssx->it_lim > 0) ssx->it_lim--;
alpar@1:          ssx->it_cnt++;
alpar@1:       }
alpar@1:       /* display final progress of the search */
alpar@1:       show_progress(ssx, 2);
alpar@1:       /* return to the calling program */
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /*----------------------------------------------------------------------
alpar@1: // ssx_driver - base driver to exact simplex method.
alpar@1: //
alpar@1: // This routine is a base driver to a version of the primal simplex
alpar@1: // method using exact (bignum) arithmetic.
alpar@1: //
alpar@1: // On exit the routine returns one of the following codes:
alpar@1: //
alpar@1: // 0 - optimal solution found;
alpar@1: // 1 - problem has no feasible solution;
alpar@1: // 2 - problem has unbounded solution;
alpar@1: // 3 - iterations limit exceeded (phase I);
alpar@1: // 4 - iterations limit exceeded (phase II);
alpar@1: // 5 - time limit exceeded (phase I);
alpar@1: // 6 - time limit exceeded (phase II);
alpar@1: // 7 - initial basis matrix is exactly singular.
alpar@1: ----------------------------------------------------------------------*/
alpar@1: 
alpar@1: int ssx_driver(SSX *ssx)
alpar@1: {     int m = ssx->m;
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 i, k, ret;
alpar@1:       ssx->tm_beg = xtime();
alpar@1:       /* factorize the initial basis matrix */
alpar@1:       if (ssx_factorize(ssx))
alpar@1:       {  xprintf("Initial basis matrix is singular\n");
alpar@1:          ret = 7;
alpar@1:          goto done;
alpar@1:       }
alpar@1:       /* compute values of basic variables */
alpar@1:       ssx_eval_bbar(ssx);
alpar@1:       /* check if the initial basic solution is primal feasible */
alpar@1:       for (i = 1; i <= m; i++)
alpar@1:       {  int t;
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:          {  /* x[k] has lower bound */
alpar@1:             if (mpq_cmp(bbar[i], lb[k]) < 0)
alpar@1:             {  /* which is violated */
alpar@1:                break;
alpar@1:             }
alpar@1:          }
alpar@1:          if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
alpar@1:          {  /* x[k] has upper bound */
alpar@1:             if (mpq_cmp(bbar[i], ub[k]) > 0)
alpar@1:             {  /* which is violated */
alpar@1:                break;
alpar@1:             }
alpar@1:          }
alpar@1:       }
alpar@1:       if (i > m)
alpar@1:       {  /* no basic variable violates its bounds */
alpar@1:          ret = 0;
alpar@1:          goto skip;
alpar@1:       }
alpar@1:       /* phase I: find primal feasible solution */
alpar@1:       ret = ssx_phase_I(ssx);
alpar@1:       switch (ret)
alpar@1:       {  case 0:
alpar@1:             ret = 0;
alpar@1:             break;
alpar@1:          case 1:
alpar@1:             xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
alpar@1:             ret = 1;
alpar@1:             break;
alpar@1:          case 2:
alpar@1:             xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@1:             ret = 3;
alpar@1:             break;
alpar@1:          case 3:
alpar@1:             xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@1:             ret = 5;
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(ret != ret);
alpar@1:       }
alpar@1:       /* compute values of basic variables (actually only the objective
alpar@1:          value needs to be computed) */
alpar@1:       ssx_eval_bbar(ssx);
alpar@1: skip: /* compute simplex multipliers */
alpar@1:       ssx_eval_pi(ssx);
alpar@1:       /* compute reduced costs of non-basic variables */
alpar@1:       ssx_eval_cbar(ssx);
alpar@1:       /* if phase I failed, do not start phase II */
alpar@1:       if (ret != 0) goto done;
alpar@1:       /* phase II: find optimal solution */
alpar@1:       ret = ssx_phase_II(ssx);
alpar@1:       switch (ret)
alpar@1:       {  case 0:
alpar@1:             xprintf("OPTIMAL SOLUTION FOUND\n");
alpar@1:             ret = 0;
alpar@1:             break;
alpar@1:          case 1:
alpar@1:             xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
alpar@1:             ret = 2;
alpar@1:             break;
alpar@1:          case 2:
alpar@1:             xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@1:             ret = 4;
alpar@1:             break;
alpar@1:          case 3:
alpar@1:             xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
alpar@1:             ret = 6;
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(ret != ret);
alpar@1:       }
alpar@1: done: /* decrease the time limit by the spent amount of time */
alpar@1:       if (ssx->tm_lim >= 0.0)
alpar@1: #if 0
alpar@1:       {  ssx->tm_lim -= utime() - ssx->tm_beg;
alpar@1: #else
alpar@1:       {  ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg);
alpar@1: #endif
alpar@1:          if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0;
alpar@1:       }
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /* eof */