praat/external/glpk/glpapi12.c

/* glpapi12.c (basis factorization and simplex tableau routines) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
*  2009, 2010 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 "glpapi.h"

/***********************************************************************
*  NAME
*
*  glp_bf_exists - check if the basis factorization exists
*
*  SYNOPSIS
*
*  int glp_bf_exists(glp_prob *lp);
*
*  RETURNS
*
*  If the basis factorization for the current basis associated with
*  the specified problem object exists and therefore is available for
*  computations, the routine glp_bf_exists returns non-zero. Otherwise
*  the routine returns zero. */

int glp_bf_exists(glp_prob *lp)
{     int ret;
      ret = (lp->m == 0 || lp->valid);
      return ret;
}

/***********************************************************************
*  NAME
*
*  glp_factorize - compute the basis factorization
*
*  SYNOPSIS
*
*  int glp_factorize(glp_prob *lp);
*
*  DESCRIPTION
*
*  The routine glp_factorize computes the basis factorization for the
*  current basis associated with the specified problem object.
*
*  RETURNS
*
*  0  The basis factorization has been successfully computed.
*
*  GLP_EBADB
*     The basis matrix is invalid, i.e. the number of basic (auxiliary
*     and structural) variables differs from the number of rows in the
*     problem object.
*
*  GLP_ESING
*     The basis matrix is singular within the working precision.
*
*  GLP_ECOND
*     The basis matrix is ill-conditioned. */

static int b_col(void *info, int j, int ind[], double val[])
{     glp_prob *lp = info;
      int m = lp->m;
      GLPAIJ *aij;
      int k, len;
      xassert(1 <= j && j <= m);
      /* determine the ordinal number of basic auxiliary or structural
         variable x[k] corresponding to basic variable xB[j] */
      k = lp->head[j];
      /* build j-th column of the basic matrix, which is k-th column of
         the scaled augmented matrix (I | -R*A*S) */
      if (k <= m)
      {  /* x[k] is auxiliary variable */
         len = 1;
         ind[1] = k;
         val[1] = 1.0;
      }
      else
      {  /* x[k] is structural variable */
         len = 0;
         for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next)
         {  len++;
            ind[len] = aij->row->i;
            val[len] = - aij->row->rii * aij->val * aij->col->sjj;
         }
      }
      return len;
}

static void copy_bfcp(glp_prob *lp);

int glp_factorize(glp_prob *lp)
{     int m = lp->m;
      int n = lp->n;
      GLPROW **row = lp->row;
      GLPCOL **col = lp->col;
      int *head = lp->head;
      int j, k, stat, ret;
      /* invalidate the basis factorization */
      lp->valid = 0;
      /* build the basis header */
      j = 0;
      for (k = 1; k <= m+n; k++)
      {  if (k <= m)
         {  stat = row[k]->stat;
            row[k]->bind = 0;
         }
         else
         {  stat = col[k-m]->stat;
            col[k-m]->bind = 0;
         }
         if (stat == GLP_BS)
         {  j++;
            if (j > m)
            {  /* too many basic variables */
               ret = GLP_EBADB;
               goto fini;
            }
            head[j] = k;
            if (k <= m)
               row[k]->bind = j;
            else
               col[k-m]->bind = j;
         }
      }
      if (j < m)
      {  /* too few basic variables */
         ret = GLP_EBADB;
         goto fini;
      }
      /* try to factorize the basis matrix */
      if (m > 0)
      {  if (lp->bfd == NULL)
         {  lp->bfd = bfd_create_it();
            copy_bfcp(lp);
         }
         switch (bfd_factorize(lp->bfd, m, lp->head, b_col, lp))
         {  case 0:
               /* ok */
               break;
            case BFD_ESING:
               /* singular matrix */
               ret = GLP_ESING;
               goto fini;
            case BFD_ECOND:
               /* ill-conditioned matrix */
               ret = GLP_ECOND;
               goto fini;
            default:
               xassert(lp != lp);
         }
         lp->valid = 1;
      }
      /* factorization successful */
      ret = 0;
fini: /* bring the return code to the calling program */
      return ret;
}

/***********************************************************************
*  NAME
*
*  glp_bf_updated - check if the basis factorization has been updated
*
*  SYNOPSIS
*
*  int glp_bf_updated(glp_prob *lp);
*
*  RETURNS
*
*  If the basis factorization has been just computed from scratch, the
*  routine glp_bf_updated returns zero. Otherwise, if the factorization
*  has been updated one or more times, the routine returns non-zero. */

int glp_bf_updated(glp_prob *lp)
{     int cnt;
      if (!(lp->m == 0 || lp->valid))
         xerror("glp_bf_update: basis factorization does not exist\n");
#if 0 /* 15/XI-2009 */
      cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt);
#else
      cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd));
#endif
      return cnt;
}

/***********************************************************************
*  NAME
*
*  glp_get_bfcp - retrieve basis factorization control parameters
*
*  SYNOPSIS
*
*  void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
*
*  DESCRIPTION
*
*  The routine glp_get_bfcp retrieves control parameters, which are
*  used on computing and updating the basis factorization associated
*  with the specified problem object.
*
*  Current values of control parameters are stored by the routine in
*  a glp_bfcp structure, which the parameter parm points to. */

void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm)
{     glp_bfcp *bfcp = lp->bfcp;
      if (bfcp == NULL)
      {  parm->type = GLP_BF_FT;
         parm->lu_size = 0;
         parm->piv_tol = 0.10;
         parm->piv_lim = 4;
         parm->suhl = GLP_ON;
         parm->eps_tol = 1e-15;
         parm->max_gro = 1e+10;
         parm->nfs_max = 100;
         parm->upd_tol = 1e-6;
         parm->nrs_max = 100;
         parm->rs_size = 0;
      }
      else
         memcpy(parm, bfcp, sizeof(glp_bfcp));
      return;
}

/***********************************************************************
*  NAME
*
*  glp_set_bfcp - change basis factorization control parameters
*
*  SYNOPSIS
*
*  void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
*
*  DESCRIPTION
*
*  The routine glp_set_bfcp changes control parameters, which are used
*  by internal GLPK routines in computing and updating the basis
*  factorization associated with the specified problem object.
*
*  New values of the control parameters should be passed in a structure
*  glp_bfcp, which the parameter parm points to.
*
*  The parameter parm can be specified as NULL, in which case all
*  control parameters are reset to their default values. */

#if 0 /* 15/XI-2009 */
static void copy_bfcp(glp_prob *lp)
{     glp_bfcp _parm, *parm = &_parm;
      BFD *bfd = lp->bfd;
      glp_get_bfcp(lp, parm);
      xassert(bfd != NULL);
      bfd->type = parm->type;
      bfd->lu_size = parm->lu_size;
      bfd->piv_tol = parm->piv_tol;
      bfd->piv_lim = parm->piv_lim;
      bfd->suhl = parm->suhl;
      bfd->eps_tol = parm->eps_tol;
      bfd->max_gro = parm->max_gro;
      bfd->nfs_max = parm->nfs_max;
      bfd->upd_tol = parm->upd_tol;
      bfd->nrs_max = parm->nrs_max;
      bfd->rs_size = parm->rs_size;
      return;
}
#else
static void copy_bfcp(glp_prob *lp)
{     glp_bfcp _parm, *parm = &_parm;
      glp_get_bfcp(lp, parm);
      bfd_set_parm(lp->bfd, parm);
      return;
}
#endif

void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm)
{     glp_bfcp *bfcp = lp->bfcp;
      if (parm == NULL)
      {  /* reset to default values */
         if (bfcp != NULL)
            xfree(bfcp), lp->bfcp = NULL;
      }
      else
      {  /* set to specified values */
         if (bfcp == NULL)
            bfcp = lp->bfcp = xmalloc(sizeof(glp_bfcp));
         memcpy(bfcp, parm, sizeof(glp_bfcp));
         if (!(bfcp->type == GLP_BF_FT || bfcp->type == GLP_BF_BG ||
               bfcp->type == GLP_BF_GR))
            xerror("glp_set_bfcp: type = %d; invalid parameter\n",
               bfcp->type);
         if (bfcp->lu_size < 0)
            xerror("glp_set_bfcp: lu_size = %d; invalid parameter\n",
               bfcp->lu_size);
         if (!(0.0 < bfcp->piv_tol && bfcp->piv_tol < 1.0))
            xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n",
               bfcp->piv_tol);
         if (bfcp->piv_lim < 1)
            xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n",
               bfcp->piv_lim);
         if (!(bfcp->suhl == GLP_ON || bfcp->suhl == GLP_OFF))
            xerror("glp_set_bfcp: suhl = %d; invalid parameter\n",
               bfcp->suhl);
         if (!(0.0 <= bfcp->eps_tol && bfcp->eps_tol <= 1e-6))
            xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n",
               bfcp->eps_tol);
         if (bfcp->max_gro < 1.0)
            xerror("glp_set_bfcp: max_gro = %g; invalid parameter\n",
               bfcp->max_gro);
         if (!(1 <= bfcp->nfs_max && bfcp->nfs_max <= 32767))
            xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n",
               bfcp->nfs_max);
         if (!(0.0 < bfcp->upd_tol && bfcp->upd_tol < 1.0))
            xerror("glp_set_bfcp: upd_tol = %g; invalid parameter\n",
               bfcp->upd_tol);
         if (!(1 <= bfcp->nrs_max && bfcp->nrs_max <= 32767))
            xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n",
               bfcp->nrs_max);
         if (bfcp->rs_size < 0)
            xerror("glp_set_bfcp: rs_size = %d; invalid parameter\n",
               bfcp->nrs_max);
         if (bfcp->rs_size == 0)
            bfcp->rs_size = 20 * bfcp->nrs_max;
      }
      if (lp->bfd != NULL) copy_bfcp(lp);
      return;
}

