/* glpspx01.c (primal simplex method) */

/***********************************************************************

*  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 "glpspx.h"

struct csa

{     /* common storage area */

      /*--------------------------------------------------------------*/

      /* LP data */

      int m;

      /* number of rows (auxiliary variables), m > 0 */

      int n;

      /* number of columns (structural variables), n > 0 */

      char *type; /* char type[1+m+n]; */

      /* type[0] is not used;

         type[k], 1 <= k <= m+n, is the type of variable x[k]:

         GLP_FR - free variable

         GLP_LO - variable with lower bound

         GLP_UP - variable with upper bound

         GLP_DB - double-bounded variable

         GLP_FX - fixed variable */

      double *lb; /* double lb[1+m+n]; */

      /* lb[0] is not used;

         lb[k], 1 <= k <= m+n, is an lower bound of variable x[k];

         if x[k] has no lower bound, lb[k] is zero */

      double *ub; /* double ub[1+m+n]; */

      /* ub[0] is not used;

         ub[k], 1 <= k <= m+n, is an upper bound of variable x[k];

         if x[k] has no upper bound, ub[k] is zero;

         if x[k] is of fixed type, ub[k] is the same as lb[k] */

      double *coef; /* double coef[1+m+n]; */

      /* coef[0] is not used;

         coef[k], 1 <= k <= m+n, is an objective coefficient at

         variable x[k] (note that on phase I auxiliary variables also

         may have non-zero objective coefficients) */

      /*--------------------------------------------------------------*/

      /* original objective function */

      double *obj; /* double obj[1+n]; */

      /* obj[0] is a constant term of the original objective function;

         obj[j], 1 <= j <= n, is an original objective coefficient at

         structural variable x[m+j] */

      double zeta;

      /* factor used to scale original objective coefficients; its

         sign defines original optimization direction: zeta > 0 means

         minimization, zeta < 0 means maximization */

      /*--------------------------------------------------------------*/

      /* constraint matrix A; it has m rows and n columns and is stored

         by columns */

      int *A_ptr; /* int A_ptr[1+n+1]; */

      /* A_ptr[0] is not used;

         A_ptr[j], 1 <= j <= n, is starting position of j-th column in

         arrays A_ind and A_val; note that A_ptr[1] is always 1;

         A_ptr[n+1] indicates the position after the last element in

         arrays A_ind and A_val */

      int *A_ind; /* int A_ind[A_ptr[n+1]]; */

      /* row indices */

      double *A_val; /* double A_val[A_ptr[n+1]]; */

      /* non-zero element values */

      /*--------------------------------------------------------------*/

      /* basis header */

      int *head; /* int head[1+m+n]; */

      /* head[0] is not used;

         head[i], 1 <= i <= m, is the ordinal number of basic variable

         xB[i]; head[i] = k means that xB[i] = x[k] and i-th column of

         matrix B is k-th column of matrix (I|-A);

         head[m+j], 1 <= j <= n, is the ordinal number of non-basic

         variable xN[j]; head[m+j] = k means that xN[j] = x[k] and j-th

         column of matrix N is k-th column of matrix (I|-A) */

      char *stat; /* char stat[1+n]; */

      /* stat[0] is not used;

         stat[j], 1 <= j <= n, is the status of non-basic variable

         xN[j], which defines its active bound:

         GLP_NL - lower bound is active

         GLP_NU - upper bound is active

         GLP_NF - free variable

         GLP_NS - fixed variable */

      /*--------------------------------------------------------------*/

      /* matrix B is the basis matrix; it is composed from columns of

         the augmented constraint matrix (I|-A) corresponding to basic

         variables and stored in a factorized (invertable) form */

      int valid;

      /* factorization is valid only if this flag is set */

      BFD *bfd; /* BFD bfd[1:m,1:m]; */

      /* factorized (invertable) form of the basis matrix */

      /*--------------------------------------------------------------*/

      /* matrix N is a matrix composed from columns of the augmented

         constraint matrix (I|-A) corresponding to non-basic variables

         except fixed ones; it is stored by rows and changes every time

         the basis changes */

      int *N_ptr; /* int N_ptr[1+m+1]; */

      /* N_ptr[0] is not used;

         N_ptr[i], 1 <= i <= m, is starting position of i-th row in

         arrays N_ind and N_val; note that N_ptr[1] is always 1;

         N_ptr[m+1] indicates the position after the last element in

         arrays N_ind and N_val */

      int *N_len; /* int N_len[1+m]; */

      /* N_len[0] is not used;

         N_len[i], 1 <= i <= m, is length of i-th row (0 to n) */

      int *N_ind; /* int N_ind[N_ptr[m+1]]; */

      /* column indices */

      double *N_val; /* double N_val[N_ptr[m+1]]; */

      /* non-zero element values */

      /*--------------------------------------------------------------*/

      /* working parameters */

      int phase;

      /* search phase:

         0 - not determined yet

         1 - search for primal feasible solution

         2 - search for optimal solution */

      glp_long tm_beg;

      /* time value at the beginning of the search */

      int it_beg;

      /* simplex iteration count at the beginning of the search */

      int it_cnt;

      /* simplex iteration count; it increases by one every time the

         basis changes (including the case when a non-basic variable

         jumps to its opposite bound) */

      int it_dpy;

      /* simplex iteration count at the most recent display output */

      /*--------------------------------------------------------------*/

      /* basic solution components */

      double *bbar; /* double bbar[1+m]; */

      /* bbar[0] is not used;

         bbar[i], 1 <= i <= m, is primal value of basic variable xB[i]

         (if xB[i] is free, its primal value is not updated) */

      double *cbar; /* double cbar[1+n]; */

      /* cbar[0] is not used;

         cbar[j], 1 <= j <= n, is reduced cost of non-basic variable

         xN[j] (if xN[j] is fixed, its reduced cost is not updated) */

      /*--------------------------------------------------------------*/

      /* the following pricing technique options may be used:

         GLP_PT_STD - standard ("textbook") pricing;

         GLP_PT_PSE - projected steepest edge;

         GLP_PT_DVX - Devex pricing (not implemented yet);

         in case of GLP_PT_STD the reference space is not used, and all

         steepest edge coefficients are set to 1 */

      int refct;

      /* this count is set to an initial value when the reference space

         is defined and decreases by one every time the basis changes;

         once this count reaches zero, the reference space is redefined

         again */

      char *refsp; /* char refsp[1+m+n]; */

      /* refsp[0] is not used;

         refsp[k], 1 <= k <= m+n, is the flag which means that variable

         x[k] belongs to the current reference space */

      double *gamma; /* double gamma[1+n]; */

      /* gamma[0] is not used;

         gamma[j], 1 <= j <= n, is the steepest edge coefficient for

         non-basic variable xN[j]; if xN[j] is fixed, gamma[j] is not

         used and just set to 1 */

      /*--------------------------------------------------------------*/

      /* non-basic variable xN[q] chosen to enter the basis */

      int q;

      /* index of the non-basic variable xN[q] chosen, 1 <= q <= n;

         if the set of eligible non-basic variables is empty and thus

         no variable has been chosen, q is set to 0 */

      /*--------------------------------------------------------------*/

      /* pivot column of the simplex table corresponding to non-basic

         variable xN[q] chosen is the following vector:

            T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],

         where B is the current basis matrix, N[q] is a column of the

         matrix (I|-A) corresponding to xN[q] */

      int tcol_nnz;

      /* number of non-zero components, 0 <= nnz <= m */

      int *tcol_ind; /* int tcol_ind[1+m]; */

      /* tcol_ind[0] is not used;

         tcol_ind[t], 1 <= t <= nnz, is an index of non-zero component,

         i.e. tcol_ind[t] = i means that tcol_vec[i] != 0 */

      double *tcol_vec; /* double tcol_vec[1+m]; */

      /* tcol_vec[0] is not used;

         tcol_vec[i], 1 <= i <= m, is a numeric value of i-th component

         of the column */

      double tcol_max;

      /* infinity (maximum) norm of the column (max |tcol_vec[i]|) */

      int tcol_num;

      /* number of significant non-zero components, which means that:

         |tcol_vec[i]| >= eps for i in tcol_ind[1,...,num],

         |tcol_vec[i]| <  eps for i in tcol_ind[num+1,...,nnz],

         where eps is a pivot tolerance */

      /*--------------------------------------------------------------*/

      /* basic variable xB[p] chosen to leave the basis */

      int p;

      /* index of the basic variable xB[p] chosen, 1 <= p <= m;

         p = 0 means that no basic variable reaches its bound;

         p < 0 means that non-basic variable xN[q] reaches its opposite

         bound before any basic variable */

      int p_stat;

      /* new status (GLP_NL, GLP_NU, or GLP_NS) to be assigned to xB[p]

         once it has left the basis */

      double teta;

      /* change of non-basic variable xN[q] (see above), on which xB[p]

         (or, if p < 0, xN[q] itself) reaches its bound */

      /*--------------------------------------------------------------*/

      /* pivot row of the simplex table corresponding to basic variable

         xB[p] chosen is the following vector:

            T' * e[p] = - N' * inv(B') * e[p] = - N' * rho,

         where B' is a matrix transposed to the current basis matrix,

         N' is a matrix, whose rows are columns of the matrix (I|-A)

         corresponding to non-basic non-fixed variables */

      int trow_nnz;

      /* number of non-zero components, 0 <= nnz <= n */

      int *trow_ind; /* int trow_ind[1+n]; */

      /* trow_ind[0] is not used;

         trow_ind[t], 1 <= t <= nnz, is an index of non-zero component,

         i.e. trow_ind[t] = j means that trow_vec[j] != 0 */

      double *trow_vec; /* int trow_vec[1+n]; */

      /* trow_vec[0] is not used;

         trow_vec[j], 1 <= j <= n, is a numeric value of j-th component

         of the row */

      /*--------------------------------------------------------------*/

      /* working arrays */

      double *work1; /* double work1[1+m]; */

      double *work2; /* double work2[1+m]; */

      double *work3; /* double work3[1+m]; */

      double *work4; /* double work4[1+m]; */

};

