/* 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 */