/***********************************************************************
*  NAME
*
*  glp_get_bhead - retrieve the basis header information
*
*  SYNOPSIS
*
*  int glp_get_bhead(glp_prob *lp, int k);
*
*  DESCRIPTION
*
*  The routine glp_get_bhead returns the basis header information for
*  the current basis associated with the specified problem object.
*
*  RETURNS
*
*  If xB[k], 1 <= k <= m, is i-th auxiliary variable (1 <= i <= m), the
*  routine returns i. Otherwise, if xB[k] is j-th structural variable
*  (1 <= j <= n), the routine returns m+j. Here m is the number of rows
*  and n is the number of columns in the problem object. */

int glp_get_bhead(glp_prob *lp, int k)
{     if (!(lp->m == 0 || lp->valid))
         xerror("glp_get_bhead: basis factorization does not exist\n");
      if (!(1 <= k && k <= lp->m))
         xerror("glp_get_bhead: k = %d; index out of range\n", k);
      return lp->head[k];
}

/***********************************************************************
*  NAME
*
*  glp_get_row_bind - retrieve row index in the basis header
*
*  SYNOPSIS
*
*  int glp_get_row_bind(glp_prob *lp, int i);
*
*  RETURNS
*
*  The routine glp_get_row_bind returns the index k of basic variable
*  xB[k], 1 <= k <= m, which is i-th auxiliary variable, 1 <= i <= m,
*  in the current basis associated with the specified problem object,
*  where m is the number of rows. However, if i-th auxiliary variable
*  is non-basic, the routine returns zero. */

int glp_get_row_bind(glp_prob *lp, int i)
{     if (!(lp->m == 0 || lp->valid))
         xerror("glp_get_row_bind: basis factorization does not exist\n"
            );
      if (!(1 <= i && i <= lp->m))
         xerror("glp_get_row_bind: i = %d; row number out of range\n",
            i);
      return lp->row[i]->bind;
}

/***********************************************************************
*  NAME
*
*  glp_get_col_bind - retrieve column index in the basis header
*
*  SYNOPSIS
*
*  int glp_get_col_bind(glp_prob *lp, int j);
*
*  RETURNS
*
*  The routine glp_get_col_bind returns the index k of basic variable
*  xB[k], 1 <= k <= m, which is j-th structural variable, 1 <= j <= n,
*  in the current basis associated with the specified problem object,
*  where m is the number of rows, n is the number of columns. However,
*  if j-th structural variable is non-basic, the routine returns zero.*/

int glp_get_col_bind(glp_prob *lp, int j)
{     if (!(lp->m == 0 || lp->valid))
         xerror("glp_get_col_bind: basis factorization does not exist\n"
            );
      if (!(1 <= j && j <= lp->n))
         xerror("glp_get_col_bind: j = %d; column number out of range\n"
            , j);
      return lp->col[j]->bind;
}

/***********************************************************************
*  NAME
*
*  glp_ftran - perform forward transformation (solve system B*x = b)
*
*  SYNOPSIS
*
*  void glp_ftran(glp_prob *lp, double x[]);
*
*  DESCRIPTION
*
*  The routine glp_ftran performs forward transformation, i.e. solves
*  the system B*x = b, where B is the basis matrix corresponding to the
*  current basis for the specified problem object, x is the vector of
*  unknowns to be computed, b is the vector of right-hand sides.
*
*  On entry elements of the vector b should be stored in dense format
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
*  the routine stores elements of the vector x in the same locations.
*
*  SCALING/UNSCALING
*
*  Let A~ = (I | -A) is the augmented constraint matrix of the original
*  (unscaled) problem. In the scaled LP problem instead the matrix A the
*  scaled matrix A" = R*A*S is actually used, so
*
*     A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) =
*                                                                    (1)
*         = R*(I | A)*S~ = R*A~*S~,
*
*  is the scaled augmented constraint matrix, where R and S are diagonal
*  scaling matrices used to scale rows and columns of the matrix A, and
*
*     S~ = diag(inv(R) | S)                                          (2)
*
*  is an augmented diagonal scaling matrix.
*
*  By definition:
*
*     A~ = (B | N),                                                  (3)
*
*  where B is the basic matrix, which consists of basic columns of the
*  augmented constraint matrix A~, and N is a matrix, which consists of
*  non-basic columns of A~. From (1) it follows that:
*
*     A~" = (B" | N") = (R*B*SB | R*N*SN),                           (4)
*
*  where SB and SN are parts of the augmented scaling matrix S~, which
*  correspond to basic and non-basic variables, respectively. Therefore
*
*     B" = R*B*SB,                                                   (5)
*
*  which is the scaled basis matrix. */

void glp_ftran(glp_prob *lp, double x[])
{     int m = lp->m;
      GLPROW **row = lp->row;
      GLPCOL **col = lp->col;
      int i, k;
      /* B*x = b ===> (R*B*SB)*(inv(SB)*x) = R*b ===>
         B"*x" = b", where b" = R*b, x = SB*x" */
      if (!(m == 0 || lp->valid))
         xerror("glp_ftran: basis factorization does not exist\n");
      /* b" := R*b */
      for (i = 1; i <= m; i++)
         x[i] *= row[i]->rii;
      /* x" := inv(B")*b" */
      if (m > 0) bfd_ftran(lp->bfd, x);
      /* x := SB*x" */
      for (i = 1; i <= m; i++)
      {  k = lp->head[i];
         if (k <= m)
            x[i] /= row[k]->rii;
         else
            x[i] *= col[k-m]->sjj;
      }
      return;
}

/***********************************************************************
*  NAME
*
*  glp_btran - perform backward transformation (solve system B'*x = b)
*
*  SYNOPSIS
*
*  void glp_btran(glp_prob *lp, double x[]);
*
*  DESCRIPTION
*
*  The routine glp_btran performs backward transformation, i.e. solves
*  the system B'*x = b, where B' is a matrix transposed to the basis
*  matrix corresponding to the current basis for the specified problem
*  problem object, x is the vector of unknowns to be computed, b is the
*  vector of right-hand sides.
*
*  On entry elements of the vector b should be stored in dense format
*  in locations x[1], ..., x[m], where m is the number of rows. On exit
*  the routine stores elements of the vector x in the same locations.
*
*  SCALING/UNSCALING
*
*  See comments to the routine glp_ftran. */

void glp_btran(glp_prob *lp, double x[])
{     int m = lp->m;
      GLPROW **row = lp->row;
      GLPCOL **col = lp->col;
      int i, k;
      /* B'*x = b ===> (SB*B'*R)*(inv(R)*x) = SB*b ===>
         (B")'*x" = b", where b" = SB*b, x = R*x" */
      if (!(m == 0 || lp->valid))
         xerror("glp_btran: basis factorization does not exist\n");
      /* b" := SB*b */
      for (i = 1; i <= m; i++)
      {  k = lp->head[i];
         if (k <= m)
            x[i] /= row[k]->rii;
         else
            x[i] *= col[k-m]->sjj;
      }
      /* x" := inv[(B")']*b" */
      if (m > 0) bfd_btran(lp->bfd, x);
      /* x := R*x" */
      for (i = 1; i <= m; i++)
         x[i] *= row[i]->rii;
      return;
}

/***********************************************************************
*  NAME
*
*  glp_warm_up - "warm up" LP basis
*
*  SYNOPSIS
*
*  int glp_warm_up(glp_prob *P);
*
*  DESCRIPTION
*
*  The routine glp_warm_up "warms up" the LP basis for the specified
*  problem object using current statuses assigned to rows and columns
*  (that is, to auxiliary and structural variables).
*
*  This operation includes computing factorization of the basis matrix
*  (if it does not exist), computing primal and dual components of basic
*  solution, and determining the solution status.
*
*  RETURNS
*
*  0  The operation has been successfully performed.
*
*  GLP_EBADB
*     The basis matrix is invalid, i.e. the number of basic (auxiliary
*     and structural) variables differs from the number of rows in the
*     problem object.
*
*  GLP_ESING
*     The basis matrix is singular within the working precision.
*
*  GLP_ECOND
*     The basis matrix is ill-conditioned. */