static const double kappa = 0.10;

/***********************************************************************

*  alloc_csa - allocate common storage area

*

*  This routine allocates all arrays in the common storage area (CSA)

*  and returns a pointer to the CSA. */

static struct csa *alloc_csa(glp_prob *lp)

{     struct csa *csa;

      int m = lp->m;

      int n = lp->n;

      int nnz = lp->nnz;

      csa = xmalloc(sizeof(struct csa));

      xassert(m > 0 && n > 0);

      csa->m = m;

      csa->n = n;

      csa->type = xcalloc(1+m+n, sizeof(char));

      csa->lb = xcalloc(1+m+n, sizeof(double));

      csa->ub = xcalloc(1+m+n, sizeof(double));

      csa->coef = xcalloc(1+m+n, sizeof(double));

      csa->obj = xcalloc(1+n, sizeof(double));

      csa->A_ptr = xcalloc(1+n+1, sizeof(int));

      csa->A_ind = xcalloc(1+nnz, sizeof(int));

      csa->A_val = xcalloc(1+nnz, sizeof(double));

      csa->head = xcalloc(1+m+n, sizeof(int));

      csa->stat = xcalloc(1+n, sizeof(char));

      csa->N_ptr = xcalloc(1+m+1, sizeof(int));

      csa->N_len = xcalloc(1+m, sizeof(int));

