/* glpssx02.c (simplex method, rational arithmetic) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
*  2009, 2010, 2011, 2013, 2017 Andrew Makhorin, Department for Applied
*  Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
*  reserved. E-mail: <mao@gnu.org>.
*
*  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 "env.h"
#include "glpssx.h"

static void show_progress(SSX *ssx, int phase)
{     /* this auxiliary routine displays information about progress of
         the search */
      int i, def = 0;
      for (i = 1; i <= ssx->m; i++)
         if (ssx->type[ssx->Q_col[i]] == SSX_FX) def++;
      xprintf("%s%6d:   %s = %22.15g   (%d)\n", phase == 1 ? " " : "*",
         ssx->it_cnt, phase == 1 ? "infsum" : "objval",
         mpq_get_d(ssx->bbar[0]), def);
#if 0
      ssx->tm_lag = utime();
#else
      ssx->tm_lag = xtime();
#endif
      return;
}

/*----------------------------------------------------------------------
// ssx_phase_I - find primal feasible solution.
//
// This routine implements phase I of the primal simplex method.
//
// On exit the routine returns one of the following codes:
//
// 0 - feasible solution found;
// 1 - problem has no feasible solution;
// 2 - iterations limit exceeded;
// 3 - time limit exceeded.
----------------------------------------------------------------------*/