int glp_warm_up(glp_prob *P)
{     GLPROW *row;
      GLPCOL *col;
      GLPAIJ *aij;
      int i, j, type, ret;
      double eps, temp, *work;
      /* invalidate basic solution */
      P->pbs_stat = P->dbs_stat = GLP_UNDEF;
      P->obj_val = 0.0;
      P->some = 0;
      for (i = 1; i <= P->m; i++)
      {  row = P->row[i];
         row->prim = row->dual = 0.0;
      }
      for (j = 1; j <= P->n; j++)
      {  col = P->col[j];
         col->prim = col->dual = 0.0;
      }
      /* compute the basis factorization, if necessary */
      if (!glp_bf_exists(P))
      {  ret = glp_factorize(P);
         if (ret != 0) goto done;
      }
      /* allocate working array */
      work = xcalloc(1+P->m, sizeof(double));
      /* determine and store values of non-basic variables, compute
         vector (- N * xN) */
      for (i = 1; i <= P->m; i++)
         work[i] = 0.0;
      for (i = 1; i <= P->m; i++)
      {  row = P->row[i];
         if (row->stat == GLP_BS)
            continue;
         else if (row->stat == GLP_NL)
            row->prim = row->lb;
         else if (row->stat == GLP_NU)
            row->prim = row->ub;
         else if (row->stat == GLP_NF)
            row->prim = 0.0;
         else if (row->stat == GLP_NS)
            row->prim = row->lb;
         else
            xassert(row != row);
         /* N[j] is i-th column of matrix (I|-A) */
         work[i] -= row->prim;
      }
      for (j = 1; j <= P->n; j++)
      {  col = P->col[j];
         if (col->stat == GLP_BS)
            continue;
         else if (col->stat == GLP_NL)
            col->prim = col->lb;
         else if (col->stat == GLP_NU)
            col->prim = col->ub;
         else if (col->stat == GLP_NF)
            col->prim = 0.0;
         else if (col->stat == GLP_NS)
            col->prim = col->lb;
         else
            xassert(col != col);
         /* N[j] is (m+j)-th column of matrix (I|-A) */
         if (col->prim != 0.0)
         {  for (aij = col->ptr; aij != NULL; aij = aij->c_next)
               work[aij->row->i] += aij->val * col->prim;
         }
      }
      /* compute vector of basic variables xB = - inv(B) * N * xN */
      glp_ftran(P, work);
      /* store values of basic variables, check primal feasibility */
      P->pbs_stat = GLP_FEAS;
      for (i = 1; i <= P->m; i++)
      {  row = P->row[i];
         if (row->stat != GLP_BS)
            continue;
         row->prim = work[row->bind];
         type = row->type;
         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
         {  eps = 1e-6 + 1e-9 * fabs(row->lb);
            if (row->prim < row->lb - eps)
               P->pbs_stat = GLP_INFEAS;
         }
         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
         {  eps = 1e-6 + 1e-9 * fabs(row->ub);
            if (row->prim > row->ub + eps)
               P->pbs_stat = GLP_INFEAS;
         }
      }
      for (j = 1; j <= P->n; j++)
      {  col = P->col[j];
         if (col->stat != GLP_BS)
            continue;
         col->prim = work[col->bind];
         type = col->type;
         if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
         {  eps = 1e-6 + 1e-9 * fabs(col->lb);
            if (col->prim < col->lb - eps)
               P->pbs_stat = GLP_INFEAS;
         }
         if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
         {  eps = 1e-6 + 1e-9 * fabs(col->ub);
            if (col->prim > col->ub + eps)
               P->pbs_stat = GLP_INFEAS;
         }
      }
      /* compute value of the objective function */
      P->obj_val = P->c0;
      for (j = 1; j <= P->n; j++)
      {  col = P->col[j];
         P->obj_val += col->coef * col->prim;
      }
      /* build vector cB of objective coefficients at basic variables */
      for (i = 1; i <= P->m; i++)
         work[i] = 0.0;
      for (j = 1; j <= P->n; j++)
      {  col = P->col[j];
         if (col->stat == GLP_BS)
            work[col->bind] = col->coef;
      }
      /* compute vector of simplex multipliers pi = inv(B') * cB */
      glp_btran(P, work);
      /* compute and store reduced costs of non-basic variables d[j] =
         c[j] - N'[j] * pi, check dual feasibility */
      P->dbs_stat = GLP_FEAS;
      for (i = 1; i <= P->m; i++)
      {  row = P->row[i];
         if (row->stat == GLP_BS)
         {  row->dual = 0.0;
            continue;
         }
         /* N[j] is i-th column of matrix (I|-A) */
         row->dual = - work[i];
         type = row->type;
         temp = (P->dir == GLP_MIN ? + row->dual : - row->dual);
         if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
             (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
            P->dbs_stat = GLP_INFEAS;
      }
      for (j = 1; j <= P->n; j++)
      {  col = P->col[j];
         if (col->stat == GLP_BS)
         {  col->dual = 0.0;
            continue;
         }
         /* N[j] is (m+j)-th column of matrix (I|-A) */
         col->dual = col->coef;
         for (aij = col->ptr; aij != NULL; aij = aij->c_next)
            col->dual += aij->val * work[aij->row->i];
         type = col->type;
         temp = (P->dir == GLP_MIN ? + col->dual : - col->dual);
         if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
             (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
            P->dbs_stat = GLP_INFEAS;
      }
      /* free working array */
      xfree(work);
      ret = 0;
done: return ret;
}

/***********************************************************************
*  NAME
*
*  glp_eval_tab_row - compute row of the simplex tableau
*
*  SYNOPSIS
*
*  int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]);
*
*  DESCRIPTION
*
*  The routine glp_eval_tab_row computes a row of the current simplex
*  tableau for the basic variable, which is specified by the number k:
*  if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
*  x[k] is (k-m)-th structural variable, where m is number of rows, and
*  n is number of columns. The current basis must be available.
*
*  The routine stores column indices and numerical values of non-zero
*  elements of the computed row using sparse format to the locations
*  ind[1], ..., ind[len] and val[1], ..., val[len], respectively, where
*  0 <= len <= n is number of non-zeros returned on exit.
*
*  Element indices stored in the array ind have the same sense as the
*  index k, i.e. indices 1 to m denote auxiliary variables and indices
*  m+1 to m+n denote structural ones (all these variables are obviously
*  non-basic by definition).
*
*  The computed row shows how the specified basic variable x[k] = xB[i]
*  depends on non-basic variables:
*
*     xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n],
*
*  where alfa[i,j] are elements of the simplex table row, xN[j] are
*  non-basic (auxiliary and structural) variables.
*
*  RETURNS
*
*  The routine returns number of non-zero elements in the simplex table
*  row stored in the arrays ind and val.
*
*  BACKGROUND
*
*  The system of equality constraints of the LP problem is:
*
*     xR = A * xS,                                                   (1)
*
*  where xR is the vector of auxliary variables, xS is the vector of
*  structural variables, A is the matrix of constraint coefficients.
*
*  The system (1) can be written in homogenous form as follows:
*
*     A~ * x = 0,                                                    (2)
*
*  where A~ = (I | -A) is the augmented constraint matrix (has m rows
*  and m+n columns), x = (xR | xS) is the vector of all (auxiliary and
*  structural) variables.
*
*  By definition for the current basis we have:
*
*     A~ = (B | N),                                                  (3)
*
*  where B is the basis matrix. Thus, the system (2) can be written as:
*
*     B * xB + N * xN = 0.                                           (4)
*
*  From (4) it follows that:
*
*     xB = A^ * xN,                                                  (5)
*
*  where the matrix
*
*     A^ = - inv(B) * N                                              (6)
*
*  is called the simplex table.
*
*  It is understood that i-th row of the simplex table is:
*
*     e * A^ = - e * inv(B) * N,                                     (7)
*
*  where e is a unity vector with e[i] = 1.
*
*  To compute i-th row of the simplex table the routine first computes
*  i-th row of the inverse:
*
*     rho = inv(B') * e,                                             (8)
*
*  where B' is a matrix transposed to B, and then computes elements of
*  i-th row of the simplex table as scalar products:
*
*     alfa[i,j] = - rho * N[j]   for all j,                          (9)
*
*  where N[j] is a column of the augmented constraint matrix A~, which
*  corresponds to some non-basic auxiliary or structural variable. */

int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[])
{     int m = lp->m;
      int n = lp->n;
      int i, t, len, lll, *iii;
      double alfa, *rho, *vvv;
      if (!(m == 0 || lp->valid))
         xerror("glp_eval_tab_row: basis factorization does not exist\n"
            );
      if (!(1 <= k && k <= m+n))
         xerror("glp_eval_tab_row: k = %d; variable number out of range"
            , k);
      /* determine xB[i] which corresponds to x[k] */
      if (k <= m)
         i = glp_get_row_bind(lp, k);
      else
         i = glp_get_col_bind(lp, k-m);
      if (i == 0)
         xerror("glp_eval_tab_row: k = %d; variable must be basic", k);
      xassert(1 <= i && i <= m);
      /* allocate working arrays */
      rho = xcalloc(1+m, sizeof(double));
      iii = xcalloc(1+m, sizeof(int));
      vvv = xcalloc(1+m, sizeof(double));
      /* compute i-th row of the inverse; see (8) */
      for (t = 1; t <= m; t++) rho[t] = 0.0;
      rho[i] = 1.0;
      glp_btran(lp, rho);
      /* compute i-th row of the simplex table */
      len = 0;
      for (k = 1; k <= m+n; k++)
      {  if (k <= m)
         {  /* x[k] is auxiliary variable, so N[k] is a unity column */
            if (glp_get_row_stat(lp, k) == GLP_BS) continue;
            /* compute alfa[i,j]; see (9) */
            alfa = - rho[k];
         }
         else
         {  /* x[k] is structural variable, so N[k] is a column of the
               original constraint matrix A with negative sign */
            if (glp_get_col_stat(lp, k-m) == GLP_BS) continue;
            /* compute alfa[i,j]; see (9) */
            lll = glp_get_mat_col(lp, k-m, iii, vvv);
            alfa = 0.0;
            for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t];
         }
         /* store alfa[i,j] */
         if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa;
      }
      xassert(len <= n);
      /* free working arrays */
      xfree(rho);
      xfree(iii);
      xfree(vvv);
      /* return to the calling program */
      return len;
}