      csa->N_ind = NULL; /* will be allocated later */

      csa->N_val = NULL; /* will be allocated later */

      csa->bbar = xcalloc(1+m, sizeof(double));

      csa->cbar = xcalloc(1+n, sizeof(double));

      csa->refsp = xcalloc(1+m+n, sizeof(char));

      csa->gamma = xcalloc(1+n, sizeof(double));

      csa->tcol_ind = xcalloc(1+m, sizeof(int));

      csa->tcol_vec = xcalloc(1+m, sizeof(double));

      csa->trow_ind = xcalloc(1+n, sizeof(int));

      csa->trow_vec = xcalloc(1+n, sizeof(double));

      csa->work1 = xcalloc(1+m, sizeof(double));

      csa->work2 = xcalloc(1+m, sizeof(double));

      csa->work3 = xcalloc(1+m, sizeof(double));

      csa->work4 = xcalloc(1+m, sizeof(double));

      return csa;

}

/***********************************************************************

*  init_csa - initialize common storage area

*

*  This routine initializes all data structures in the common storage

*  area (CSA). */

static void alloc_N(struct csa *csa);

static void build_N(struct csa *csa);

static void init_csa(struct csa *csa, glp_prob *lp)

{     int m = csa->m;

      int n = csa->n;

      char *type = csa->type;

      double *lb = csa->lb;

      double *ub = csa->ub;

      double *coef = csa->coef;

      double *obj = csa->obj;

      int *A_ptr = csa->A_ptr;

      int *A_ind = csa->A_ind;

      double *A_val = csa->A_val;

      int *head = csa->head;

      char *stat = csa->stat;

      char *refsp = csa->refsp;

      double *gamma = csa->gamma;

      int i, j, k, loc;

      double cmax;

      /* auxiliary variables */

      for (i = 1; i <= m; i++)

      {  GLPROW *row = lp->row[i];

         type[i] = (char)row->type;

         lb[i] = row->lb * row->rii;

         ub[i] = row->ub * row->rii;

         coef[i] = 0.0;

      }

      /* structural variables */

      for (j = 1; j <= n; j++)

      {  GLPCOL *col = lp->col[j];

         type[m+j] = (char)col->type;

         lb[m+j] = col->lb / col->sjj;

         ub[m+j] = col->ub / col->sjj;

         coef[m+j] = col->coef * col->sjj;

      }

      /* original objective function */

      obj[0] = lp->c0;

      memcpy(&obj[1], &coef[m+1], n * sizeof(double));

      /* factor used to scale original objective coefficients */

      cmax = 0.0;

      for (j = 1; j <= n; j++)

         if (cmax < fabs(obj[j])) cmax = fabs(obj[j]);

      if (cmax == 0.0) cmax = 1.0;

      switch (lp->dir)

      {  case GLP_MIN:

            csa->zeta = + 1.0 / cmax;

            break;

         case GLP_MAX:

            csa->zeta = - 1.0 / cmax;

            break;

         default:

            xassert(lp != lp);

      }

#if 1

      if (fabs(csa->zeta) < 1.0) csa->zeta *= 1000.0;

#endif

      /* matrix A (by columns) */

      loc = 1;

      for (j = 1; j <= n; j++)

      {  GLPAIJ *aij;

         A_ptr[j] = loc;

         for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)

         {  A_ind[loc] = aij->row->i;

            A_val[loc] = aij->row->rii * aij->val * aij->col->sjj;

            loc++;

         }

      }

      A_ptr[n+1] = loc;

      xassert(loc == lp->nnz+1);

      /* basis header */

      xassert(lp->valid);

      memcpy(&head[1], &lp->head[1], m * sizeof(int));

      k = 0;

      for (i = 1; i <= m; i++)

      {  GLPROW *row = lp->row[i];

         if (row->stat != GLP_BS)

         {  k++;

            xassert(k <= n);

            head[m+k] = i;

            stat[k] = (char)row->stat;

         }

      }

      for (j = 1; j <= n; j++)

      {  GLPCOL *col = lp->col[j];

         if (col->stat != GLP_BS)

         {  k++;

            xassert(k <= n);

            head[m+k] = m + j;

            stat[k] = (char)col->stat;

         }

      }

      xassert(k == n);

      /* factorization of matrix B */

      csa->valid = 1, lp->valid = 0;

      csa->bfd = lp->bfd, lp->bfd = NULL;

      /* matrix N (by rows) */

      alloc_N(csa);

      build_N(csa);

      /* working parameters */

      csa->phase = 0;

      csa->tm_beg = xtime();

      csa->it_beg = csa->it_cnt = lp->it_cnt;

      csa->it_dpy = -1;

      /* reference space and steepest edge coefficients */

      csa->refct = 0;

      memset(&refsp[1], 0, (m+n) * sizeof(char));

      for (j = 1; j <= n; j++) gamma[j] = 1.0;

      return;

}