int ssx_phase_I(SSX *ssx)
{     int m = ssx->m;
      int n = ssx->n;
      int *type = ssx->type;
      mpq_t *lb = ssx->lb;
      mpq_t *ub = ssx->ub;
      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;
      mpq_t *pi = ssx->pi;
      mpq_t *cbar = ssx->cbar;
      int *orig_type, orig_dir;
      mpq_t *orig_lb, *orig_ub, *orig_coef;
      int i, k, ret;
      /* save components of the original LP problem, which are changed
         by the routine */
      orig_type = xcalloc(1+m+n, sizeof(int));
      orig_lb = xcalloc(1+m+n, sizeof(mpq_t));
      orig_ub = xcalloc(1+m+n, sizeof(mpq_t));
      orig_coef = xcalloc(1+m+n, sizeof(mpq_t));
      for (k = 1; k <= m+n; k++)
      {  orig_type[k] = type[k];
         mpq_init(orig_lb[k]);
         mpq_set(orig_lb[k], lb[k]);
         mpq_init(orig_ub[k]);
         mpq_set(orig_ub[k], ub[k]);
      }
      orig_dir = ssx->dir;
      for (k = 0; k <= m+n; k++)
      {  mpq_init(orig_coef[k]);
         mpq_set(orig_coef[k], coef[k]);
      }
      /* build an artificial basic solution, which is primal feasible,
         and also build an auxiliary objective function to minimize the
         sum of infeasibilities for the original problem */
      ssx->dir = SSX_MIN;
      for (k = 0; k <= m+n; k++) mpq_set_si(coef[k], 0, 1);
      mpq_set_si(bbar[0], 0, 1);
      for (i = 1; i <= m; i++)
      {  int t;
         k = Q_col[i]; /* x[k] = xB[i] */
         t = type[k];
         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
         {  /* in the original problem x[k] has lower bound */
            if (mpq_cmp(bbar[i], lb[k]) < 0)
            {  /* which is violated */
               type[k] = SSX_UP;
               mpq_set(ub[k], lb[k]);
               mpq_set_si(lb[k], 0, 1);
               mpq_set_si(coef[k], -1, 1);
               mpq_add(bbar[0], bbar[0], ub[k]);
               mpq_sub(bbar[0], bbar[0], bbar[i]);
            }
         }
         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
         {  /* in the original problem x[k] has upper bound */
            if (mpq_cmp(bbar[i], ub[k]) > 0)
            {  /* which is violated */
               type[k] = SSX_LO;
               mpq_set(lb[k], ub[k]);
               mpq_set_si(ub[k], 0, 1);
               mpq_set_si(coef[k], +1, 1);
               mpq_add(bbar[0], bbar[0], bbar[i]);
               mpq_sub(bbar[0], bbar[0], lb[k]);
            }
         }
      }
      /* now the initial basic solution should be primal feasible due
         to changes of bounds of some basic variables, which turned to
         implicit artifical variables */
      /* compute simplex multipliers and reduced costs */
      ssx_eval_pi(ssx);
      ssx_eval_cbar(ssx);
      /* display initial progress of the search */
#if 1 /* 25/XI-2017 */
      if (ssx->msg_lev >= GLP_MSG_ON)
#endif
      show_progress(ssx, 1);
      /* main loop starts here */
      for (;;)
      {  /* display current progress of the search */
#if 1 /* 25/XI-2017 */
         if (ssx->msg_lev >= GLP_MSG_ON)
#endif
#if 0
         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
#else
         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
#endif
            show_progress(ssx, 1);
         /* we do not need to wait until all artificial variables have
            left the basis */
         if (mpq_sgn(bbar[0]) == 0)
         {  /* the sum of infeasibilities is zero, therefore the current
               solution is primal feasible for the original problem */
            ret = 0;
            break;
         }
         /* check if the iterations limit has been exhausted */
         if (ssx->it_lim == 0)
         {  ret = 2;
            break;
         }
         /* check if the time limit has been exhausted */
#if 0
         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
#else
         if (ssx->tm_lim >= 0.0 &&
             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
#endif
         {  ret = 3;
            break;
         }
         /* choose non-basic variable xN[q] */
         ssx_chuzc(ssx);
         /* if xN[q] cannot be chosen, the sum of infeasibilities is
            minimal but non-zero; therefore the original problem has no
            primal feasible solution */
         if (ssx->q == 0)
         {  ret = 1;
            break;
         }
         /* compute q-th column of the simplex table */
         ssx_eval_col(ssx);
         /* choose basic variable xB[p] */
         ssx_chuzr(ssx);
         /* the sum of infeasibilities cannot be negative, therefore
            the auxiliary lp problem cannot have unbounded solution */
         xassert(ssx->p != 0);
         /* update values of basic variables */
         ssx_update_bbar(ssx);
         if (ssx->p > 0)
         {  /* compute p-th row of the inverse inv(B) */
            ssx_eval_rho(ssx);
            /* compute p-th row of the simplex table */
            ssx_eval_row(ssx);
            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
            /* update simplex multipliers */
            ssx_update_pi(ssx);
            /* update reduced costs of non-basic variables */
            ssx_update_cbar(ssx);
         }
         /* xB[p] is leaving the basis; if it is implicit artificial
            variable, the corresponding residual vanishes; therefore
            bounds of this variable should be restored to the original
            values */
         if (ssx->p > 0)
         {  k = Q_col[ssx->p]; /* x[k] = xB[p] */
            if (type[k] != orig_type[k])
            {  /* x[k] is implicit artificial variable */
               type[k] = orig_type[k];
               mpq_set(lb[k], orig_lb[k]);
               mpq_set(ub[k], orig_ub[k]);
               xassert(ssx->p_stat == SSX_NL || ssx->p_stat == SSX_NU);
               ssx->p_stat = (ssx->p_stat == SSX_NL ? SSX_NU : SSX_NL);
               if (type[k] == SSX_FX) ssx->p_stat = SSX_NS;
               /* nullify the objective coefficient at x[k] */
               mpq_set_si(coef[k], 0, 1);
               /* since coef[k] has been changed, we need to compute
                  new reduced cost of x[k], which it will have in the
                  adjacent basis */
               /* the formula d[j] = cN[j] - pi' * N[j] is used (note
                  that the vector pi is not changed, because it depends
                  on objective coefficients at basic variables, but in
                  the adjacent basis, for which the vector pi has been
                  just recomputed, x[k] is non-basic) */
               if (k <= m)
               {  /* x[k] is auxiliary variable */
                  mpq_neg(cbar[ssx->q], pi[k]);
               }
               else
               {  /* x[k] is structural variable */
                  int ptr;
                  mpq_t temp;
                  mpq_init(temp);
                  mpq_set_si(cbar[ssx->q], 0, 1);
                  for (ptr = A_ptr[k-m]; ptr < A_ptr[k-m+1]; ptr++)
                  {  mpq_mul(temp, pi[A_ind[ptr]], A_val[ptr]);
                     mpq_add(cbar[ssx->q], cbar[ssx->q], temp);
                  }
                  mpq_clear(temp);
               }
            }
         }
         /* jump to the adjacent vertex of the polyhedron */
         ssx_change_basis(ssx);
         /* one simplex iteration has been performed */
         if (ssx->it_lim > 0) ssx->it_lim--;
         ssx->it_cnt++;
      }
      /* display final progress of the search */
#if 1 /* 25/XI-2017 */
      if (ssx->msg_lev >= GLP_MSG_ON)
#endif
      show_progress(ssx, 1);
      /* restore components of the original problem, which were changed
         by the routine */
      for (k = 1; k <= m+n; k++)
      {  type[k] = orig_type[k];
         mpq_set(lb[k], orig_lb[k]);
         mpq_clear(orig_lb[k]);
         mpq_set(ub[k], orig_ub[k]);
         mpq_clear(orig_ub[k]);
      }
      ssx->dir = orig_dir;
      for (k = 0; k <= m+n; k++)
      {  mpq_set(coef[k], orig_coef[k]);
         mpq_clear(orig_coef[k]);
      }
      xfree(orig_type);
      xfree(orig_lb);
      xfree(orig_ub);
      xfree(orig_coef);
      /* return to the calling program */
      return ret;
}

/*----------------------------------------------------------------------
// ssx_phase_II - find optimal solution.
//
// This routine implements phase II of the primal simplex method.
//
// On exit the routine returns one of the following codes:
//
// 0 - optimal solution found;
// 1 - problem has unbounded solution;
// 2 - iterations limit exceeded;
// 3 - time limit exceeded.
----------------------------------------------------------------------*/

int ssx_phase_II(SSX *ssx)
{     int ret;
      /* display initial progress of the search */
#if 1 /* 25/XI-2017 */
      if (ssx->msg_lev >= GLP_MSG_ON)
#endif
      show_progress(ssx, 2);
      /* main loop starts here */
      for (;;)
      {  /* display current progress of the search */
#if 1 /* 25/XI-2017 */
         if (ssx->msg_lev >= GLP_MSG_ON)
#endif
#if 0
         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
#else
         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
#endif
            show_progress(ssx, 2);
         /* check if the iterations limit has been exhausted */
         if (ssx->it_lim == 0)
         {  ret = 2;
            break;
         }
         /* check if the time limit has been exhausted */
#if 0
         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
#else
         if (ssx->tm_lim >= 0.0 &&
             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
#endif
         {  ret = 3;
            break;
         }
         /* choose non-basic variable xN[q] */
         ssx_chuzc(ssx);
         /* if xN[q] cannot be chosen, the current basic solution is
            dual feasible and therefore optimal */
         if (ssx->q == 0)
         {  ret = 0;
            break;
         }
         /* compute q-th column of the simplex table */
         ssx_eval_col(ssx);
         /* choose basic variable xB[p] */
         ssx_chuzr(ssx);
         /* if xB[p] cannot be chosen, the problem has no dual feasible
            solution (i.e. unbounded) */
         if (ssx->p == 0)
         {  ret = 1;
            break;
         }
         /* update values of basic variables */
         ssx_update_bbar(ssx);
         if (ssx->p > 0)
         {  /* compute p-th row of the inverse inv(B) */
            ssx_eval_rho(ssx);
            /* compute p-th row of the simplex table */
            ssx_eval_row(ssx);
            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
#if 0
            /* update simplex multipliers */
            ssx_update_pi(ssx);
#endif
            /* update reduced costs of non-basic variables */
            ssx_update_cbar(ssx);
         }
         /* jump to the adjacent vertex of the polyhedron */
         ssx_change_basis(ssx);
         /* one simplex iteration has been performed */
         if (ssx->it_lim > 0) ssx->it_lim--;
         ssx->it_cnt++;
      }
      /* display final progress of the search */
#if 1 /* 25/XI-2017 */
      if (ssx->msg_lev >= GLP_MSG_ON)
#endif
      show_progress(ssx, 2);
      /* return to the calling program */
      return ret;
}

/*----------------------------------------------------------------------
// ssx_driver - base driver to exact simplex method.
//
// This routine is a base driver to a version of the primal simplex
// method using exact (bignum) arithmetic.
//
// On exit the routine returns one of the following codes:
//
// 0 - optimal solution found;
// 1 - problem has no feasible solution;
// 2 - problem has unbounded solution;
// 3 - iterations limit exceeded (phase I);
// 4 - iterations limit exceeded (phase II);
// 5 - time limit exceeded (phase I);
// 6 - time limit exceeded (phase II);
// 7 - initial basis matrix is exactly singular.
----------------------------------------------------------------------*/

int ssx_driver(SSX *ssx)
{     int m = ssx->m;
      int *type = ssx->type;
      mpq_t *lb = ssx->lb;
      mpq_t *ub = ssx->ub;
      int *Q_col = ssx->Q_col;
      mpq_t *bbar = ssx->bbar;
      int i, k, ret;
      ssx->tm_beg = xtime();
      /* factorize the initial basis matrix */
      if (ssx_factorize(ssx))
#if 0 /* 25/XI-2017 */
      {  xprintf("Initial basis matrix is singular\n");
#else
      {  if (ssx->msg_lev >= GLP_MSG_ERR)
            xprintf("Initial basis matrix is singular\n");
#endif
         ret = 7;
         goto done;
      }
      /* compute values of basic variables */
      ssx_eval_bbar(ssx);
      /* check if the initial basic solution is primal feasible */
      for (i = 1; i <= m; i++)
      {  int t;
         k = Q_col[i]; /* x[k] = xB[i] */
         t = type[k];
         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
         {  /* x[k] has lower bound */
            if (mpq_cmp(bbar[i], lb[k]) < 0)
            {  /* which is violated */
               break;
            }
         }
         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
         {  /* x[k] has upper bound */
            if (mpq_cmp(bbar[i], ub[k]) > 0)
            {  /* which is violated */
               break;
            }
         }
      }
      if (i > m)
      {  /* no basic variable violates its bounds */
         ret = 0;
         goto skip;
      }
      /* phase I: find primal feasible solution */
      ret = ssx_phase_I(ssx);
      switch (ret)
      {  case 0:
            ret = 0;
            break;
         case 1:
#if 1 /* 25/XI-2017 */
            if (ssx->msg_lev >= GLP_MSG_ALL)
#endif
            xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
            ret = 1;
            break;
         case 2:
#if 1 /* 25/XI-2017 */
            if (ssx->msg_lev >= GLP_MSG_ALL)
#endif
            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 3;
            break;
         case 3:
#if 1 /* 25/XI-2017 */
            if (ssx->msg_lev >= GLP_MSG_ALL)
#endif
            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 5;
            break;
         default:
            xassert(ret != ret);
      }
      /* compute values of basic variables (actually only the objective
         value needs to be computed) */
      ssx_eval_bbar(ssx);
skip: /* compute simplex multipliers */
      ssx_eval_pi(ssx);
      /* compute reduced costs of non-basic variables */
      ssx_eval_cbar(ssx);
      /* if phase I failed, do not start phase II */
      if (ret != 0) goto done;
      /* phase II: find optimal solution */
      ret = ssx_phase_II(ssx);
      switch (ret)
      {  case 0:
#if 1 /* 25/XI-2017 */
            if (ssx->msg_lev >= GLP_MSG_ALL)
#endif
            xprintf("OPTIMAL SOLUTION FOUND\n");
            ret = 0;
            break;
         case 1:
#if 1 /* 25/XI-2017 */
            if (ssx->msg_lev >= GLP_MSG_ALL)
#endif
            xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
            ret = 2;
            break;
         case 2:
#if 1 /* 25/XI-2017 */
            if (ssx->msg_lev >= GLP_MSG_ALL)
#endif
            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 4;
            break;
         case 3:
#if 1 /* 25/XI-2017 */
            if (ssx->msg_lev >= GLP_MSG_ALL)
#endif
            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 6;
            break;
         default:
            xassert(ret != ret);
      }
done: /* decrease the time limit by the spent amount of time */
      if (ssx->tm_lim >= 0.0)
#if 0
      {  ssx->tm_lim -= utime() - ssx->tm_beg;
#else
      {  ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg);
#endif
         if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0;
      }
      return ret;
}

/* eof */