/***********************************************************************
*  NAME
*
*  glp_eval_tab_col - compute column of the simplex tableau
*
*  SYNOPSIS
*
*  int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]);
*
*  DESCRIPTION
*
*  The routine glp_eval_tab_col computes a column of the current simplex
*  table for the non-basic variable, which is specified by the number k:
*  if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
*  x[k] is (k-m)-th structural variable, where m is number of rows, and
*  n is number of columns. The current basis must be available.
*
*  The routine stores row indices and numerical values of non-zero
*  elements of the computed column using sparse format to the locations
*  ind[1], ..., ind[len] and val[1], ..., val[len] respectively, where
*  0 <= len <= m is number of non-zeros returned on exit.
*
*  Element indices stored in the array ind have the same sense as the
*  index k, i.e. indices 1 to m denote auxiliary variables and indices
*  m+1 to m+n denote structural ones (all these variables are obviously
*  basic by the definition).
*
*  The computed column shows how basic variables depend on the specified
*  non-basic variable x[k] = xN[j]:
*
*     xB[1] = ... + alfa[1,j]*xN[j] + ...
*     xB[2] = ... + alfa[2,j]*xN[j] + ...
*              . . . . . .
*     xB[m] = ... + alfa[m,j]*xN[j] + ...
*
*  where alfa[i,j] are elements of the simplex table column, xB[i] are
*  basic (auxiliary and structural) variables.
*
*  RETURNS
*
*  The routine returns number of non-zero elements in the simplex table
*  column stored in the arrays ind and val.
*
*  BACKGROUND
*
*  As it was explained in comments to the routine glp_eval_tab_row (see
*  above) the simplex table is the following matrix:
*
*     A^ = - inv(B) * N.                                             (1)
*
*  Therefore j-th column of the simplex table is:
*
*     A^ * e = - inv(B) * N * e = - inv(B) * N[j],                   (2)
*
*  where e is a unity vector with e[j] = 1, B is the basis matrix, N[j]
*  is a column of the augmented constraint matrix A~, which corresponds
*  to the given non-basic auxiliary or structural variable. */

int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[])
{     int m = lp->m;
      int n = lp->n;
      int t, len, stat;
      double *col;
      if (!(m == 0 || lp->valid))
         xerror("glp_eval_tab_col: basis factorization does not exist\n"
            );
      if (!(1 <= k && k <= m+n))
         xerror("glp_eval_tab_col: k = %d; variable number out of range"
            , k);
      if (k <= m)
         stat = glp_get_row_stat(lp, k);
      else
         stat = glp_get_col_stat(lp, k-m);
      if (stat == GLP_BS)
         xerror("glp_eval_tab_col: k = %d; variable must be non-basic",
            k);
      /* obtain column N[k] with negative sign */
      col = xcalloc(1+m, sizeof(double));
      for (t = 1; t <= m; t++) col[t] = 0.0;
      if (k <= m)
      {  /* x[k] is auxiliary variable, so N[k] is a unity column */
         col[k] = -1.0;
      }
      else
      {  /* x[k] is structural variable, so N[k] is a column of the
            original constraint matrix A with negative sign */
         len = glp_get_mat_col(lp, k-m, ind, val);
         for (t = 1; t <= len; t++) col[ind[t]] = val[t];
      }
      /* compute column of the simplex table, which corresponds to the
         specified non-basic variable x[k] */
      glp_ftran(lp, col);
      len = 0;
      for (t = 1; t <= m; t++)
      {  if (col[t] != 0.0)
         {  len++;
            ind[len] = glp_get_bhead(lp, t);
            val[len] = col[t];
         }
      }
      xfree(col);
      /* return to the calling program */
      return len;
}

/***********************************************************************
*  NAME
*
*  glp_transform_row - transform explicitly specified row
*
*  SYNOPSIS
*
*  int glp_transform_row(glp_prob *P, int len, int ind[], double val[]);
*
*  DESCRIPTION
*
*  The routine glp_transform_row performs the same operation as the
*  routine glp_eval_tab_row with exception that the row to be
*  transformed is specified explicitly as a sparse vector.
*
*  The explicitly specified row may be thought as a linear form:
*
*     x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n],             (1)
*
*  where x is an auxiliary variable for this row, a[j] are coefficients
*  of the linear form, x[m+j] are structural variables.
*
*  On entry column indices and numerical values of non-zero elements of
*  the row should be stored in locations ind[1], ..., ind[len] and
*  val[1], ..., val[len], where len is the number of non-zero elements.
*
*  This routine uses the system of equality constraints and the current
*  basis in order to express the auxiliary variable x in (1) through the
*  current non-basic variables (as if the transformed row were added to
*  the problem object and its auxiliary variable were basic), i.e. the
*  resultant row has the form:
*
*     x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n],       (2)
*
*  where xN[j] are non-basic (auxiliary or structural) variables, n is
*  the number of columns in the LP problem object.
*
*  On exit the routine stores indices and numerical values of non-zero
*  elements of the resultant row (2) in locations ind[1], ..., ind[len']
*  and val[1], ..., val[len'], where 0 <= len' <= n is the number of
*  non-zero elements in the resultant row returned by the routine. Note
*  that indices (numbers) of non-basic variables stored in the array ind
*  correspond to original ordinal numbers of variables: indices 1 to m
*  mean auxiliary variables and indices m+1 to m+n mean structural ones.
*
*  RETURNS
*
*  The routine returns len', which is the number of non-zero elements in
*  the resultant row stored in the arrays ind and val.
*
*  BACKGROUND
*
*  The explicitly specified row (1) is transformed in the same way as it
*  were the objective function row.
*
*  From (1) it follows that:
*
*     x = aB * xB + aN * xN,                                         (3)
*
*  where xB is the vector of basic variables, xN is the vector of
*  non-basic variables.
*
*  The simplex table, which corresponds to the current basis, is:
*
*     xB = [-inv(B) * N] * xN.                                       (4)
*
*  Therefore substituting xB from (4) to (3) we have:
*
*     x = aB * [-inv(B) * N] * xN + aN * xN =
*                                                                    (5)
*       = rho * (-N) * xN + aN * xN = alfa * xN,
*
*  where:
*
*     rho = inv(B') * aB,                                            (6)
*
*  and
*
*     alfa = aN + rho * (-N)                                         (7)
*
*  is the resultant row computed by the routine. */

int glp_transform_row(glp_prob *P, int len, int ind[], double val[])
{     int i, j, k, m, n, t, lll, *iii;
      double alfa, *a, *aB, *rho, *vvv;
      if (!glp_bf_exists(P))
         xerror("glp_transform_row: basis factorization does not exist "
            "\n");
      m = glp_get_num_rows(P);
      n = glp_get_num_cols(P);
      /* unpack the row to be transformed to the array a */
      a = xcalloc(1+n, sizeof(double));
      for (j = 1; j <= n; j++) a[j] = 0.0;
      if (!(0 <= len && len <= n))
         xerror("glp_transform_row: len = %d; invalid row length\n",
            len);
      for (t = 1; t <= len; t++)
      {  j = ind[t];
         if (!(1 <= j && j <= n))
            xerror("glp_transform_row: ind[%d] = %d; column index out o"
               "f range\n", t, j);
         if (val[t] == 0.0)
            xerror("glp_transform_row: val[%d] = 0; zero coefficient no"
               "t allowed\n", t);
         if (a[j] != 0.0)
            xerror("glp_transform_row: ind[%d] = %d; duplicate column i"
               "ndices not allowed\n", t, j);
         a[j] = val[t];
      }
      /* construct the vector aB */
      aB = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++)
      {  k = glp_get_bhead(P, i);
         /* xB[i] is k-th original variable */
         xassert(1 <= k && k <= m+n);
         aB[i] = (k <= m ? 0.0 : a[k-m]);
      }
      /* solve the system B'*rho = aB to compute the vector rho */
      rho = aB, glp_btran(P, rho);
      /* compute coefficients at non-basic auxiliary variables */
      len = 0;
      for (i = 1; i <= m; i++)
      {  if (glp_get_row_stat(P, i) != GLP_BS)
         {  alfa = - rho[i];
            if (alfa != 0.0)
            {  len++;
               ind[len] = i;
               val[len] = alfa;
            }
         }
      }
      /* compute coefficients at non-basic structural variables */
      iii = xcalloc(1+m, sizeof(int));
      vvv = xcalloc(1+m, sizeof(double));
      for (j = 1; j <= n; j++)
      {  if (glp_get_col_stat(P, j) != GLP_BS)
         {  alfa = a[j];
            lll = glp_get_mat_col(P, j, iii, vvv);
            for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
            if (alfa != 0.0)
            {  len++;
               ind[len] = m+j;
               val[len] = alfa;
            }
         }
      }
      xassert(len <= n);
      xfree(iii);
      xfree(vvv);
      xfree(aB);
      xfree(a);
      return len;
}

/***********************************************************************
*  NAME
*
*  glp_transform_col - transform explicitly specified column
*
*  SYNOPSIS
*
*  int glp_transform_col(glp_prob *P, int len, int ind[], double val[]);
*
*  DESCRIPTION
*
*  The routine glp_transform_col performs the same operation as the
*  routine glp_eval_tab_col with exception that the column to be
*  transformed is specified explicitly as a sparse vector.
*
*  The explicitly specified column may be thought as if it were added
*  to the original system of equality constraints:
*
*     x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x
*     x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x            (1)
*        .  .  .  .  .  .  .  .  .  .  .  .  .  .  .
*     x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x
*
*  where x[i] are auxiliary variables, x[m+j] are structural variables,
*  x is a structural variable for the explicitly specified column, a[i]
*  are constraint coefficients for x.
*
*  On entry row indices and numerical values of non-zero elements of
*  the column should be stored in locations ind[1], ..., ind[len] and
*  val[1], ..., val[len], where len is the number of non-zero elements.
*
*  This routine uses the system of equality constraints and the current
*  basis in order to express the current basic variables through the
*  structural variable x in (1) (as if the transformed column were added
*  to the problem object and the variable x were non-basic), i.e. the
*  resultant column has the form:
*
*     xB[1] = ... + alfa[1]*x
*     xB[2] = ... + alfa[2]*x                                        (2)
*        .  .  .  .  .  .
*     xB[m] = ... + alfa[m]*x
*
*  where xB are basic (auxiliary and structural) variables, m is the
*  number of rows in the problem object.
*
*  On exit the routine stores indices and numerical values of non-zero
*  elements of the resultant column (2) in locations ind[1], ...,
*  ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the
*  number of non-zero element in the resultant column returned by the
*  routine. Note that indices (numbers) of basic variables stored in
*  the array ind correspond to original ordinal numbers of variables:
*  indices 1 to m mean auxiliary variables and indices m+1 to m+n mean
*  structural ones.
*
*  RETURNS
*
*  The routine returns len', which is the number of non-zero elements
*  in the resultant column stored in the arrays ind and val.
*
*  BACKGROUND
*
*  The explicitly specified column (1) is transformed in the same way
*  as any other column of the constraint matrix using the formula:
*
*     alfa = inv(B) * a,                                             (3)
*
*  where alfa is the resultant column computed by the routine. */