/***********************************************************************

*  invert_B - compute factorization of the basis matrix

*

*  This routine computes factorization of the current basis matrix B.

*

*  If the operation is successful, the routine returns zero, otherwise

*  non-zero. */

static int inv_col(void *info, int i, int ind[], double val[])

{     /* this auxiliary routine returns row indices and numeric values

         of non-zero elements of i-th column of the basis matrix */

      struct csa *csa = info;

      int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      int *A_ptr = csa->A_ptr;

      int *A_ind = csa->A_ind;

      double *A_val = csa->A_val;

      int *head = csa->head;

      int k, len, ptr, t;

#ifdef GLP_DEBUG

      xassert(1 <= i && i <= m);

#endif

      k = head[i]; /* B[i] is k-th column of (I|-A) */

#ifdef GLP_DEBUG

      xassert(1 <= k && k <= m+n);

#endif

      if (k <= m)

      {  /* B[i] is k-th column of submatrix I */

         len = 1;

         ind[1] = k;

         val[1] = 1.0;

      }

      else

      {  /* B[i] is (k-m)-th column of submatrix (-A) */

         ptr = A_ptr[k-m];

         len = A_ptr[k-m+1] - ptr;

         memcpy(&ind[1], &A_ind[ptr], len * sizeof(int));

         memcpy(&val[1], &A_val[ptr], len * sizeof(double));

         for (t = 1; t <= len; t++) val[t] = - val[t];

      }

      return len;

}

