text: /* glpapi11.c (basis factorization and simplex tableau routines) */
  
  /***********************************************************************
  *  This code is part of GLPK (GNU Linear Programming Kit).
  *
  *  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
  *  Department for Applied Informatics, Moscow Aviation Institute,
  *  Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
  *
  *  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/>.
  ***********************************************************************/
  
  #define _GLPBFD_PRIVATE
  #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");
        cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt);
        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 = 50;
           parm->upd_tol = 1e-6;
           parm->nrs_max = 50;
           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. */
  
  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;
  }
  
  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_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;
  }
  
  /* eof */
  