int glp_transform_col(glp_prob *P, int len, int ind[], double val[])
{     int i, m, t;
      double *a, *alfa;
      if (!glp_bf_exists(P))
         xerror("glp_transform_col: basis factorization does not exist "
            "\n");
      m = glp_get_num_rows(P);
      /* unpack the column to be transformed to the array a */
      a = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++) a[i] = 0.0;
      if (!(0 <= len && len <= m))
         xerror("glp_transform_col: len = %d; invalid column length\n",
            len);
      for (t = 1; t <= len; t++)
      {  i = ind[t];
         if (!(1 <= i && i <= m))
            xerror("glp_transform_col: ind[%d] = %d; row index out of r"
               "ange\n", t, i);
         if (val[t] == 0.0)
            xerror("glp_transform_col: val[%d] = 0; zero coefficient no"
               "t allowed\n", t);
         if (a[i] != 0.0)
            xerror("glp_transform_col: ind[%d] = %d; duplicate row indi"
               "ces not allowed\n", t, i);
         a[i] = val[t];
      }
      /* solve the system B*a = alfa to compute the vector alfa */
      alfa = a, glp_ftran(P, alfa);
      /* store resultant coefficients */
      len = 0;
      for (i = 1; i <= m; i++)
      {  if (alfa[i] != 0.0)
         {  len++;
            ind[len] = glp_get_bhead(P, i);
            val[len] = alfa[i];
         }
      }
      xfree(a);
      return len;
}

/***********************************************************************
*  NAME
*
*  glp_prim_rtest - perform primal ratio test
*
*  SYNOPSIS
*
*  int glp_prim_rtest(glp_prob *P, int len, const int ind[],
*     const double val[], int dir, double eps);
*
*  DESCRIPTION
*
*  The routine glp_prim_rtest performs the primal ratio test using an
*  explicitly specified column of the simplex table.
*
*  The current basic solution associated with the LP problem object
*  must be primal feasible.
*
*  The explicitly specified column of the simplex table shows how the
*  basic variables xB depend on some non-basic variable x (which is not
*  necessarily presented in the problem object):
*
*     xB[1] = ... + alfa[1] * x + ...
*     xB[2] = ... + alfa[2] * x + ...                                (*)
*         .  .  .  .  .  .  .  .
*     xB[m] = ... + alfa[m] * x + ...
*
*  The column (*) is specifed on entry to the routine using the sparse
*  format. Ordinal numbers of basic variables xB[i] should be placed in
*  locations ind[1], ..., ind[len], where ordinal number 1 to m denote
*  auxiliary variables, and ordinal numbers m+1 to m+n denote structural
*  variables. The corresponding non-zero coefficients alfa[i] should be
*  placed in locations val[1], ..., val[len]. The arrays ind and val are
*  not changed on exit.
*
*  The parameter dir specifies direction in which the variable x changes
*  on entering the basis: +1 means increasing, -1 means decreasing.
*
*  The parameter eps is an absolute tolerance (small positive number)
*  used by the routine to skip small alfa[j] of the row (*).
*
*  The routine determines which basic variable (among specified in
*  ind[1], ..., ind[len]) should leave the basis in order to keep primal
*  feasibility.
*
*  RETURNS
*
*  The routine glp_prim_rtest returns the index piv in the arrays ind
*  and val corresponding to the pivot element chosen, 1 <= piv <= len.
*  If the adjacent basic solution is primal unbounded and therefore the
*  choice cannot be made, the routine returns zero.
*
*  COMMENTS
*
*  If the non-basic variable x is presented in the LP problem object,
*  the column (*) can be computed with the routine glp_eval_tab_col;
*  otherwise it can be computed with the routine glp_transform_col. */

int glp_prim_rtest(glp_prob *P, int len, const int ind[],
      const double val[], int dir, double eps)
{     int k, m, n, piv, t, type, stat;
      double alfa, big, beta, lb, ub, temp, teta;
      if (glp_get_prim_stat(P) != GLP_FEAS)
         xerror("glp_prim_rtest: basic solution is not primal feasible "
            "\n");
      if (!(dir == +1 || dir == -1))
         xerror("glp_prim_rtest: dir = %d; invalid parameter\n", dir);
      if (!(0.0 < eps && eps < 1.0))
         xerror("glp_prim_rtest: eps = %g; invalid parameter\n", eps);
      m = glp_get_num_rows(P);
      n = glp_get_num_cols(P);
      /* initial settings */
      piv = 0, teta = DBL_MAX, big = 0.0;
      /* walk through the entries of the specified column */
      for (t = 1; t <= len; t++)
      {  /* get the ordinal number of basic variable */
         k = ind[t];
         if (!(1 <= k && k <= m+n))
            xerror("glp_prim_rtest: ind[%d] = %d; variable number out o"
               "f range\n", t, k);
         /* determine type, bounds, status and primal value of basic
            variable xB[i] = x[k] in the current basic solution */
         if (k <= m)
         {  type = glp_get_row_type(P, k);
            lb = glp_get_row_lb(P, k);
            ub = glp_get_row_ub(P, k);
            stat = glp_get_row_stat(P, k);
            beta = glp_get_row_prim(P, k);
         }
         else
         {  type = glp_get_col_type(P, k-m);
            lb = glp_get_col_lb(P, k-m);
            ub = glp_get_col_ub(P, k-m);
            stat = glp_get_col_stat(P, k-m);
            beta = glp_get_col_prim(P, k-m);
         }
         if (stat != GLP_BS)
            xerror("glp_prim_rtest: ind[%d] = %d; non-basic variable no"
               "t allowed\n", t, k);
         /* determine influence coefficient at basic variable xB[i]
            in the explicitly specified column and turn to the case of
            increasing the variable x in order to simplify the program
            logic */
         alfa = (dir > 0 ? + val[t] : - val[t]);
         /* analyze main cases */
         if (type == GLP_FR)
         {  /* xB[i] is free variable */
            continue;
         }
         else if (type == GLP_LO)
lo:      {  /* xB[i] has an lower bound */
            if (alfa > - eps) continue;
            temp = (lb - beta) / alfa;
         }
         else if (type == GLP_UP)
up:      {  /* xB[i] has an upper bound */
            if (alfa < + eps) continue;
            temp = (ub - beta) / alfa;
         }
         else if (type == GLP_DB)
         {  /* xB[i] has both lower and upper bounds */
            if (alfa < 0.0) goto lo; else goto up;
         }
         else if (type == GLP_FX)
         {  /* xB[i] is fixed variable */
            if (- eps < alfa && alfa < + eps) continue;
            temp = 0.0;
         }
         else
            xassert(type != type);
         /* if the value of the variable xB[i] violates its lower or
            upper bound (slightly, because the current basis is assumed
            to be primal feasible), temp is negative; we can think this
            happens due to round-off errors and the value is exactly on
            the bound; this allows replacing temp by zero */
         if (temp < 0.0) temp = 0.0;
         /* apply the minimal ratio test */
         if (teta > temp || teta == temp && big < fabs(alfa))
            piv = t, teta = temp, big = fabs(alfa);
      }
      /* return index of the pivot element chosen */
      return piv;
}

/***********************************************************************
*  NAME
*
*  glp_dual_rtest - perform dual ratio test
*
*  SYNOPSIS
*
*  int glp_dual_rtest(glp_prob *P, int len, const int ind[],
*     const double val[], int dir, double eps);
*
*  DESCRIPTION
*
*  The routine glp_dual_rtest performs the dual ratio test using an
*  explicitly specified row of the simplex table.
*
*  The current basic solution associated with the LP problem object
*  must be dual feasible.
*
*  The explicitly specified row of the simplex table is a linear form
*  that shows how some basic variable x (which is not necessarily
*  presented in the problem object) depends on non-basic variables xN:
*
*     x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*)
*
*  The row (*) is specified on entry to the routine using the sparse
*  format. Ordinal numbers of non-basic variables xN[j] should be placed
*  in locations ind[1], ..., ind[len], where ordinal numbers 1 to m
*  denote auxiliary variables, and ordinal numbers m+1 to m+n denote
*  structural variables. The corresponding non-zero coefficients alfa[j]
*  should be placed in locations val[1], ..., val[len]. The arrays ind
*  and val are not changed on exit.
*
*  The parameter dir specifies direction in which the variable x changes
*  on leaving the basis: +1 means that x goes to its lower bound, and -1
*  means that x goes to its upper bound.
*
*  The parameter eps is an absolute tolerance (small positive number)
*  used by the routine to skip small alfa[j] of the row (*).
*
*  The routine determines which non-basic variable (among specified in
*  ind[1], ..., ind[len]) should enter the basis in order to keep dual
*  feasibility.
*
*  RETURNS
*
*  The routine glp_dual_rtest returns the index piv in the arrays ind
*  and val corresponding to the pivot element chosen, 1 <= piv <= len.
*  If the adjacent basic solution is dual unbounded and therefore the
*  choice cannot be made, the routine returns zero.
*
*  COMMENTS
*
*  If the basic variable x is presented in the LP problem object, the
*  row (*) can be computed with the routine glp_eval_tab_row; otherwise
*  it can be computed with the routine glp_transform_row. */