static int invert_B(struct csa *csa)

{     int ret;

      ret = bfd_factorize(csa->bfd, csa->m, NULL, inv_col, csa);

      csa->valid = (ret == 0);

      return ret;

}

/***********************************************************************

*  update_B - update factorization of the basis matrix

*

*  This routine replaces i-th column of the basis matrix B by k-th

*  column of the augmented constraint matrix (I|-A) and then updates

*  the factorization of B.

*

*  If the factorization has been successfully updated, the routine

*  returns zero, otherwise non-zero. */

static int update_B(struct csa *csa, int i, int k)

{     int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      int ret;

#ifdef GLP_DEBUG

      xassert(1 <= i && i <= m);

      xassert(1 <= k && k <= m+n);

#endif

      if (k <= m)

      {  /* new i-th column of B is k-th column of I */

         int ind[1+1];

         double val[1+1];

         ind[1] = k;

         val[1] = 1.0;

         xassert(csa->valid);

         ret = bfd_update_it(csa->bfd, i, 0, 1, ind, val);

      }

      else

      {  /* new i-th column of B is (k-m)-th column of (-A) */

         int *A_ptr = csa->A_ptr;

         int *A_ind = csa->A_ind;

         double *A_val = csa->A_val;

         double *val = csa->work1;

         int beg, end, ptr, len;

         beg = A_ptr[k-m];

         end = A_ptr[k-m+1];

         len = 0;

         for (ptr = beg; ptr < end; ptr++)

            val[++len] = - A_val[ptr];

         xassert(csa->valid);

         ret = bfd_update_it(csa->bfd, i, 0, len, &A_ind[beg-1], val);

      }

      csa->valid = (ret == 0);

      return ret;

}

/***********************************************************************

*  error_ftran - compute residual vector r = h - B * x

*

*  This routine computes the residual vector r = h - B * x, where B is

*  the current basis matrix, h is the vector of right-hand sides, x is

*  the solution vector. */

static void error_ftran(struct csa *csa, double h[], double x[],

      double r[])

{     int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      int *A_ptr = csa->A_ptr;

      int *A_ind = csa->A_ind;

      double *A_val = csa->A_val;

      int *head = csa->head;

      int i, k, beg, end, ptr;

      double temp;

      /* compute the residual vector:

         r = h - B * x = h - B[1] * x[1] - ... - B[m] * x[m],

         where B[1], ..., B[m] are columns of matrix B */

      memcpy(&r[1], &h[1], m * sizeof(double));

      for (i = 1; i <= m; i++)

      {  temp = x[i];

         if (temp == 0.0) continue;

         k = head[i]; /* B[i] is k-th column of (I|-A) */

#ifdef GLP_DEBUG

         xassert(1 <= k && k <= m+n);

#endif

         if (k <= m)

         {  /* B[i] is k-th column of submatrix I */

            r[k] -= temp;

         }

         else

         {  /* B[i] is (k-m)-th column of submatrix (-A) */

            beg = A_ptr[k-m];

            end = A_ptr[k-m+1];

            for (ptr = beg; ptr < end; ptr++)

               r[A_ind[ptr]] += A_val[ptr] * temp;

         }

      }

      return;

}

/***********************************************************************

*  refine_ftran - refine solution of B * x = h

*

*  This routine performs one iteration to refine the solution of

*  the system B * x = h, where B is the current basis matrix, h is the

*  vector of right-hand sides, x is the solution vector. */

static void refine_ftran(struct csa *csa, double h[], double x[])

{     int m = csa->m;

      double *r = csa->work1;

      double *d = csa->work1;

      int i;

      /* compute the residual vector r = h - B * x */

      error_ftran(csa, h, x, r);

      /* compute the correction vector d = inv(B) * r */

      xassert(csa->valid);

      bfd_ftran(csa->bfd, d);

      /* refine the solution vector (new x) = (old x) + d */

      for (i = 1; i <= m; i++) x[i] += d[i];

      return;

}

/***********************************************************************

*  error_btran - compute residual vector r = h - B'* x

*

*  This routine computes the residual vector r = h - B'* x, where B'

*  is a matrix transposed to the current basis matrix, h is the vector

*  of right-hand sides, x is the solution vector. */

static void error_btran(struct csa *csa, double h[], double x[],

      double r[])