int glp_dual_rtest(glp_prob *P, int len, const int ind[],
      const double val[], int dir, double eps)
{     int k, m, n, piv, t, stat;
      double alfa, big, cost, obj, temp, teta;
      if (glp_get_dual_stat(P) != GLP_FEAS)
         xerror("glp_dual_rtest: basic solution is not dual feasible\n")
            ;
      if (!(dir == +1 || dir == -1))
         xerror("glp_dual_rtest: dir = %d; invalid parameter\n", dir);
      if (!(0.0 < eps && eps < 1.0))
         xerror("glp_dual_rtest: eps = %g; invalid parameter\n", eps);
      m = glp_get_num_rows(P);
      n = glp_get_num_cols(P);
      /* take into account optimization direction */
      obj = (glp_get_obj_dir(P) == GLP_MIN ? +1.0 : -1.0);
      /* initial settings */
      piv = 0, teta = DBL_MAX, big = 0.0;
      /* walk through the entries of the specified row */
      for (t = 1; t <= len; t++)
      {  /* get ordinal number of non-basic variable */
         k = ind[t];
         if (!(1 <= k && k <= m+n))
            xerror("glp_dual_rtest: ind[%d] = %d; variable number out o"
               "f range\n", t, k);
         /* determine status and reduced cost of non-basic variable
            x[k] = xN[j] in the current basic solution */
         if (k <= m)
         {  stat = glp_get_row_stat(P, k);
            cost = glp_get_row_dual(P, k);
         }
         else
         {  stat = glp_get_col_stat(P, k-m);
            cost = glp_get_col_dual(P, k-m);
         }
         if (stat == GLP_BS)
            xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al"
               "lowed\n", t, k);
         /* determine influence coefficient at non-basic variable xN[j]
            in the explicitly specified row and turn to the case of
            increasing the variable x in order to simplify the program
            logic */
         alfa = (dir > 0 ? + val[t] : - val[t]);
         /* analyze main cases */
         if (stat == GLP_NL)
         {  /* xN[j] is on its lower bound */
            if (alfa < + eps) continue;
            temp = (obj * cost) / alfa;
         }
         else if (stat == GLP_NU)
         {  /* xN[j] is on its upper bound */
            if (alfa > - eps) continue;
            temp = (obj * cost) / alfa;
         }
         else if (stat == GLP_NF)
         {  /* xN[j] is non-basic free variable */
            if (- eps < alfa && alfa < + eps) continue;
            temp = 0.0;
         }
         else if (stat == GLP_NS)
         {  /* xN[j] is non-basic fixed variable */
            continue;
         }
         else
            xassert(stat != stat);
         /* if the reduced cost of the variable xN[j] violates its zero
            bound (slightly, because the current basis is assumed to be
            dual feasible), temp is negative; we can think this happens
            due to round-off errors and the reduced cost is exact zero;
            this allows replacing temp by zero */
         if (temp < 0.0) temp = 0.0;
         /* apply the minimal ratio test */
         if (teta > temp || teta == temp && big < fabs(alfa))
            piv = t, teta = temp, big = fabs(alfa);
      }
      /* return index of the pivot element chosen */
      return piv;
}

/***********************************************************************
*  NAME
*
*  glp_analyze_row - simulate one iteration of dual simplex method
*
*  SYNOPSIS
*
*  int glp_analyze_row(glp_prob *P, int len, const int ind[],
*     const double val[], int type, double rhs, double eps, int *piv,
*     double *x, double *dx, double *y, double *dy, double *dz);
*
*  DESCRIPTION
*
*  Let the current basis be optimal or dual feasible, and there be
*  specified a row (constraint), which is violated by the current basic
*  solution. The routine glp_analyze_row simulates one iteration of the
*  dual simplex method to determine some information on the adjacent
*  basis (see below), where the specified row becomes active constraint
*  (i.e. its auxiliary variable becomes non-basic).
*
*  The current basic solution associated with the problem object passed
*  to the routine must be dual feasible, and its primal components must
*  be defined.
*
*  The row to be analyzed must be previously transformed either with
*  the routine glp_eval_tab_row (if the row is in the problem object)
*  or with the routine glp_transform_row (if the row is external, i.e.
*  not in the problem object). This is needed to express the row only
*  through (auxiliary and structural) variables, which are non-basic in
*  the current basis:
*
*     y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n],
*
*  where y is an auxiliary variable of the row, alfa[j] is an influence
*  coefficient, xN[j] is a non-basic variable.
*
*  The row is passed to the routine in sparse format. Ordinal numbers
*  of non-basic variables are stored in locations ind[1], ..., ind[len],
*  where numbers 1 to m denote auxiliary variables while numbers m+1 to
*  m+n denote structural variables. Corresponding non-zero coefficients
*  alfa[j] are stored in locations val[1], ..., val[len]. The arrays
*  ind and val are ot changed on exit.
*
*  The parameters type and rhs specify the row type and its right-hand
*  side as follows:
*
*     type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs
*
*     type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs
*
*  The parameter eps is an absolute tolerance (small positive number)
*  used by the routine to skip small coefficients alfa[j] on performing
*  the dual ratio test.
*
*  If the operation was successful, the routine stores the following
*  information to corresponding location (if some parameter is NULL,
*  its value is not stored):
*
*  piv   index in the array ind and val, 1 <= piv <= len, determining
*        the non-basic variable, which would enter the adjacent basis;
*
*  x     value of the non-basic variable in the current basis;
*
*  dx    difference between values of the non-basic variable in the
*        adjacent and current bases, dx = x.new - x.old;
*
*  y     value of the row (i.e. of its auxiliary variable) in the
*        current basis;
*
*  dy    difference between values of the row in the adjacent and
*        current bases, dy = y.new - y.old;
*
*  dz    difference between values of the objective function in the
*        adjacent and current bases, dz = z.new - z.old. Note that in
*        case of minimization dz >= 0, and in case of maximization
*        dz <= 0, i.e. in the adjacent basis the objective function
*        always gets worse (degrades). */

int _glp_analyze_row(glp_prob *P, int len, const int ind[],
      const double val[], int type, double rhs, double eps, int *_piv,
      double *_x, double *_dx, double *_y, double *_dy, double *_dz)
{     int t, k, dir, piv, ret = 0;
      double x, dx, y, dy, dz;
      if (P->pbs_stat == GLP_UNDEF)
         xerror("glp_analyze_row: primal basic solution components are "
            "undefined\n");
      if (P->dbs_stat != GLP_FEAS)
         xerror("glp_analyze_row: basic solution is not dual feasible\n"
            );
      /* compute the row value y = sum alfa[j] * xN[j] in the current
         basis */
      if (!(0 <= len && len <= P->n))
         xerror("glp_analyze_row: len = %d; invalid row length\n", len);
      y = 0.0;
      for (t = 1; t <= len; t++)
      {  /* determine value of x[k] = xN[j] in the current basis */
         k = ind[t];
         if (!(1 <= k && k <= P->m+P->n))
            xerror("glp_analyze_row: ind[%d] = %d; row/column index out"
               " of range\n", t, k);
         if (k <= P->m)
         {  /* x[k] is auxiliary variable */
            if (P->row[k]->stat == GLP_BS)
               xerror("glp_analyze_row: ind[%d] = %d; basic auxiliary v"
                  "ariable is not allowed\n", t, k);
            x = P->row[k]->prim;
         }
         else
         {  /* x[k] is structural variable */
            if (P->col[k-P->m]->stat == GLP_BS)
               xerror("glp_analyze_row: ind[%d] = %d; basic structural "
                  "variable is not allowed\n", t, k);
            x = P->col[k-P->m]->prim;
         }
         y += val[t] * x;
      }
      /* check if the row is primal infeasible in the current basis,
         i.e. the constraint is violated at the current point */
      if (type == GLP_LO)
      {  if (y >= rhs)
         {  /* the constraint is not violated */
            ret = 1;
            goto done;
         }
         /* in the adjacent basis y goes to its lower bound */
         dir = +1;
      }
      else if (type == GLP_UP)
      {  if (y <= rhs)
         {  /* the constraint is not violated */
            ret = 1;
            goto done;
         }
         /* in the adjacent basis y goes to its upper bound */
         dir = -1;
      }
      else
         xerror("glp_analyze_row: type = %d; invalid parameter\n",
            type);
      /* compute dy = y.new - y.old */
      dy = rhs - y;
      /* perform dual ratio test to determine which non-basic variable
         should enter the adjacent basis to keep it dual feasible */
      piv = glp_dual_rtest(P, len, ind, val, dir, eps);
      if (piv == 0)
      {  /* no dual feasible adjacent basis exists */
         ret = 2;
         goto done;
      }
      /* non-basic variable x[k] = xN[j] should enter the basis */
      k = ind[piv];
      xassert(1 <= k && k <= P->m+P->n);
      /* determine its value in the current basis */
      if (k <= P->m)
         x = P->row[k]->prim;
      else
         x = P->col[k-P->m]->prim;
      /* compute dx = x.new - x.old = dy / alfa[j] */
      xassert(val[piv] != 0.0);
      dx = dy / val[piv];
      /* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced
         cost of xN[j] in the current basis */
      if (k <= P->m)
         dz = P->row[k]->dual * dx;
      else
         dz = P->col[k-P->m]->dual * dx;
      /* store the analysis results */
      if (_piv != NULL) *_piv = piv;
      if (_x   != NULL) *_x   = x;
      if (_dx  != NULL) *_dx  = dx;
      if (_y   != NULL) *_y   = y;
      if (_dy  != NULL) *_dy  = dy;
      if (_dz  != NULL) *_dz  = dz;
done: return ret;
}

#if 0
int main(void)
{     /* example program for the routine glp_analyze_row */
      glp_prob *P;
      glp_smcp parm;
      int i, k, len, piv, ret, ind[1+100];
      double rhs, x, dx, y, dy, dz, val[1+100];
      P = glp_create_prob();
      /* read plan.mps (see glpk/examples) */
      ret = glp_read_mps(P, GLP_MPS_DECK, NULL, "plan.mps");
      glp_assert(ret == 0);
      /* and solve it to optimality */
      ret = glp_simplex(P, NULL);
      glp_assert(ret == 0);
      glp_assert(glp_get_status(P) == GLP_OPT);
      /* the optimal objective value is 296.217 */
      /* we would like to know what happens if we would add a new row
         (constraint) to plan.mps:
         .01 * bin1 + .01 * bin2 + .02 * bin4 + .02 * bin5 <= 12 */
      /* first, we specify this new row */
      glp_create_index(P);
      len = 0;
      ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
      ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
      ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
      ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
      rhs = 12;
      /* then we can compute value of the row (i.e. of its auxiliary
         variable) in the current basis to see if the constraint is
         violated */
      y = 0.0;
      for (k = 1; k <= len; k++)
         y += val[k] * glp_get_col_prim(P, ind[k]);
      glp_printf("y = %g\n", y);
      /* this prints y = 15.1372, so the constraint is violated, since
         we require that y <= rhs = 12 */
      /* now we transform the row to express it only through non-basic
         (auxiliary and artificial) variables */
      len = glp_transform_row(P, len, ind, val);
      /* finally, we simulate one step of the dual simplex method to
         obtain necessary information for the adjacent basis */
      ret = _glp_analyze_row(P, len, ind, val, GLP_UP, rhs, 1e-9, &piv,
         &x, &dx, &y, &dy, &dz);
      glp_assert(ret == 0);
      glp_printf("k = %d, x = %g; dx = %g; y = %g; dy = %g; dz = %g\n",
         ind[piv], x, dx, y, dy, dz);
      /* this prints dz = 5.64418 and means that in the adjacent basis
         the objective function would be 296.217 + 5.64418 = 301.861 */
      /* now we actually include the row into the problem object; note
         that the arrays ind and val are clobbered, so we need to build
         them once again */
      len = 0;
      ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
      ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
      ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
      ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
      rhs = 12;
      i = glp_add_rows(P, 1);
      glp_set_row_bnds(P, i, GLP_UP, 0, rhs);
      glp_set_mat_row(P, i, len, ind, val);
      /* and perform one dual simplex iteration */
      glp_init_smcp(&parm);
      parm.meth = GLP_DUAL;
      parm.it_lim = 1;
      glp_simplex(P, &parm);
      /* the current objective value is 301.861 */
      return 0;
}
#endif

/***********************************************************************
*  NAME
*
*  glp_analyze_bound - analyze active bound of non-basic variable
*
*  SYNOPSIS
*
*  void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1,
*     double *limit2, int *var2);
*
*  DESCRIPTION
*
*  The routine glp_analyze_bound analyzes the effect of varying the
*  active bound of specified non-basic variable.
*
*  The non-basic variable is specified by the parameter k, where
*  1 <= k <= m means auxiliary variable of corresponding row while
*  m+1 <= k <= m+n means structural variable (column).
*
*  Note that the current basic solution must be optimal, and the basis
*  factorization must exist.
*
*  Results of the analysis have the following meaning.
*
*  value1 is the minimal value of the active bound, at which the basis
*  still remains primal feasible and thus optimal. -DBL_MAX means that
*  the active bound has no lower limit.
*
*  var1 is the ordinal number of an auxiliary (1 to m) or structural
*  (m+1 to n) basic variable, which reaches its bound first and thereby
*  limits further decreasing the active bound being analyzed.
*  if value1 = -DBL_MAX, var1 is set to 0.
*
*  value2 is the maximal value of the active bound, at which the basis
*  still remains primal feasible and thus optimal. +DBL_MAX means that
*  the active bound has no upper limit.
*
*  var2 is the ordinal number of an auxiliary (1 to m) or structural
*  (m+1 to n) basic variable, which reaches its bound first and thereby
*  limits further increasing the active bound being analyzed.
*  if value2 = +DBL_MAX, var2 is set to 0. */

void glp_analyze_bound(glp_prob *P, int k, double *value1, int *var1,
      double *value2, int *var2)
{     GLPROW *row;
      GLPCOL *col;
      int m, n, stat, kase, p, len, piv, *ind;
      double x, new_x, ll, uu, xx, delta, *val;
      /* sanity checks */
      if (P == NULL || P->magic != GLP_PROB_MAGIC)
         xerror("glp_analyze_bound: P = %p; invalid problem object\n",
            P);
      m = P->m, n = P->n;
      if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
         xerror("glp_analyze_bound: optimal basic solution required\n");
      if (!(m == 0 || P->valid))
         xerror("glp_analyze_bound: basis factorization required\n");
      if (!(1 <= k && k <= m+n))
         xerror("glp_analyze_bound: k = %d; variable number out of rang"
            "e\n", k);
      /* retrieve information about the specified non-basic variable
         x[k] whose active bound is to be analyzed */
      if (k <= m)
      {  row = P->row[k];
         stat = row->stat;
         x = row->prim;
      }
      else
      {  col = P->col[k-m];
         stat = col->stat;
         x = col->prim;
      }
      if (stat == GLP_BS)
         xerror("glp_analyze_bound: k = %d; basic variable not allowed "
            "\n", k);
      /* allocate working arrays */
      ind = xcalloc(1+m, sizeof(int));
      val = xcalloc(1+m, sizeof(double));
      /* compute column of the simplex table corresponding to the
         non-basic variable x[k] */
      len = glp_eval_tab_col(P, k, ind, val);
      xassert(0 <= len && len <= m);
      /* perform analysis */
      for (kase = -1; kase <= +1; kase += 2)
      {  /* kase < 0 means active bound of x[k] is decreasing;
            kase > 0 means active bound of x[k] is increasing */
         /* use the primal ratio test to determine some basic variable
            x[p] which reaches its bound first */
         piv = glp_prim_rtest(P, len, ind, val, kase, 1e-9);
         if (piv == 0)
         {  /* nothing limits changing the active bound of x[k] */
            p = 0;
            new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX);
            goto store;
         }
         /* basic variable x[p] limits changing the active bound of
            x[k]; determine its value in the current basis */
         xassert(1 <= piv && piv <= len);
         p = ind[piv];
         if (p <= m)
         {  row = P->row[p];
            ll = glp_get_row_lb(P, row->i);
            uu = glp_get_row_ub(P, row->i);
            stat = row->stat;
            xx = row->prim;
         }
         else
         {  col = P->col[p-m];
            ll = glp_get_col_lb(P, col->j);
            uu = glp_get_col_ub(P, col->j);
            stat = col->stat;
            xx = col->prim;
         }
         xassert(stat == GLP_BS);
         /* determine delta x[p] = bound of x[p] - value of x[p] */
         if (kase < 0 && val[piv] > 0.0 ||
             kase > 0 && val[piv] < 0.0)
         {  /* delta x[p] < 0, so x[p] goes toward its lower bound */
            xassert(ll != -DBL_MAX);
            delta = ll - xx;
         }
         else
         {  /* delta x[p] > 0, so x[p] goes toward its upper bound */
            xassert(uu != +DBL_MAX);
            delta = uu - xx;
         }
         /* delta x[p] = alfa[p,k] * delta x[k], so new x[k] = x[k] +
            delta x[k] = x[k] + delta x[p] / alfa[p,k] is the value of
            x[k] in the adjacent basis */
         xassert(val[piv] != 0.0);
         new_x = x + delta / val[piv];
store:   /* store analysis results */
         if (kase < 0)
         {  if (value1 != NULL) *value1 = new_x;
            if (var1 != NULL) *var1 = p;
         }
         else
         {  if (value2 != NULL) *value2 = new_x;
            if (var2 != NULL) *var2 = p;
         }
      }
      /* free working arrays */
      xfree(ind);
      xfree(val);
      return;
}