{     int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      int *A_ptr = csa->A_ptr;

      int *A_ind = csa->A_ind;

      double *A_val = csa->A_val;

      int *head = csa->head;

      int i, k, beg, end, ptr;

      double temp;

      /* compute the residual vector r = b - B'* x */

      for (i = 1; i <= m; i++)

      {  /* r[i] := b[i] - (i-th column of B)'* x */

         k = head[i]; /* B[i] is k-th column of (I|-A) */

#ifdef GLP_DEBUG

         xassert(1 <= k && k <= m+n);

#endif

         temp = h[i];

         if (k <= m)

         {  /* B[i] is k-th column of submatrix I */

            temp -= x[k];

         }

         else

         {  /* B[i] is (k-m)-th column of submatrix (-A) */

            beg = A_ptr[k-m];

            end = A_ptr[k-m+1];

            for (ptr = beg; ptr < end; ptr++)

               temp += A_val[ptr] * x[A_ind[ptr]];

         }

         r[i] = temp;

      }

      return;

}

/***********************************************************************

*  refine_btran - refine solution of B'* x = h

*

*  This routine performs one iteration to refine the solution of the

*  system B'* x = h, where B' is a matrix transposed to the current

*  basis matrix, h is the vector of right-hand sides, x is the solution

*  vector. */

static void refine_btran(struct csa *csa, double h[], double x[])

{     int m = csa->m;

      double *r = csa->work1;

      double *d = csa->work1;

      int i;

      /* compute the residual vector r = h - B'* x */

      error_btran(csa, h, x, r);

      /* compute the correction vector d = inv(B') * r */

      xassert(csa->valid);

      bfd_btran(csa->bfd, d);

      /* refine the solution vector (new x) = (old x) + d */

      for (i = 1; i <= m; i++) x[i] += d[i];

      return;

}

/***********************************************************************

*  alloc_N - allocate matrix N

*

*  This routine determines maximal row lengths of matrix N, sets its

*  row pointers, and then allocates arrays N_ind and N_val.

*

*  Note that some fixed structural variables may temporarily become

*  double-bounded, so corresponding columns of matrix A should not be

*  ignored on calculating maximal row lengths of matrix N. */

static void alloc_N(struct csa *csa)

{     int m = csa->m;

      int n = csa->n;

      int *A_ptr = csa->A_ptr;

      int *A_ind = csa->A_ind;

      int *N_ptr = csa->N_ptr;

      int *N_len = csa->N_len;

      int i, j, beg, end, ptr;

      /* determine number of non-zeros in each row of the augmented

         constraint matrix (I|-A) */

      for (i = 1; i <= m; i++)

         N_len[i] = 1;

      for (j = 1; j <= n; j++)

      {  beg = A_ptr[j];

         end = A_ptr[j+1];

         for (ptr = beg; ptr < end; ptr++)

            N_len[A_ind[ptr]]++;

      }

      /* determine maximal row lengths of matrix N and set its row

         pointers */

      N_ptr[1] = 1;

      for (i = 1; i <= m; i++)

      {  /* row of matrix N cannot have more than n non-zeros */

         if (N_len[i] > n) N_len[i] = n;

         N_ptr[i+1] = N_ptr[i] + N_len[i];

      }

      /* now maximal number of non-zeros in matrix N is known */

      csa->N_ind = xcalloc(N_ptr[m+1], sizeof(int));

      csa->N_val = xcalloc(N_ptr[m+1], sizeof(double));

      return;

}

/***********************************************************************

*  add_N_col - add column of matrix (I|-A) to matrix N

*

*  This routine adds j-th column to matrix N which is k-th column of

*  the augmented constraint matrix (I|-A). (It is assumed that old j-th

*  column was previously removed from matrix N.) */

static void add_N_col(struct csa *csa, int j, int k)

{     int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      int *N_ptr = csa->N_ptr;

      int *N_len = csa->N_len;

      int *N_ind = csa->N_ind;

      double *N_val = csa->N_val;

      int pos;

#ifdef GLP_DEBUG

      xassert(1 <= j && j <= n);

      xassert(1 <= k && k <= m+n);

#endif

      if (k <= m)

      {  /* N[j] is k-th column of submatrix I */

         pos = N_ptr[k] + (N_len[k]++);

#ifdef GLP_DEBUG

         xassert(pos < N_ptr[k+1]);

#endif

         N_ind[pos] = j;

         N_val[pos] = 1.0;

      }

      else

      {  /* N[j] is (k-m)-th column of submatrix (-A) */

         int *A_ptr = csa->A_ptr;

         int *A_ind = csa->A_ind;

         double *A_val = csa->A_val;

         int i, beg, end, ptr;

         beg = A_ptr[k-m];

         end = A_ptr[k-m+1];

         for (ptr = beg; ptr < end; ptr++)

         {  i = A_ind[ptr]; /* row number */

            pos = N_ptr[i] + (N_len[i]++);

#ifdef GLP_DEBUG

            xassert(pos < N_ptr[i+1]);

#endif

            N_ind[pos] = j;

            N_val[pos] = - A_val[ptr];

         }

      }

      return;

}

/***********************************************************************

*  del_N_col - remove column of matrix (I|-A) from matrix N

*

*  This routine removes j-th column from matrix N which is k-th column

*  of the augmented constraint matrix (I|-A). */

static void del_N_col(struct csa *csa, int j, int k)

{     int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      int *N_ptr = csa->N_ptr;

      int *N_len = csa->N_len;

      int *N_ind = csa->N_ind;

      double *N_val = csa->N_val;

      int pos, head, tail;

#ifdef GLP_DEBUG

      xassert(1 <= j && j <= n);

      xassert(1 <= k && k <= m+n);

#endif

      if (k <= m)

      {  /* N[j] is k-th column of submatrix I */

         /* find element in k-th row of N */

         head = N_ptr[k];

         for (pos = head; N_ind[pos] != j; pos++) /* nop */;

         /* and remove it from the row list */

         tail = head + (--N_len[k]);

#ifdef GLP_DEBUG

         xassert(pos <= tail);

#endif

         N_ind[pos] = N_ind[tail];

         N_val[pos] = N_val[tail];

      }

      else

      {  /* N[j] is (k-m)-th column of submatrix (-A) */

         int *A_ptr = csa->A_ptr;

         int *A_ind = csa->A_ind;

         int i, beg, end, ptr;

         beg = A_ptr[k-m];

         end = A_ptr[k-m+1];

         for (ptr = beg; ptr < end; ptr++)

         {  i = A_ind[ptr]; /* row number */

            /* find element in i-th row of N */

            head = N_ptr[i];

            for (pos = head; N_ind[pos] != j; pos++) /* nop */;

            /* and remove it from the row list */

            tail = head + (--N_len[i]);

#ifdef GLP_DEBUG

            xassert(pos <= tail);

#endif

            N_ind[pos] = N_ind[tail];

            N_val[pos] = N_val[tail];

         }

      }

      return;

}

/***********************************************************************

*  build_N - build matrix N for current basis

*

*  This routine builds matrix N for the current basis from columns

*  of the augmented constraint matrix (I|-A) corresponding to non-basic

*  non-fixed variables. */

static void build_N(struct csa *csa)

{     int m = csa->m;

      int n = csa->n;

      int *head = csa->head;

      char *stat = csa->stat;

      int *N_len = csa->N_len;

      int j, k;

      /* N := empty matrix */

      memset(&N_len[1], 0, m * sizeof(int));

      /* go through non-basic columns of matrix (I|-A) */

      for (j = 1; j <= n; j++)

      {  if (stat[j] != GLP_NS)

         {  /* xN[j] is non-fixed; add j-th column to matrix N which is

               k-th column of matrix (I|-A) */

            k = head[m+j]; /* x[k] = xN[j] */

#ifdef GLP_DEBUG

            xassert(1 <= k && k <= m+n);

#endif

            add_N_col(csa, j, k);

         }

      }

      return;

}

/***********************************************************************

*  get_xN - determine current value of non-basic variable xN[j]

*

*  This routine returns the current value of non-basic variable xN[j],

*  which is a value of its active bound. */

static double get_xN(struct csa *csa, int j)

{     int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      double *lb = csa->lb;

      double *ub = csa->ub;

      int *head = csa->head;

      char *stat = csa->stat;

      int k;

      double xN;

#ifdef GLP_DEBUG

      xassert(1 <= j && j <= n);

#endif

      k = head[m+j]; /* x[k] = xN[j] */

#ifdef GLP_DEBUG

      xassert(1 <= k && k <= m+n);

#endif

      switch (stat[j])

      {  case GLP_NL:

            /* x[k] is on its lower bound */

            xN = lb[k]; break;

         case GLP_NU:

            /* x[k] is on its upper bound */

            xN = ub[k]; break;

         case GLP_NF:

            /* x[k] is free non-basic variable */

            xN = 0.0; break;

         case GLP_NS:

            /* x[k] is fixed non-basic variable */

            xN = lb[k]; break;

         default:

            xassert(stat != stat);

      }

      return xN;

}