/***********************************************************************
*  NAME
*
*  glp_analyze_coef - analyze objective coefficient at basic variable
*
*  SYNOPSIS
*
*  void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
*     double *value1, double *coef2, int *var2, double *value2);
*
*  DESCRIPTION
*
*  The routine glp_analyze_coef analyzes the effect of varying the
*  objective coefficient at specified basic variable.
*
*  The basic variable is specified by the parameter k, where
*  1 <= k <= m means auxiliary variable of corresponding row while
*  m+1 <= k <= m+n means structural variable (column).
*
*  Note that the current basic solution must be optimal, and the basis
*  factorization must exist.
*
*  Results of the analysis have the following meaning.
*
*  coef1 is the minimal value of the objective coefficient, at which
*  the basis still remains dual feasible and thus optimal. -DBL_MAX
*  means that the objective coefficient has no lower limit.
*
*  var1 is the ordinal number of an auxiliary (1 to m) or structural
*  (m+1 to n) non-basic variable, whose reduced cost reaches its zero
*  bound first and thereby limits further decreasing the objective
*  coefficient being analyzed. If coef1 = -DBL_MAX, var1 is set to 0.
*
*  value1 is value of the basic variable being analyzed in an adjacent
*  basis, which is defined as follows. Let the objective coefficient
*  reaches its minimal value (coef1) and continues decreasing. Then the
*  reduced cost of the limiting non-basic variable (var1) becomes dual
*  infeasible and the current basis becomes non-optimal that forces the
*  limiting non-basic variable to enter the basis replacing there some
*  basic variable that leaves the basis to keep primal feasibility.
*  Should note that on determining the adjacent basis current bounds
*  of the basic variable being analyzed are ignored as if it were free
*  (unbounded) variable, so it cannot leave the basis. It may happen
*  that no dual feasible adjacent basis exists, in which case value1 is
*  set to -DBL_MAX or +DBL_MAX.
*
*  coef2 is the maximal value of the objective coefficient, at which
*  the basis still remains dual feasible and thus optimal. +DBL_MAX
*  means that the objective coefficient has no upper limit.
*
*  var2 is the ordinal number of an auxiliary (1 to m) or structural
*  (m+1 to n) non-basic variable, whose reduced cost reaches its zero
*  bound first and thereby limits further increasing the objective
*  coefficient being analyzed. If coef2 = +DBL_MAX, var2 is set to 0.
*
*  value2 is value of the basic variable being analyzed in an adjacent
*  basis, which is defined exactly in the same way as value1 above with
*  exception that now the objective coefficient is increasing. */

void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
      double *value1, double *coef2, int *var2, double *value2)
{     GLPROW *row; GLPCOL *col;
      int m, n, type, stat, kase, p, q, dir, clen, cpiv, rlen, rpiv,
         *cind, *rind;
      double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx,
         *rval, *cval;
      /* sanity checks */
      if (P == NULL || P->magic != GLP_PROB_MAGIC)
         xerror("glp_analyze_coef: P = %p; invalid problem object\n",
            P);
      m = P->m, n = P->n;
      if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
         xerror("glp_analyze_coef: optimal basic solution required\n");
      if (!(m == 0 || P->valid))
         xerror("glp_analyze_coef: basis factorization required\n");
      if (!(1 <= k && k <= m+n))
         xerror("glp_analyze_coef: k = %d; variable number out of range"
            "\n", k);
      /* retrieve information about the specified basic variable x[k]
         whose objective coefficient c[k] is to be analyzed */
      if (k <= m)
      {  row = P->row[k];
         type = row->type;
         lb = row->lb;
         ub = row->ub;
         coef = 0.0;
         stat = row->stat;
         x = row->prim;
      }
      else
      {  col = P->col[k-m];
         type = col->type;
         lb = col->lb;
         ub = col->ub;
         coef = col->coef;
         stat = col->stat;
         x = col->prim;
      }
      if (stat != GLP_BS)
         xerror("glp_analyze_coef: k = %d; non-basic variable not allow"
            "ed\n", k);
      /* allocate working arrays */
      cind = xcalloc(1+m, sizeof(int));
      cval = xcalloc(1+m, sizeof(double));
      rind = xcalloc(1+n, sizeof(int));
      rval = xcalloc(1+n, sizeof(double));
      /* compute row of the simplex table corresponding to the basic
         variable x[k] */
      rlen = glp_eval_tab_row(P, k, rind, rval);
      xassert(0 <= rlen && rlen <= n);
      /* perform analysis */
      for (kase = -1; kase <= +1; kase += 2)
      {  /* kase < 0 means objective coefficient c[k] is decreasing;
            kase > 0 means objective coefficient c[k] is increasing */
         /* note that decreasing c[k] is equivalent to increasing dual
            variable lambda[k] and vice versa; we need to correctly set
            the dir flag as required by the routine glp_dual_rtest */
         if (P->dir == GLP_MIN)
            dir = - kase;
         else if (P->dir == GLP_MAX)
            dir = + kase;
         else
            xassert(P != P);
         /* use the dual ratio test to determine non-basic variable
            x[q] whose reduced cost d[q] reaches zero bound first */
         rpiv = glp_dual_rtest(P, rlen, rind, rval, dir, 1e-9);
         if (rpiv == 0)
         {  /* nothing limits changing c[k] */
            lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX);
            q = 0;
            /* x[k] keeps its current value */
            new_x = x;
            goto store;
         }
         /* non-basic variable x[q] limits changing coefficient c[k];
            determine its status and reduced cost d[k] in the current
            basis */
         xassert(1 <= rpiv && rpiv <= rlen);
         q = rind[rpiv];
         xassert(1 <= q && q <= m+n);
         if (q <= m)
         {  row = P->row[q];
            stat = row->stat;
            d = row->dual;
         }
         else
         {  col = P->col[q-m];
            stat = col->stat;
            d = col->dual;
         }
         /* note that delta d[q] = new d[q] - d[q] = - d[q], because
            new d[q] = 0; delta d[q] = alfa[k,q] * delta c[k], so
            delta c[k] = delta d[q] / alfa[k,q] = - d[q] / alfa[k,q] */
         xassert(rval[rpiv] != 0.0);
         delta = - d / rval[rpiv];
         /* compute new c[k] = c[k] + delta c[k], which is the limiting
            value of the objective coefficient c[k] */
         lim_coef = coef + delta;
         /* let c[k] continue decreasing/increasing that makes d[q]
            dual infeasible and forces x[q] to enter the basis;
            to perform the primal ratio test we need to know in which
            direction x[q] changes on entering the basis; we determine
            that analyzing the sign of delta d[q] (see above), since
            d[q] may be close to zero having wrong sign */
         /* let, for simplicity, the problem is minimization */
         if (kase < 0 && rval[rpiv] > 0.0 ||
             kase > 0 && rval[rpiv] < 0.0)
         {  /* delta d[q] < 0, so d[q] being non-negative will become
               negative, so x[q] will increase */
            dir = +1;
         }
         else
         {  /* delta d[q] > 0, so d[q] being non-positive will become
               positive, so x[q] will decrease */
            dir = -1;
         }
         /* if the problem is maximization, correct the direction */
         if (P->dir == GLP_MAX) dir = - dir;
         /* check that we didn't make a silly mistake */
         if (dir > 0)
            xassert(stat == GLP_NL || stat == GLP_NF);
         else
            xassert(stat == GLP_NU || stat == GLP_NF);
         /* compute column of the simplex table corresponding to the
            non-basic variable x[q] */
         clen = glp_eval_tab_col(P, q, cind, cval);
         /* make x[k] temporarily free (unbounded) */
         if (k <= m)
         {  row = P->row[k];
            row->type = GLP_FR;
            row->lb = row->ub = 0.0;
         }
         else
         {  col = P->col[k-m];
            col->type = GLP_FR;
            col->lb = col->ub = 0.0;
         }
         /* use the primal ratio test to determine some basic variable
            which leaves the basis */
         cpiv = glp_prim_rtest(P, clen, cind, cval, dir, 1e-9);
         /* restore original bounds of the basic variable x[k] */
         if (k <= m)
         {  row = P->row[k];
            row->type = type;
            row->lb = lb, row->ub = ub;
         }
         else
         {  col = P->col[k-m];
            col->type = type;
            col->lb = lb, col->ub = ub;
         }
         if (cpiv == 0)
         {  /* non-basic variable x[q] can change unlimitedly */
            if (dir < 0 && rval[rpiv] > 0.0 ||
                dir > 0 && rval[rpiv] < 0.0)
            {  /* delta x[k] = alfa[k,q] * delta x[q] < 0 */
               new_x = -DBL_MAX;
            }
            else
            {  /* delta x[k] = alfa[k,q] * delta x[q] > 0 */
               new_x = +DBL_MAX;
            }
            goto store;
         }
         /* some basic variable x[p] limits changing non-basic variable
            x[q] in the adjacent basis */
         xassert(1 <= cpiv && cpiv <= clen);
         p = cind[cpiv];
         xassert(1 <= p && p <= m+n);
         xassert(p != k);
         if (p <= m)
         {  row = P->row[p];
            xassert(row->stat == GLP_BS);
            ll = glp_get_row_lb(P, row->i);
            uu = glp_get_row_ub(P, row->i);
            xx = row->prim;
         }
         else
         {  col = P->col[p-m];
            xassert(col->stat == GLP_BS);
            ll = glp_get_col_lb(P, col->j);
            uu = glp_get_col_ub(P, col->j);
            xx = col->prim;
         }
         /* determine delta x[p] = new x[p] - x[p] */
         if (dir < 0 && cval[cpiv] > 0.0 ||
             dir > 0 && cval[cpiv] < 0.0)
         {  /* delta x[p] < 0, so x[p] goes toward its lower bound */
            xassert(ll != -DBL_MAX);
            delta = ll - xx;
         }
         else
         {  /* delta x[p] > 0, so x[p] goes toward its upper bound */
            xassert(uu != +DBL_MAX);
            delta = uu - xx;
         }
         /* compute new x[k] = x[k] + alfa[k,q] * delta x[q], where
            delta x[q] = delta x[p] / alfa[p,q] */
         xassert(cval[cpiv] != 0.0);
         new_x = x + (rval[rpiv] / cval[cpiv]) * delta;
store:   /* store analysis results */
         if (kase < 0)
         {  if (coef1 != NULL) *coef1 = lim_coef;
            if (var1 != NULL) *var1 = q;
            if (value1 != NULL) *value1 = new_x;
         }
         else
         {  if (coef2 != NULL) *coef2 = lim_coef;
            if (var2 != NULL) *var2 = q;
            if (value2 != NULL) *value2 = new_x;
         }
      }
      /* free working arrays */
      xfree(cind);
      xfree(cval);
      xfree(rind);
      xfree(rval);
      return;
}

/* eof */