/***********************************************************************

*  eval_beta - compute primal values of basic variables

*

*  This routine computes current primal values of all basic variables:

*

*     beta = - inv(B) * N * xN,

*

*  where B is the current basis matrix, N is a matrix built of columns

*  of matrix (I|-A) corresponding to non-basic variables, and xN is the

*  vector of current values of non-basic variables. */

static void eval_beta(struct csa *csa, double beta[])

{     int m = csa->m;

      int n = csa->n;

      int *A_ptr = csa->A_ptr;

      int *A_ind = csa->A_ind;

      double *A_val = csa->A_val;

      int *head = csa->head;

      double *h = csa->work2;

      int i, j, k, beg, end, ptr;

      double xN;

      /* compute the right-hand side vector:

         h := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n],

         where N[1], ..., N[n] are columns of matrix N */

      for (i = 1; i <= m; i++)

         h[i] = 0.0;

      for (j = 1; j <= n; j++)

      {  k = head[m+j]; /* x[k] = xN[j] */

#ifdef GLP_DEBUG

         xassert(1 <= k && k <= m+n);

#endif

         /* determine current value of xN[j] */

         xN = get_xN(csa, j);

         if (xN == 0.0) continue;

         if (k <= m)

         {  /* N[j] is k-th column of submatrix I */

            h[k] -= xN;

         }

         else

         {  /* N[j] is (k-m)-th column of submatrix (-A) */

            beg = A_ptr[k-m];

            end = A_ptr[k-m+1];

            for (ptr = beg; ptr < end; ptr++)

               h[A_ind[ptr]] += xN * A_val[ptr];

         }

      }

      /* solve system B * beta = h */

      memcpy(&beta[1], &h[1], m * sizeof(double));

      xassert(csa->valid);

      bfd_ftran(csa->bfd, beta);

      /* and refine the solution */

      refine_ftran(csa, h, beta);

      return;

}

/***********************************************************************

*  eval_pi - compute vector of simplex multipliers

*

*  This routine computes the vector of current simplex multipliers:

*

*     pi = inv(B') * cB,

*

*  where B' is a matrix transposed to the current basis matrix, cB is

*  a subvector of objective coefficients at basic variables. */

static void eval_pi(struct csa *csa, double pi[])

{     int m = csa->m;

      double *c = csa->coef;

      int *head = csa->head;

      double *cB = csa->work2;

      int i;

      /* construct the right-hand side vector cB */

      for (i = 1; i <= m; i++)

         cB[i] = c[head[i]];

      /* solve system B'* pi = cB */

      memcpy(&pi[1], &cB[1], m * sizeof(double));

      xassert(csa->valid);

      bfd_btran(csa->bfd, pi);

      /* and refine the solution */

      refine_btran(csa, cB, pi);

      return;

}

/***********************************************************************

*  eval_cost - compute reduced cost of non-basic variable xN[j]

*

*  This routine computes the current reduced cost of non-basic variable

*  xN[j]:

*

*     d[j] = cN[j] - N'[j] * pi,

*

*  where cN[j] is the objective coefficient at variable xN[j], N[j] is

*  a column of the augmented constraint matrix (I|-A) corresponding to

*  xN[j], pi is the vector of simplex multipliers. */

static double eval_cost(struct csa *csa, double pi[], int j)

{     int m = csa->m;

#ifdef GLP_DEBUG

      int n = csa->n;

#endif

      double *coef = csa->coef;

      int *head = csa->head;

      int k;

      double dj;

#ifdef GLP_DEBUG

      xassert(1 <= j && j <= n);

#endif

      k = head[m+j]; /* x[k] = xN[j] */

#ifdef GLP_DEBUG

      xassert(1 <= k && k <= m+n);

#endif

      dj = coef[k];

      if (k <= m)

      {  /* N[j] is k-th column of submatrix I */

         dj -= pi[k];

      }

      else

      {  /* N[j] is (k-m)-th column of submatrix (-A) */

         int *A_ptr = csa->A_ptr;

         int *A_ind = csa->A_ind;

         double *A_val = csa->A_val;

         int beg, end, ptr;

         beg = A_ptr[k-m];

         end = A_ptr[k-m+1];

         for (ptr = beg; ptr < end; ptr++)

            dj += A_val[ptr] * pi[A_ind[ptr]];

      }

      return dj;

}

/***********************************************************************

*  eval_bbar - compute and store primal values of basic variables

*

*  This routine computes primal values of all basic variables and then

*  stores them in the solution array. */