COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpspx01.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 96.6 KB
Line 
1/* glpspx01.c (primal simplex method) */
2
3/***********************************************************************
4*  This code is part of GLPK (GNU Linear Programming Kit).
5*
6*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7*  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
8*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9*  E-mail: <mao@gnu.org>.
10*
11*  GLPK is free software: you can redistribute it and/or modify it
12*  under the terms of the GNU General Public License as published by
13*  the Free Software Foundation, either version 3 of the License, or
14*  (at your option) any later version.
15*
16*  GLPK is distributed in the hope that it will be useful, but WITHOUT
17*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19*  License for more details.
20*
21*  You should have received a copy of the GNU General Public License
22*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23***********************************************************************/
24
25#include "glpspx.h"
26
27struct csa
28{     /* common storage area */
29      /*--------------------------------------------------------------*/
30      /* LP data */
31      int m;
32      /* number of rows (auxiliary variables), m > 0 */
33      int n;
34      /* number of columns (structural variables), n > 0 */
35      char *type; /* char type[1+m+n]; */
36      /* type[0] is not used;
37         type[k], 1 <= k <= m+n, is the type of variable x[k]:
38         GLP_FR - free variable
39         GLP_LO - variable with lower bound
40         GLP_UP - variable with upper bound
41         GLP_DB - double-bounded variable
42         GLP_FX - fixed variable */
43      double *lb; /* double lb[1+m+n]; */
44      /* lb[0] is not used;
45         lb[k], 1 <= k <= m+n, is an lower bound of variable x[k];
46         if x[k] has no lower bound, lb[k] is zero */
47      double *ub; /* double ub[1+m+n]; */
48      /* ub[0] is not used;
49         ub[k], 1 <= k <= m+n, is an upper bound of variable x[k];
50         if x[k] has no upper bound, ub[k] is zero;
51         if x[k] is of fixed type, ub[k] is the same as lb[k] */
52      double *coef; /* double coef[1+m+n]; */
53      /* coef[0] is not used;
54         coef[k], 1 <= k <= m+n, is an objective coefficient at
55         variable x[k] (note that on phase I auxiliary variables also
56         may have non-zero objective coefficients) */
57      /*--------------------------------------------------------------*/
58      /* original objective function */
59      double *obj; /* double obj[1+n]; */
60      /* obj[0] is a constant term of the original objective function;
61         obj[j], 1 <= j <= n, is an original objective coefficient at
62         structural variable x[m+j] */
63      double zeta;
64      /* factor used to scale original objective coefficients; its
65         sign defines original optimization direction: zeta > 0 means
66         minimization, zeta < 0 means maximization */
67      /*--------------------------------------------------------------*/
68      /* constraint matrix A; it has m rows and n columns and is stored
69         by columns */
70      int *A_ptr; /* int A_ptr[1+n+1]; */
71      /* A_ptr[0] is not used;
72         A_ptr[j], 1 <= j <= n, is starting position of j-th column in
73         arrays A_ind and A_val; note that A_ptr[1] is always 1;
74         A_ptr[n+1] indicates the position after the last element in
75         arrays A_ind and A_val */
76      int *A_ind; /* int A_ind[A_ptr[n+1]]; */
77      /* row indices */
78      double *A_val; /* double A_val[A_ptr[n+1]]; */
79      /* non-zero element values */
80      /*--------------------------------------------------------------*/
81      /* basis header */
82      int *head; /* int head[1+m+n]; */
83      /* head[0] is not used;
84         head[i], 1 <= i <= m, is the ordinal number of basic variable
85         xB[i]; head[i] = k means that xB[i] = x[k] and i-th column of
86         matrix B is k-th column of matrix (I|-A);
87         head[m+j], 1 <= j <= n, is the ordinal number of non-basic
88         variable xN[j]; head[m+j] = k means that xN[j] = x[k] and j-th
89         column of matrix N is k-th column of matrix (I|-A) */
90      char *stat; /* char stat[1+n]; */
91      /* stat[0] is not used;
92         stat[j], 1 <= j <= n, is the status of non-basic variable
93         xN[j], which defines its active bound:
94         GLP_NL - lower bound is active
95         GLP_NU - upper bound is active
96         GLP_NF - free variable
97         GLP_NS - fixed variable */
98      /*--------------------------------------------------------------*/
99      /* matrix B is the basis matrix; it is composed from columns of
100         the augmented constraint matrix (I|-A) corresponding to basic
101         variables and stored in a factorized (invertable) form */
102      int valid;
103      /* factorization is valid only if this flag is set */
104      BFD *bfd; /* BFD bfd[1:m,1:m]; */
105      /* factorized (invertable) form of the basis matrix */
106      /*--------------------------------------------------------------*/
107      /* matrix N is a matrix composed from columns of the augmented
108         constraint matrix (I|-A) corresponding to non-basic variables
109         except fixed ones; it is stored by rows and changes every time
110         the basis changes */
111      int *N_ptr; /* int N_ptr[1+m+1]; */
112      /* N_ptr[0] is not used;
113         N_ptr[i], 1 <= i <= m, is starting position of i-th row in
114         arrays N_ind and N_val; note that N_ptr[1] is always 1;
115         N_ptr[m+1] indicates the position after the last element in
116         arrays N_ind and N_val */
117      int *N_len; /* int N_len[1+m]; */
118      /* N_len[0] is not used;
119         N_len[i], 1 <= i <= m, is length of i-th row (0 to n) */
120      int *N_ind; /* int N_ind[N_ptr[m+1]]; */
121      /* column indices */
122      double *N_val; /* double N_val[N_ptr[m+1]]; */
123      /* non-zero element values */
124      /*--------------------------------------------------------------*/
125      /* working parameters */
126      int phase;
127      /* search phase:
128         0 - not determined yet
129         1 - search for primal feasible solution
130         2 - search for optimal solution */
131      glp_long tm_beg;
132      /* time value at the beginning of the search */
133      int it_beg;
134      /* simplex iteration count at the beginning of the search */
135      int it_cnt;
136      /* simplex iteration count; it increases by one every time the
137         basis changes (including the case when a non-basic variable
138         jumps to its opposite bound) */
139      int it_dpy;
140      /* simplex iteration count at the most recent display output */
141      /*--------------------------------------------------------------*/
142      /* basic solution components */
143      double *bbar; /* double bbar[1+m]; */
144      /* bbar[0] is not used;
145         bbar[i], 1 <= i <= m, is primal value of basic variable xB[i]
146         (if xB[i] is free, its primal value is not updated) */
147      double *cbar; /* double cbar[1+n]; */
148      /* cbar[0] is not used;
149         cbar[j], 1 <= j <= n, is reduced cost of non-basic variable
150         xN[j] (if xN[j] is fixed, its reduced cost is not updated) */
151      /*--------------------------------------------------------------*/
152      /* the following pricing technique options may be used:
153         GLP_PT_STD - standard ("textbook") pricing;
154         GLP_PT_PSE - projected steepest edge;
155         GLP_PT_DVX - Devex pricing (not implemented yet);
156         in case of GLP_PT_STD the reference space is not used, and all
157         steepest edge coefficients are set to 1 */
158      int refct;
159      /* this count is set to an initial value when the reference space
160         is defined and decreases by one every time the basis changes;
161         once this count reaches zero, the reference space is redefined
162         again */
163      char *refsp; /* char refsp[1+m+n]; */
164      /* refsp[0] is not used;
165         refsp[k], 1 <= k <= m+n, is the flag which means that variable
166         x[k] belongs to the current reference space */
167      double *gamma; /* double gamma[1+n]; */
168      /* gamma[0] is not used;
169         gamma[j], 1 <= j <= n, is the steepest edge coefficient for
170         non-basic variable xN[j]; if xN[j] is fixed, gamma[j] is not
171         used and just set to 1 */
172      /*--------------------------------------------------------------*/
173      /* non-basic variable xN[q] chosen to enter the basis */
174      int q;
175      /* index of the non-basic variable xN[q] chosen, 1 <= q <= n;
176         if the set of eligible non-basic variables is empty and thus
177         no variable has been chosen, q is set to 0 */
178      /*--------------------------------------------------------------*/
179      /* pivot column of the simplex table corresponding to non-basic
180         variable xN[q] chosen is the following vector:
181            T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
182         where B is the current basis matrix, N[q] is a column of the
183         matrix (I|-A) corresponding to xN[q] */
184      int tcol_nnz;
185      /* number of non-zero components, 0 <= nnz <= m */
186      int *tcol_ind; /* int tcol_ind[1+m]; */
187      /* tcol_ind[0] is not used;
188         tcol_ind[t], 1 <= t <= nnz, is an index of non-zero component,
189         i.e. tcol_ind[t] = i means that tcol_vec[i] != 0 */
190      double *tcol_vec; /* double tcol_vec[1+m]; */
191      /* tcol_vec[0] is not used;
192         tcol_vec[i], 1 <= i <= m, is a numeric value of i-th component
193         of the column */
194      double tcol_max;
195      /* infinity (maximum) norm of the column (max |tcol_vec[i]|) */
196      int tcol_num;
197      /* number of significant non-zero components, which means that:
198         |tcol_vec[i]| >= eps for i in tcol_ind[1,...,num],
199         |tcol_vec[i]| <  eps for i in tcol_ind[num+1,...,nnz],
200         where eps is a pivot tolerance */
201      /*--------------------------------------------------------------*/
202      /* basic variable xB[p] chosen to leave the basis */
203      int p;
204      /* index of the basic variable xB[p] chosen, 1 <= p <= m;
205         p = 0 means that no basic variable reaches its bound;
206         p < 0 means that non-basic variable xN[q] reaches its opposite
207         bound before any basic variable */
208      int p_stat;
209      /* new status (GLP_NL, GLP_NU, or GLP_NS) to be assigned to xB[p]
210         once it has left the basis */
211      double teta;
212      /* change of non-basic variable xN[q] (see above), on which xB[p]
213         (or, if p < 0, xN[q] itself) reaches its bound */
214      /*--------------------------------------------------------------*/
215      /* pivot row of the simplex table corresponding to basic variable
216         xB[p] chosen is the following vector:
217            T' * e[p] = - N' * inv(B') * e[p] = - N' * rho,
218         where B' is a matrix transposed to the current basis matrix,
219         N' is a matrix, whose rows are columns of the matrix (I|-A)
220         corresponding to non-basic non-fixed variables */
221      int trow_nnz;
222      /* number of non-zero components, 0 <= nnz <= n */
223      int *trow_ind; /* int trow_ind[1+n]; */
224      /* trow_ind[0] is not used;
225         trow_ind[t], 1 <= t <= nnz, is an index of non-zero component,
226         i.e. trow_ind[t] = j means that trow_vec[j] != 0 */
227      double *trow_vec; /* int trow_vec[1+n]; */
228      /* trow_vec[0] is not used;
229         trow_vec[j], 1 <= j <= n, is a numeric value of j-th component
230         of the row */
231      /*--------------------------------------------------------------*/
232      /* working arrays */
233      double *work1; /* double work1[1+m]; */
234      double *work2; /* double work2[1+m]; */
235      double *work3; /* double work3[1+m]; */
236      double *work4; /* double work4[1+m]; */
237};
238
239static const double kappa = 0.10;
240
241/***********************************************************************
242*  alloc_csa - allocate common storage area
243*
244*  This routine allocates all arrays in the common storage area (CSA)
245*  and returns a pointer to the CSA. */
246
247static struct csa *alloc_csa(glp_prob *lp)
248{     struct csa *csa;
249      int m = lp->m;
250      int n = lp->n;
251      int nnz = lp->nnz;
252      csa = xmalloc(sizeof(struct csa));
253      xassert(m > 0 && n > 0);
254      csa->m = m;
255      csa->n = n;
256      csa->type = xcalloc(1+m+n, sizeof(char));
257      csa->lb = xcalloc(1+m+n, sizeof(double));
258      csa->ub = xcalloc(1+m+n, sizeof(double));
259      csa->coef = xcalloc(1+m+n, sizeof(double));
260      csa->obj = xcalloc(1+n, sizeof(double));
261      csa->A_ptr = xcalloc(1+n+1, sizeof(int));
262      csa->A_ind = xcalloc(1+nnz, sizeof(int));
263      csa->A_val = xcalloc(1+nnz, sizeof(double));
264      csa->head = xcalloc(1+m+n, sizeof(int));
265      csa->stat = xcalloc(1+n, sizeof(char));
266      csa->N_ptr = xcalloc(1+m+1, sizeof(int));
267      csa->N_len = xcalloc(1+m, sizeof(int));
268      csa->N_ind = NULL; /* will be allocated later */
269      csa->N_val = NULL; /* will be allocated later */
270      csa->bbar = xcalloc(1+m, sizeof(double));
271      csa->cbar = xcalloc(1+n, sizeof(double));
272      csa->refsp = xcalloc(1+m+n, sizeof(char));
273      csa->gamma = xcalloc(1+n, sizeof(double));
274      csa->tcol_ind = xcalloc(1+m, sizeof(int));
275      csa->tcol_vec = xcalloc(1+m, sizeof(double));
276      csa->trow_ind = xcalloc(1+n, sizeof(int));
277      csa->trow_vec = xcalloc(1+n, sizeof(double));
278      csa->work1 = xcalloc(1+m, sizeof(double));
279      csa->work2 = xcalloc(1+m, sizeof(double));
280      csa->work3 = xcalloc(1+m, sizeof(double));
281      csa->work4 = xcalloc(1+m, sizeof(double));
282      return csa;
283}
284
285/***********************************************************************
286*  init_csa - initialize common storage area
287*
288*  This routine initializes all data structures in the common storage
289*  area (CSA). */
290
291static void alloc_N(struct csa *csa);
292static void build_N(struct csa *csa);
293
294static void init_csa(struct csa *csa, glp_prob *lp)
295{     int m = csa->m;
296      int n = csa->n;
297      char *type = csa->type;
298      double *lb = csa->lb;
299      double *ub = csa->ub;
300      double *coef = csa->coef;
301      double *obj = csa->obj;
302      int *A_ptr = csa->A_ptr;
303      int *A_ind = csa->A_ind;
304      double *A_val = csa->A_val;
305      int *head = csa->head;
306      char *stat = csa->stat;
307      char *refsp = csa->refsp;
308      double *gamma = csa->gamma;
309      int i, j, k, loc;
310      double cmax;
311      /* auxiliary variables */
312      for (i = 1; i <= m; i++)
313      {  GLPROW *row = lp->row[i];
314         type[i] = (char)row->type;
315         lb[i] = row->lb * row->rii;
316         ub[i] = row->ub * row->rii;
317         coef[i] = 0.0;
318      }
319      /* structural variables */
320      for (j = 1; j <= n; j++)
321      {  GLPCOL *col = lp->col[j];
322         type[m+j] = (char)col->type;
323         lb[m+j] = col->lb / col->sjj;
324         ub[m+j] = col->ub / col->sjj;
325         coef[m+j] = col->coef * col->sjj;
326      }
327      /* original objective function */
328      obj[0] = lp->c0;
329      memcpy(&obj[1], &coef[m+1], n * sizeof(double));
330      /* factor used to scale original objective coefficients */
331      cmax = 0.0;
332      for (j = 1; j <= n; j++)
333         if (cmax < fabs(obj[j])) cmax = fabs(obj[j]);
334      if (cmax == 0.0) cmax = 1.0;
335      switch (lp->dir)
336      {  case GLP_MIN:
337            csa->zeta = + 1.0 / cmax;
338            break;
339         case GLP_MAX:
340            csa->zeta = - 1.0 / cmax;
341            break;
342         default:
343            xassert(lp != lp);
344      }
345#if 1
346      if (fabs(csa->zeta) < 1.0) csa->zeta *= 1000.0;
347#endif
348      /* matrix A (by columns) */
349      loc = 1;
350      for (j = 1; j <= n; j++)
351      {  GLPAIJ *aij;
352         A_ptr[j] = loc;
353         for (aij = lp->col[j]->ptr; aij != NULL; aij = aij->c_next)
354         {  A_ind[loc] = aij->row->i;
355            A_val[loc] = aij->row->rii * aij->val * aij->col->sjj;
356            loc++;
357         }
358      }
359      A_ptr[n+1] = loc;
360      xassert(loc == lp->nnz+1);
361      /* basis header */
362      xassert(lp->valid);
363      memcpy(&head[1], &lp->head[1], m * sizeof(int));
364      k = 0;
365      for (i = 1; i <= m; i++)
366      {  GLPROW *row = lp->row[i];
367         if (row->stat != GLP_BS)
368         {  k++;
369            xassert(k <= n);
370            head[m+k] = i;
371            stat[k] = (char)row->stat;
372         }
373      }
374      for (j = 1; j <= n; j++)
375      {  GLPCOL *col = lp->col[j];
376         if (col->stat != GLP_BS)
377         {  k++;
378            xassert(k <= n);
379            head[m+k] = m + j;
380            stat[k] = (char)col->stat;
381         }
382      }
383      xassert(k == n);
384      /* factorization of matrix B */
385      csa->valid = 1, lp->valid = 0;
386      csa->bfd = lp->bfd, lp->bfd = NULL;
387      /* matrix N (by rows) */
388      alloc_N(csa);
389      build_N(csa);
390      /* working parameters */
391      csa->phase = 0;
392      csa->tm_beg = xtime();
393      csa->it_beg = csa->it_cnt = lp->it_cnt;
394      csa->it_dpy = -1;
395      /* reference space and steepest edge coefficients */
396      csa->refct = 0;
397      memset(&refsp[1], 0, (m+n) * sizeof(char));
398      for (j = 1; j <= n; j++) gamma[j] = 1.0;
399      return;
400}
401
402/***********************************************************************
403*  invert_B - compute factorization of the basis matrix
404*
405*  This routine computes factorization of the current basis matrix B.
406*
407*  If the operation is successful, the routine returns zero, otherwise
408*  non-zero. */
409
410static int inv_col(void *info, int i, int ind[], double val[])
411{     /* this auxiliary routine returns row indices and numeric values
412         of non-zero elements of i-th column of the basis matrix */
413      struct csa *csa = info;
414      int m = csa->m;
415#ifdef GLP_DEBUG
416      int n = csa->n;
417#endif
418      int *A_ptr = csa->A_ptr;
419      int *A_ind = csa->A_ind;
420      double *A_val = csa->A_val;
421      int *head = csa->head;
422      int k, len, ptr, t;
423#ifdef GLP_DEBUG
424      xassert(1 <= i && i <= m);
425#endif
426      k = head[i]; /* B[i] is k-th column of (I|-A) */
427#ifdef GLP_DEBUG
428      xassert(1 <= k && k <= m+n);
429#endif
430      if (k <= m)
431      {  /* B[i] is k-th column of submatrix I */
432         len = 1;
433         ind[1] = k;
434         val[1] = 1.0;
435      }
436      else
437      {  /* B[i] is (k-m)-th column of submatrix (-A) */
438         ptr = A_ptr[k-m];
439         len = A_ptr[k-m+1] - ptr;
440         memcpy(&ind[1], &A_ind[ptr], len * sizeof(int));
441         memcpy(&val[1], &A_val[ptr], len * sizeof(double));
442         for (t = 1; t <= len; t++) val[t] = - val[t];
443      }
444      return len;
445}
446
447static int invert_B(struct csa *csa)
448{     int ret;
449      ret = bfd_factorize(csa->bfd, csa->m, NULL, inv_col, csa);
450      csa->valid = (ret == 0);
451      return ret;
452}
453
454/***********************************************************************
455*  update_B - update factorization of the basis matrix
456*
457*  This routine replaces i-th column of the basis matrix B by k-th
458*  column of the augmented constraint matrix (I|-A) and then updates
459*  the factorization of B.
460*
461*  If the factorization has been successfully updated, the routine
462*  returns zero, otherwise non-zero. */
463
464static int update_B(struct csa *csa, int i, int k)
465{     int m = csa->m;
466#ifdef GLP_DEBUG
467      int n = csa->n;
468#endif
469      int ret;
470#ifdef GLP_DEBUG
471      xassert(1 <= i && i <= m);
472      xassert(1 <= k && k <= m+n);
473#endif
474      if (k <= m)
475      {  /* new i-th column of B is k-th column of I */
476         int ind[1+1];
477         double val[1+1];
478         ind[1] = k;
479         val[1] = 1.0;
480         xassert(csa->valid);
481         ret = bfd_update_it(csa->bfd, i, 0, 1, ind, val);
482      }
483      else
484      {  /* new i-th column of B is (k-m)-th column of (-A) */
485         int *A_ptr = csa->A_ptr;
486         int *A_ind = csa->A_ind;
487         double *A_val = csa->A_val;
488         double *val = csa->work1;
489         int beg, end, ptr, len;
490         beg = A_ptr[k-m];
491         end = A_ptr[k-m+1];
492         len = 0;
493         for (ptr = beg; ptr < end; ptr++)
494            val[++len] = - A_val[ptr];
495         xassert(csa->valid);
496         ret = bfd_update_it(csa->bfd, i, 0, len, &A_ind[beg-1], val);
497      }
498      csa->valid = (ret == 0);
499      return ret;
500}
501
502/***********************************************************************
503*  error_ftran - compute residual vector r = h - B * x
504*
505*  This routine computes the residual vector r = h - B * x, where B is
506*  the current basis matrix, h is the vector of right-hand sides, x is
507*  the solution vector. */
508
509static void error_ftran(struct csa *csa, double h[], double x[],
510      double r[])
511{     int m = csa->m;
512#ifdef GLP_DEBUG
513      int n = csa->n;
514#endif
515      int *A_ptr = csa->A_ptr;
516      int *A_ind = csa->A_ind;
517      double *A_val = csa->A_val;
518      int *head = csa->head;
519      int i, k, beg, end, ptr;
520      double temp;
521      /* compute the residual vector:
522         r = h - B * x = h - B[1] * x[1] - ... - B[m] * x[m],
523         where B[1], ..., B[m] are columns of matrix B */
524      memcpy(&r[1], &h[1], m * sizeof(double));
525      for (i = 1; i <= m; i++)
526      {  temp = x[i];
527         if (temp == 0.0) continue;
528         k = head[i]; /* B[i] is k-th column of (I|-A) */
529#ifdef GLP_DEBUG
530         xassert(1 <= k && k <= m+n);
531#endif
532         if (k <= m)
533         {  /* B[i] is k-th column of submatrix I */
534            r[k] -= temp;
535         }
536         else
537         {  /* B[i] is (k-m)-th column of submatrix (-A) */
538            beg = A_ptr[k-m];
539            end = A_ptr[k-m+1];
540            for (ptr = beg; ptr < end; ptr++)
541               r[A_ind[ptr]] += A_val[ptr] * temp;
542         }
543      }
544      return;
545}
546
547/***********************************************************************
548*  refine_ftran - refine solution of B * x = h
549*
550*  This routine performs one iteration to refine the solution of
551*  the system B * x = h, where B is the current basis matrix, h is the
552*  vector of right-hand sides, x is the solution vector. */
553
554static void refine_ftran(struct csa *csa, double h[], double x[])
555{     int m = csa->m;
556      double *r = csa->work1;
557      double *d = csa->work1;
558      int i;
559      /* compute the residual vector r = h - B * x */
560      error_ftran(csa, h, x, r);
561      /* compute the correction vector d = inv(B) * r */
562      xassert(csa->valid);
563      bfd_ftran(csa->bfd, d);
564      /* refine the solution vector (new x) = (old x) + d */
565      for (i = 1; i <= m; i++) x[i] += d[i];
566      return;
567}
568
569/***********************************************************************
570*  error_btran - compute residual vector r = h - B'* x
571*
572*  This routine computes the residual vector r = h - B'* x, where B'
573*  is a matrix transposed to the current basis matrix, h is the vector
574*  of right-hand sides, x is the solution vector. */
575
576static void error_btran(struct csa *csa, double h[], double x[],
577      double r[])
578{     int m = csa->m;
579#ifdef GLP_DEBUG
580      int n = csa->n;
581#endif
582      int *A_ptr = csa->A_ptr;
583      int *A_ind = csa->A_ind;
584      double *A_val = csa->A_val;
585      int *head = csa->head;
586      int i, k, beg, end, ptr;
587      double temp;
588      /* compute the residual vector r = b - B'* x */
589      for (i = 1; i <= m; i++)
590      {  /* r[i] := b[i] - (i-th column of B)'* x */
591         k = head[i]; /* B[i] is k-th column of (I|-A) */
592#ifdef GLP_DEBUG
593         xassert(1 <= k && k <= m+n);
594#endif
595         temp = h[i];
596         if (k <= m)
597         {  /* B[i] is k-th column of submatrix I */
598            temp -= x[k];
599         }
600         else
601         {  /* B[i] is (k-m)-th column of submatrix (-A) */
602            beg = A_ptr[k-m];
603            end = A_ptr[k-m+1];
604            for (ptr = beg; ptr < end; ptr++)
605               temp += A_val[ptr] * x[A_ind[ptr]];
606         }
607         r[i] = temp;
608      }
609      return;
610}
611
612/***********************************************************************
613*  refine_btran - refine solution of B'* x = h
614*
615*  This routine performs one iteration to refine the solution of the
616*  system B'* x = h, where B' is a matrix transposed to the current
617*  basis matrix, h is the vector of right-hand sides, x is the solution
618*  vector. */
619
620static void refine_btran(struct csa *csa, double h[], double x[])
621{     int m = csa->m;
622      double *r = csa->work1;
623      double *d = csa->work1;
624      int i;
625      /* compute the residual vector r = h - B'* x */
626      error_btran(csa, h, x, r);
627      /* compute the correction vector d = inv(B') * r */
628      xassert(csa->valid);
629      bfd_btran(csa->bfd, d);
630      /* refine the solution vector (new x) = (old x) + d */
631      for (i = 1; i <= m; i++) x[i] += d[i];
632      return;
633}
634
635/***********************************************************************
636*  alloc_N - allocate matrix N
637*
638*  This routine determines maximal row lengths of matrix N, sets its
639*  row pointers, and then allocates arrays N_ind and N_val.
640*
641*  Note that some fixed structural variables may temporarily become
642*  double-bounded, so corresponding columns of matrix A should not be
643*  ignored on calculating maximal row lengths of matrix N. */
644
645static void alloc_N(struct csa *csa)
646{     int m = csa->m;
647      int n = csa->n;
648      int *A_ptr = csa->A_ptr;
649      int *A_ind = csa->A_ind;
650      int *N_ptr = csa->N_ptr;
651      int *N_len = csa->N_len;
652      int i, j, beg, end, ptr;
653      /* determine number of non-zeros in each row of the augmented
654         constraint matrix (I|-A) */
655      for (i = 1; i <= m; i++)
656         N_len[i] = 1;
657      for (j = 1; j <= n; j++)
658      {  beg = A_ptr[j];
659         end = A_ptr[j+1];
660         for (ptr = beg; ptr < end; ptr++)
661            N_len[A_ind[ptr]]++;
662      }
663      /* determine maximal row lengths of matrix N and set its row
664         pointers */
665      N_ptr[1] = 1;
666      for (i = 1; i <= m; i++)
667      {  /* row of matrix N cannot have more than n non-zeros */
668         if (N_len[i] > n) N_len[i] = n;
669         N_ptr[i+1] = N_ptr[i] + N_len[i];
670      }
671      /* now maximal number of non-zeros in matrix N is known */
672      csa->N_ind = xcalloc(N_ptr[m+1], sizeof(int));
673      csa->N_val = xcalloc(N_ptr[m+1], sizeof(double));
674      return;
675}
676
677/***********************************************************************
678*  add_N_col - add column of matrix (I|-A) to matrix N
679*
680*  This routine adds j-th column to matrix N which is k-th column of
681*  the augmented constraint matrix (I|-A). (It is assumed that old j-th
682*  column was previously removed from matrix N.) */
683
684static void add_N_col(struct csa *csa, int j, int k)
685{     int m = csa->m;
686#ifdef GLP_DEBUG
687      int n = csa->n;
688#endif
689      int *N_ptr = csa->N_ptr;
690      int *N_len = csa->N_len;
691      int *N_ind = csa->N_ind;
692      double *N_val = csa->N_val;
693      int pos;
694#ifdef GLP_DEBUG
695      xassert(1 <= j && j <= n);
696      xassert(1 <= k && k <= m+n);
697#endif
698      if (k <= m)
699      {  /* N[j] is k-th column of submatrix I */
700         pos = N_ptr[k] + (N_len[k]++);
701#ifdef GLP_DEBUG
702         xassert(pos < N_ptr[k+1]);
703#endif
704         N_ind[pos] = j;
705         N_val[pos] = 1.0;
706      }
707      else
708      {  /* N[j] is (k-m)-th column of submatrix (-A) */
709         int *A_ptr = csa->A_ptr;
710         int *A_ind = csa->A_ind;
711         double *A_val = csa->A_val;
712         int i, beg, end, ptr;
713         beg = A_ptr[k-m];
714         end = A_ptr[k-m+1];
715         for (ptr = beg; ptr < end; ptr++)
716         {  i = A_ind[ptr]; /* row number */
717            pos = N_ptr[i] + (N_len[i]++);
718#ifdef GLP_DEBUG
719            xassert(pos < N_ptr[i+1]);
720#endif
721            N_ind[pos] = j;
722            N_val[pos] = - A_val[ptr];
723         }
724      }
725      return;
726}
727
728/***********************************************************************
729*  del_N_col - remove column of matrix (I|-A) from matrix N
730*
731*  This routine removes j-th column from matrix N which is k-th column
732*  of the augmented constraint matrix (I|-A). */
733
734static void del_N_col(struct csa *csa, int j, int k)
735{     int m = csa->m;
736#ifdef GLP_DEBUG
737      int n = csa->n;
738#endif
739      int *N_ptr = csa->N_ptr;
740      int *N_len = csa->N_len;
741      int *N_ind = csa->N_ind;
742      double *N_val = csa->N_val;
743      int pos, head, tail;
744#ifdef GLP_DEBUG
745      xassert(1 <= j && j <= n);
746      xassert(1 <= k && k <= m+n);
747#endif
748      if (k <= m)
749      {  /* N[j] is k-th column of submatrix I */
750         /* find element in k-th row of N */
751         head = N_ptr[k];
752         for (pos = head; N_ind[pos] != j; pos++) /* nop */;
753         /* and remove it from the row list */
754         tail = head + (--N_len[k]);
755#ifdef GLP_DEBUG
756         xassert(pos <= tail);
757#endif
758         N_ind[pos] = N_ind[tail];
759         N_val[pos] = N_val[tail];
760      }
761      else
762      {  /* N[j] is (k-m)-th column of submatrix (-A) */
763         int *A_ptr = csa->A_ptr;
764         int *A_ind = csa->A_ind;
765         int i, beg, end, ptr;
766         beg = A_ptr[k-m];
767         end = A_ptr[k-m+1];
768         for (ptr = beg; ptr < end; ptr++)
769         {  i = A_ind[ptr]; /* row number */
770            /* find element in i-th row of N */
771            head = N_ptr[i];
772            for (pos = head; N_ind[pos] != j; pos++) /* nop */;
773            /* and remove it from the row list */
774            tail = head + (--N_len[i]);
775#ifdef GLP_DEBUG
776            xassert(pos <= tail);
777#endif
778            N_ind[pos] = N_ind[tail];
779            N_val[pos] = N_val[tail];
780         }
781      }
782      return;
783}
784
785/***********************************************************************
786*  build_N - build matrix N for current basis
787*
788*  This routine builds matrix N for the current basis from columns
789*  of the augmented constraint matrix (I|-A) corresponding to non-basic
790*  non-fixed variables. */
791
792static void build_N(struct csa *csa)
793{     int m = csa->m;
794      int n = csa->n;
795      int *head = csa->head;
796      char *stat = csa->stat;
797      int *N_len = csa->N_len;
798      int j, k;
799      /* N := empty matrix */
800      memset(&N_len[1], 0, m * sizeof(int));
801      /* go through non-basic columns of matrix (I|-A) */
802      for (j = 1; j <= n; j++)
803      {  if (stat[j] != GLP_NS)
804         {  /* xN[j] is non-fixed; add j-th column to matrix N which is
805               k-th column of matrix (I|-A) */
806            k = head[m+j]; /* x[k] = xN[j] */
807#ifdef GLP_DEBUG
808            xassert(1 <= k && k <= m+n);
809#endif
810            add_N_col(csa, j, k);
811         }
812      }
813      return;
814}
815
816/***********************************************************************
817*  get_xN - determine current value of non-basic variable xN[j]
818*
819*  This routine returns the current value of non-basic variable xN[j],
820*  which is a value of its active bound. */
821
822static double get_xN(struct csa *csa, int j)
823{     int m = csa->m;
824#ifdef GLP_DEBUG
825      int n = csa->n;
826#endif
827      double *lb = csa->lb;
828      double *ub = csa->ub;
829      int *head = csa->head;
830      char *stat = csa->stat;
831      int k;
832      double xN;
833#ifdef GLP_DEBUG
834      xassert(1 <= j && j <= n);
835#endif
836      k = head[m+j]; /* x[k] = xN[j] */
837#ifdef GLP_DEBUG
838      xassert(1 <= k && k <= m+n);
839#endif
840      switch (stat[j])
841      {  case GLP_NL:
842            /* x[k] is on its lower bound */
843            xN = lb[k]; break;
844         case GLP_NU:
845            /* x[k] is on its upper bound */
846            xN = ub[k]; break;
847         case GLP_NF:
848            /* x[k] is free non-basic variable */
849            xN = 0.0; break;
850         case GLP_NS:
851            /* x[k] is fixed non-basic variable */
852            xN = lb[k]; break;
853         default:
854            xassert(stat != stat);
855      }
856      return xN;
857}
858
859/***********************************************************************
860*  eval_beta - compute primal values of basic variables
861*
862*  This routine computes current primal values of all basic variables:
863*
864*     beta = - inv(B) * N * xN,
865*
866*  where B is the current basis matrix, N is a matrix built of columns
867*  of matrix (I|-A) corresponding to non-basic variables, and xN is the
868*  vector of current values of non-basic variables. */
869
870static void eval_beta(struct csa *csa, double beta[])
871{     int m = csa->m;
872      int n = csa->n;
873      int *A_ptr = csa->A_ptr;
874      int *A_ind = csa->A_ind;
875      double *A_val = csa->A_val;
876      int *head = csa->head;
877      double *h = csa->work2;
878      int i, j, k, beg, end, ptr;
879      double xN;
880      /* compute the right-hand side vector:
881         h := - N * xN = - N[1] * xN[1] - ... - N[n] * xN[n],
882         where N[1], ..., N[n] are columns of matrix N */
883      for (i = 1; i <= m; i++)
884         h[i] = 0.0;
885      for (j = 1; j <= n; j++)
886      {  k = head[m+j]; /* x[k] = xN[j] */
887#ifdef GLP_DEBUG
888         xassert(1 <= k && k <= m+n);
889#endif
890         /* determine current value of xN[j] */
891         xN = get_xN(csa, j);
892         if (xN == 0.0) continue;
893         if (k <= m)
894         {  /* N[j] is k-th column of submatrix I */
895            h[k] -= xN;
896         }
897         else
898         {  /* N[j] is (k-m)-th column of submatrix (-A) */
899            beg = A_ptr[k-m];
900            end = A_ptr[k-m+1];
901            for (ptr = beg; ptr < end; ptr++)
902               h[A_ind[ptr]] += xN * A_val[ptr];
903         }
904      }
905      /* solve system B * beta = h */
906      memcpy(&beta[1], &h[1], m * sizeof(double));
907      xassert(csa->valid);
908      bfd_ftran(csa->bfd, beta);
909      /* and refine the solution */
910      refine_ftran(csa, h, beta);
911      return;
912}
913
914/***********************************************************************
915*  eval_pi - compute vector of simplex multipliers
916*
917*  This routine computes the vector of current simplex multipliers:
918*
919*     pi = inv(B') * cB,
920*
921*  where B' is a matrix transposed to the current basis matrix, cB is
922*  a subvector of objective coefficients at basic variables. */
923
924static void eval_pi(struct csa *csa, double pi[])
925{     int m = csa->m;
926      double *c = csa->coef;
927      int *head = csa->head;
928      double *cB = csa->work2;
929      int i;
930      /* construct the right-hand side vector cB */
931      for (i = 1; i <= m; i++)
932         cB[i] = c[head[i]];
933      /* solve system B'* pi = cB */
934      memcpy(&pi[1], &cB[1], m * sizeof(double));
935      xassert(csa->valid);
936      bfd_btran(csa->bfd, pi);
937      /* and refine the solution */
938      refine_btran(csa, cB, pi);
939      return;
940}
941
942/***********************************************************************
943*  eval_cost - compute reduced cost of non-basic variable xN[j]
944*
945*  This routine computes the current reduced cost of non-basic variable
946*  xN[j]:
947*
948*     d[j] = cN[j] - N'[j] * pi,
949*
950*  where cN[j] is the objective coefficient at variable xN[j], N[j] is
951*  a column of the augmented constraint matrix (I|-A) corresponding to
952*  xN[j], pi is the vector of simplex multipliers. */
953
954static double eval_cost(struct csa *csa, double pi[], int j)
955{     int m = csa->m;
956#ifdef GLP_DEBUG
957      int n = csa->n;
958#endif
959      double *coef = csa->coef;
960      int *head = csa->head;
961      int k;
962      double dj;
963#ifdef GLP_DEBUG
964      xassert(1 <= j && j <= n);
965#endif
966      k = head[m+j]; /* x[k] = xN[j] */
967#ifdef GLP_DEBUG
968      xassert(1 <= k && k <= m+n);
969#endif
970      dj = coef[k];
971      if (k <= m)
972      {  /* N[j] is k-th column of submatrix I */
973         dj -= pi[k];
974      }
975      else
976      {  /* N[j] is (k-m)-th column of submatrix (-A) */
977         int *A_ptr = csa->A_ptr;
978         int *A_ind = csa->A_ind;
979         double *A_val = csa->A_val;
980         int beg, end, ptr;
981         beg = A_ptr[k-m];
982         end = A_ptr[k-m+1];
983         for (ptr = beg; ptr < end; ptr++)
984            dj += A_val[ptr] * pi[A_ind[ptr]];
985      }
986      return dj;
987}
988
989/***********************************************************************
990*  eval_bbar - compute and store primal values of basic variables
991*
992*  This routine computes primal values of all basic variables and then
993*  stores them in the solution array. */
994
995static void eval_bbar(struct csa *csa)
996{     eval_beta(csa, csa->bbar);
997      return;
998}
999
1000/***********************************************************************
1001*  eval_cbar - compute and store reduced costs of non-basic variables
1002*
1003*  This routine computes reduced costs of all non-basic variables and
1004*  then stores them in the solution array. */
1005
1006static void eval_cbar(struct csa *csa)
1007{
1008#ifdef GLP_DEBUG
1009      int m = csa->m;
1010#endif
1011      int n = csa->n;
1012#ifdef GLP_DEBUG
1013      int *head = csa->head;
1014#endif
1015      double *cbar = csa->cbar;
1016      double *pi = csa->work3;
1017      int j;
1018#ifdef GLP_DEBUG
1019      int k;
1020#endif
1021      /* compute simplex multipliers */
1022      eval_pi(csa, pi);
1023      /* compute and store reduced costs */
1024      for (j = 1; j <= n; j++)
1025      {
1026#ifdef GLP_DEBUG
1027         k = head[m+j]; /* x[k] = xN[j] */
1028         xassert(1 <= k && k <= m+n);
1029#endif
1030         cbar[j] = eval_cost(csa, pi, j);
1031      }
1032      return;
1033}
1034
1035/***********************************************************************
1036*  reset_refsp - reset the reference space
1037*
1038*  This routine resets (redefines) the reference space used in the
1039*  projected steepest edge pricing algorithm. */
1040
1041static void reset_refsp(struct csa *csa)
1042{     int m = csa->m;
1043      int n = csa->n;
1044      int *head = csa->head;
1045      char *refsp = csa->refsp;
1046      double *gamma = csa->gamma;
1047      int j, k;
1048      xassert(csa->refct == 0);
1049      csa->refct = 1000;
1050      memset(&refsp[1], 0, (m+n) * sizeof(char));
1051      for (j = 1; j <= n; j++)
1052      {  k = head[m+j]; /* x[k] = xN[j] */
1053         refsp[k] = 1;
1054         gamma[j] = 1.0;
1055      }
1056      return;
1057}
1058
1059/***********************************************************************
1060*  eval_gamma - compute steepest edge coefficient
1061*
1062*  This routine computes the steepest edge coefficient for non-basic
1063*  variable xN[j] using its direct definition:
1064*
1065*     gamma[j] = delta[j] +  sum   alfa[i,j]^2,
1066*                           i in R
1067*
1068*  where delta[j] = 1, if xN[j] is in the current reference space,
1069*  and 0 otherwise; R is a set of basic variables xB[i], which are in
1070*  the current reference space; alfa[i,j] are elements of the current
1071*  simplex table.
1072*
1073*  NOTE: The routine is intended only for debugginig purposes. */
1074
1075static double eval_gamma(struct csa *csa, int j)
1076{     int m = csa->m;
1077#ifdef GLP_DEBUG
1078      int n = csa->n;
1079#endif
1080      int *head = csa->head;
1081      char *refsp = csa->refsp;
1082      double *alfa = csa->work3;
1083      double *h = csa->work3;
1084      int i, k;
1085      double gamma;
1086#ifdef GLP_DEBUG
1087      xassert(1 <= j && j <= n);
1088#endif
1089      k = head[m+j]; /* x[k] = xN[j] */
1090#ifdef GLP_DEBUG
1091      xassert(1 <= k && k <= m+n);
1092#endif
1093      /* construct the right-hand side vector h = - N[j] */
1094      for (i = 1; i <= m; i++)
1095         h[i] = 0.0;
1096      if (k <= m)
1097      {  /* N[j] is k-th column of submatrix I */
1098         h[k] = -1.0;
1099      }
1100      else
1101      {  /* N[j] is (k-m)-th column of submatrix (-A) */
1102         int *A_ptr = csa->A_ptr;
1103         int *A_ind = csa->A_ind;
1104         double *A_val = csa->A_val;
1105         int beg, end, ptr;
1106         beg = A_ptr[k-m];
1107         end = A_ptr[k-m+1];
1108         for (ptr = beg; ptr < end; ptr++)
1109            h[A_ind[ptr]] = A_val[ptr];
1110      }
1111      /* solve system B * alfa = h */
1112      xassert(csa->valid);
1113      bfd_ftran(csa->bfd, alfa);
1114      /* compute gamma */
1115      gamma = (refsp[k] ? 1.0 : 0.0);
1116      for (i = 1; i <= m; i++)
1117      {  k = head[i];
1118#ifdef GLP_DEBUG
1119         xassert(1 <= k && k <= m+n);
1120#endif
1121         if (refsp[k]) gamma += alfa[i] * alfa[i];
1122      }
1123      return gamma;
1124}
1125
1126/***********************************************************************
1127*  chuzc - choose non-basic variable (column of the simplex table)
1128*
1129*  This routine chooses non-basic variable xN[q], which has largest
1130*  weighted reduced cost:
1131*
1132*     |d[q]| / sqrt(gamma[q]) = max  |d[j]| / sqrt(gamma[j]),
1133*                              j in J
1134*
1135*  where J is a subset of eligible non-basic variables xN[j], d[j] is
1136*  reduced cost of xN[j], gamma[j] is the steepest edge coefficient.
1137*
1138*  The working objective function is always minimized, so the sign of
1139*  d[q] determines direction, in which xN[q] has to change:
1140*
1141*     if d[q] < 0, xN[q] has to increase;
1142*
1143*     if d[q] > 0, xN[q] has to decrease.
1144*
1145*  If |d[j]| <= tol_dj, where tol_dj is a specified tolerance, xN[j]
1146*  is not included in J and therefore ignored. (It is assumed that the
1147*  working objective row is appropriately scaled, i.e. max|c[k]| = 1.)
1148*
1149*  If J is empty and no variable has been chosen, q is set to 0. */
1150
1151static void chuzc(struct csa *csa, double tol_dj)
1152{     int n = csa->n;
1153      char *stat = csa->stat;
1154      double *cbar = csa->cbar;
1155      double *gamma = csa->gamma;
1156      int j, q;
1157      double dj, best, temp;
1158      /* nothing is chosen so far */
1159      q = 0, best = 0.0;
1160      /* look through the list of non-basic variables */
1161      for (j = 1; j <= n; j++)
1162      {  dj = cbar[j];
1163         switch (stat[j])
1164         {  case GLP_NL:
1165               /* xN[j] can increase */
1166               if (dj >= - tol_dj) continue;
1167               break;
1168            case GLP_NU:
1169               /* xN[j] can decrease */
1170               if (dj <= + tol_dj) continue;
1171               break;
1172            case GLP_NF:
1173               /* xN[j] can change in any direction */
1174               if (- tol_dj <= dj && dj <= + tol_dj) continue;
1175               break;
1176            case GLP_NS:
1177               /* xN[j] cannot change at all */
1178               continue;
1179            default:
1180               xassert(stat != stat);
1181         }
1182         /* xN[j] is eligible non-basic variable; choose one which has
1183            largest weighted reduced cost */
1184#ifdef GLP_DEBUG
1185         xassert(gamma[j] > 0.0);
1186#endif
1187         temp = (dj * dj) / gamma[j];
1188         if (best < temp)
1189            q = j, best = temp;
1190      }
1191      /* store the index of non-basic variable xN[q] chosen */
1192      csa->q = q;
1193      return;
1194}
1195
1196/***********************************************************************
1197*  eval_tcol - compute pivot column of the simplex table
1198*
1199*  This routine computes the pivot column of the simplex table, which
1200*  corresponds to non-basic variable xN[q] chosen.
1201*
1202*  The pivot column is the following vector:
1203*
1204*     tcol = T * e[q] = - inv(B) * N * e[q] = - inv(B) * N[q],
1205*
1206*  where B is the current basis matrix, N[q] is a column of the matrix
1207*  (I|-A) corresponding to variable xN[q]. */
1208
1209static void eval_tcol(struct csa *csa)
1210{     int m = csa->m;
1211#ifdef GLP_DEBUG
1212      int n = csa->n;
1213#endif
1214      int *head = csa->head;
1215      int q = csa->q;
1216      int *tcol_ind = csa->tcol_ind;
1217      double *tcol_vec = csa->tcol_vec;
1218      double *h = csa->tcol_vec;
1219      int i, k, nnz;
1220#ifdef GLP_DEBUG
1221      xassert(1 <= q && q <= n);
1222#endif
1223      k = head[m+q]; /* x[k] = xN[q] */
1224#ifdef GLP_DEBUG
1225      xassert(1 <= k && k <= m+n);
1226#endif
1227      /* construct the right-hand side vector h = - N[q] */
1228      for (i = 1; i <= m; i++)
1229         h[i] = 0.0;
1230      if (k <= m)
1231      {  /* N[q] is k-th column of submatrix I */
1232         h[k] = -1.0;
1233      }
1234      else
1235      {  /* N[q] is (k-m)-th column of submatrix (-A) */
1236         int *A_ptr = csa->A_ptr;
1237         int *A_ind = csa->A_ind;
1238         double *A_val = csa->A_val;
1239         int beg, end, ptr;
1240         beg = A_ptr[k-m];
1241         end = A_ptr[k-m+1];
1242         for (ptr = beg; ptr < end; ptr++)
1243            h[A_ind[ptr]] = A_val[ptr];
1244      }
1245      /* solve system B * tcol = h */
1246      xassert(csa->valid);
1247      bfd_ftran(csa->bfd, tcol_vec);
1248      /* construct sparse pattern of the pivot column */
1249      nnz = 0;
1250      for (i = 1; i <= m; i++)
1251      {  if (tcol_vec[i] != 0.0)
1252            tcol_ind[++nnz] = i;
1253      }
1254      csa->tcol_nnz = nnz;
1255      return;
1256}
1257
1258/***********************************************************************
1259*  refine_tcol - refine pivot column of the simplex table
1260*
1261*  This routine refines the pivot column of the simplex table assuming
1262*  that it was previously computed by the routine eval_tcol. */
1263
1264static void refine_tcol(struct csa *csa)
1265{     int m = csa->m;
1266#ifdef GLP_DEBUG
1267      int n = csa->n;
1268#endif
1269      int *head = csa->head;
1270      int q = csa->q;
1271      int *tcol_ind = csa->tcol_ind;
1272      double *tcol_vec = csa->tcol_vec;
1273      double *h = csa->work3;
1274      int i, k, nnz;
1275#ifdef GLP_DEBUG
1276      xassert(1 <= q && q <= n);
1277#endif
1278      k = head[m+q]; /* x[k] = xN[q] */
1279#ifdef GLP_DEBUG
1280      xassert(1 <= k && k <= m+n);
1281#endif
1282      /* construct the right-hand side vector h = - N[q] */
1283      for (i = 1; i <= m; i++)
1284         h[i] = 0.0;
1285      if (k <= m)
1286      {  /* N[q] is k-th column of submatrix I */
1287         h[k] = -1.0;
1288      }
1289      else
1290      {  /* N[q] is (k-m)-th column of submatrix (-A) */
1291         int *A_ptr = csa->A_ptr;
1292         int *A_ind = csa->A_ind;
1293         double *A_val = csa->A_val;
1294         int beg, end, ptr;
1295         beg = A_ptr[k-m];
1296         end = A_ptr[k-m+1];
1297         for (ptr = beg; ptr < end; ptr++)
1298            h[A_ind[ptr]] = A_val[ptr];
1299      }
1300      /* refine solution of B * tcol = h */
1301      refine_ftran(csa, h, tcol_vec);
1302      /* construct sparse pattern of the pivot column */
1303      nnz = 0;
1304      for (i = 1; i <= m; i++)
1305      {  if (tcol_vec[i] != 0.0)
1306            tcol_ind[++nnz] = i;
1307      }
1308      csa->tcol_nnz = nnz;
1309      return;
1310}
1311
1312/***********************************************************************
1313*  sort_tcol - sort pivot column of the simplex table
1314*
1315*  This routine reorders the list of non-zero elements of the pivot
1316*  column to put significant elements, whose magnitude is not less than
1317*  a specified tolerance, in front of the list, and stores the number
1318*  of significant elements in tcol_num. */
1319
1320static void sort_tcol(struct csa *csa, double tol_piv)
1321{
1322#ifdef GLP_DEBUG
1323      int m = csa->m;
1324#endif
1325      int nnz = csa->tcol_nnz;
1326      int *tcol_ind = csa->tcol_ind;
1327      double *tcol_vec = csa->tcol_vec;
1328      int i, num, pos;
1329      double big, eps, temp;
1330      /* compute infinity (maximum) norm of the column */
1331      big = 0.0;
1332      for (pos = 1; pos <= nnz; pos++)
1333      {
1334#ifdef GLP_DEBUG
1335         i = tcol_ind[pos];
1336         xassert(1 <= i && i <= m);
1337#endif
1338         temp = fabs(tcol_vec[tcol_ind[pos]]);
1339         if (big < temp) big = temp;
1340      }
1341      csa->tcol_max = big;
1342      /* determine absolute pivot tolerance */
1343      eps = tol_piv * (1.0 + 0.01 * big);
1344      /* move significant column components to front of the list */
1345      for (num = 0; num < nnz; )
1346      {  i = tcol_ind[nnz];
1347         if (fabs(tcol_vec[i]) < eps)
1348            nnz--;
1349         else
1350         {  num++;
1351            tcol_ind[nnz] = tcol_ind[num];
1352            tcol_ind[num] = i;
1353         }
1354      }
1355      csa->tcol_num = num;
1356      return;
1357}
1358
1359/***********************************************************************
1360*  chuzr - choose basic variable (row of the simplex table)
1361*
1362*  This routine chooses basic variable xB[p], which reaches its bound
1363*  first on changing non-basic variable xN[q] in valid direction.
1364*
1365*  The parameter rtol is a relative tolerance used to relax bounds of
1366*  basic variables. If rtol = 0, the routine implements the standard
1367*  ratio test. Otherwise, if rtol > 0, the routine implements Harris'
1368*  two-pass ratio test. In the latter case rtol should be about three
1369*  times less than a tolerance used to check primal feasibility. */
1370
1371static void chuzr(struct csa *csa, double rtol)
1372{     int m = csa->m;
1373#ifdef GLP_DEBUG
1374      int n = csa->n;
1375#endif
1376      char *type = csa->type;
1377      double *lb = csa->lb;
1378      double *ub = csa->ub;
1379      double *coef = csa->coef;
1380      int *head = csa->head;
1381      int phase = csa->phase;
1382      double *bbar = csa->bbar;
1383      double *cbar = csa->cbar;
1384      int q = csa->q;
1385      int *tcol_ind = csa->tcol_ind;
1386      double *tcol_vec = csa->tcol_vec;
1387      int tcol_num = csa->tcol_num;
1388      int i, i_stat, k, p, p_stat, pos;
1389      double alfa, big, delta, s, t, teta, tmax;
1390#ifdef GLP_DEBUG
1391      xassert(1 <= q && q <= n);
1392#endif
1393      /* s := - sign(d[q]), where d[q] is reduced cost of xN[q] */
1394#ifdef GLP_DEBUG
1395      xassert(cbar[q] != 0.0);
1396#endif
1397      s = (cbar[q] > 0.0 ? -1.0 : +1.0);
1398      /*** FIRST PASS ***/
1399      k = head[m+q]; /* x[k] = xN[q] */
1400#ifdef GLP_DEBUG
1401      xassert(1 <= k && k <= m+n);
1402#endif
1403      if (type[k] == GLP_DB)
1404      {  /* xN[q] has both lower and upper bounds */
1405         p = -1, p_stat = 0, teta = ub[k] - lb[k], big = 1.0;
1406      }
1407      else
1408      {  /* xN[q] has no opposite bound */
1409         p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
1410      }
1411      /* walk through significant elements of the pivot column */
1412      for (pos = 1; pos <= tcol_num; pos++)
1413      {  i = tcol_ind[pos];
1414#ifdef GLP_DEBUG
1415         xassert(1 <= i && i <= m);
1416#endif
1417         k = head[i]; /* x[k] = xB[i] */
1418#ifdef GLP_DEBUG
1419         xassert(1 <= k && k <= m+n);
1420#endif
1421         alfa = s * tcol_vec[i];
1422#ifdef GLP_DEBUG
1423         xassert(alfa != 0.0);
1424#endif
1425         /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to
1426            consider the only case when xN[q] is increasing */
1427         if (alfa > 0.0)
1428         {  /* xB[i] is increasing */
1429            if (phase == 1 && coef[k] < 0.0)
1430            {  /* xB[i] violates its lower bound, which plays the role
1431                  of an upper bound on phase I */
1432               delta = rtol * (1.0 + kappa * fabs(lb[k]));
1433               t = ((lb[k] + delta) - bbar[i]) / alfa;
1434               i_stat = GLP_NL;
1435            }
1436            else if (phase == 1 && coef[k] > 0.0)
1437            {  /* xB[i] violates its upper bound, which plays the role
1438                  of an lower bound on phase I */
1439               continue;
1440            }
1441            else if (type[k] == GLP_UP || type[k] == GLP_DB ||
1442                     type[k] == GLP_FX)
1443            {  /* xB[i] is within its bounds and has an upper bound */
1444               delta = rtol * (1.0 + kappa * fabs(ub[k]));
1445               t = ((ub[k] + delta) - bbar[i]) / alfa;
1446               i_stat = GLP_NU;
1447            }
1448            else
1449            {  /* xB[i] is within its bounds and has no upper bound */
1450               continue;
1451            }
1452         }
1453         else
1454         {  /* xB[i] is decreasing */
1455            if (phase == 1 && coef[k] > 0.0)
1456            {  /* xB[i] violates its upper bound, which plays the role
1457                  of an lower bound on phase I */
1458               delta = rtol * (1.0 + kappa * fabs(ub[k]));
1459               t = ((ub[k] - delta) - bbar[i]) / alfa;
1460               i_stat = GLP_NU;
1461            }
1462            else if (phase == 1 && coef[k] < 0.0)
1463            {  /* xB[i] violates its lower bound, which plays the role
1464                  of an upper bound on phase I */
1465               continue;
1466            }
1467            else if (type[k] == GLP_LO || type[k] == GLP_DB ||
1468                     type[k] == GLP_FX)
1469            {  /* xB[i] is within its bounds and has an lower bound */
1470               delta = rtol * (1.0 + kappa * fabs(lb[k]));
1471               t = ((lb[k] - delta) - bbar[i]) / alfa;
1472               i_stat = GLP_NL;
1473            }
1474            else
1475            {  /* xB[i] is within its bounds and has no lower bound */
1476               continue;
1477            }
1478         }
1479         /* t is a change of xN[q], on which xB[i] reaches its bound
1480            (possibly relaxed); since the basic solution is assumed to
1481            be primal feasible (or pseudo feasible on phase I), t has
1482            to be non-negative by definition; however, it may happen
1483            that xB[i] slightly (i.e. within a tolerance) violates its
1484            bound, that leads to negative t; in the latter case, if
1485            xB[i] is chosen, negative t means that xN[q] changes in
1486            wrong direction; if pivot alfa[i,q] is close to zero, even
1487            small bound violation of xB[i] may lead to a large change
1488            of xN[q] in wrong direction; let, for example, xB[i] >= 0
1489            and in the current basis its value be -5e-9; let also xN[q]
1490            be on its zero bound and should increase; from the ratio
1491            test rule it follows that the pivot alfa[i,q] < 0; however,
1492            if alfa[i,q] is, say, -1e-9, the change of xN[q] in wrong
1493            direction is 5e-9 / (-1e-9) = -5, and using it for updating
1494            values of other basic variables will give absolutely wrong
1495            results; therefore, if t is negative, we should replace it
1496            by exact zero assuming that xB[i] is exactly on its bound,
1497            and the violation appears due to round-off errors */
1498         if (t < 0.0) t = 0.0;
1499         /* apply minimal ratio test */
1500         if (teta > t || teta == t && big < fabs(alfa))
1501            p = i, p_stat = i_stat, teta = t, big = fabs(alfa);
1502      }
1503      /* the second pass is skipped in the following cases: */
1504      /* if the standard ratio test is used */
1505      if (rtol == 0.0) goto done;
1506      /* if xN[q] reaches its opposite bound or if no basic variable
1507         has been chosen on the first pass */
1508      if (p <= 0) goto done;
1509      /* if xB[p] is a blocking variable, i.e. if it prevents xN[q]
1510         from any change */
1511      if (teta == 0.0) goto done;
1512      /*** SECOND PASS ***/
1513      /* here tmax is a maximal change of xN[q], on which the solution
1514         remains primal feasible (or pseudo feasible on phase I) within
1515         a tolerance */
1516#if 0
1517      tmax = (1.0 + 10.0 * DBL_EPSILON) * teta;
1518#else
1519      tmax = teta;
1520#endif
1521      /* nothing is chosen so far */
1522      p = 0, p_stat = 0, teta = DBL_MAX, big = 0.0;
1523      /* walk through significant elements of the pivot column */
1524      for (pos = 1; pos <= tcol_num; pos++)
1525      {  i = tcol_ind[pos];
1526#ifdef GLP_DEBUG
1527         xassert(1 <= i && i <= m);
1528#endif
1529         k = head[i]; /* x[k] = xB[i] */
1530#ifdef GLP_DEBUG
1531         xassert(1 <= k && k <= m+n);
1532#endif
1533         alfa = s * tcol_vec[i];
1534#ifdef GLP_DEBUG
1535         xassert(alfa != 0.0);
1536#endif
1537         /* xB[i] = ... + alfa * xN[q] + ..., and due to s we need to
1538            consider the only case when xN[q] is increasing */
1539         if (alfa > 0.0)
1540         {  /* xB[i] is increasing */
1541            if (phase == 1 && coef[k] < 0.0)
1542            {  /* xB[i] violates its lower bound, which plays the role
1543                  of an upper bound on phase I */
1544               t = (lb[k] - bbar[i]) / alfa;
1545               i_stat = GLP_NL;
1546            }
1547            else if (phase == 1 && coef[k] > 0.0)
1548            {  /* xB[i] violates its upper bound, which plays the role
1549                  of an lower bound on phase I */
1550               continue;
1551            }
1552            else if (type[k] == GLP_UP || type[k] == GLP_DB ||
1553                     type[k] == GLP_FX)
1554            {  /* xB[i] is within its bounds and has an upper bound */
1555               t = (ub[k] - bbar[i]) / alfa;
1556               i_stat = GLP_NU;
1557            }
1558            else
1559            {  /* xB[i] is within its bounds and has no upper bound */
1560               continue;
1561            }
1562         }
1563         else
1564         {  /* xB[i] is decreasing */
1565            if (phase == 1 && coef[k] > 0.0)
1566            {  /* xB[i] violates its upper bound, which plays the role
1567                  of an lower bound on phase I */
1568               t = (ub[k] - bbar[i]) / alfa;
1569               i_stat = GLP_NU;
1570            }
1571            else if (phase == 1 && coef[k] < 0.0)
1572            {  /* xB[i] violates its lower bound, which plays the role
1573                  of an upper bound on phase I */
1574               continue;
1575            }
1576            else if (type[k] == GLP_LO || type[k] == GLP_DB ||
1577                     type[k] == GLP_FX)
1578            {  /* xB[i] is within its bounds and has an lower bound */
1579               t = (lb[k] - bbar[i]) / alfa;
1580               i_stat = GLP_NL;
1581            }
1582            else
1583            {  /* xB[i] is within its bounds and has no lower bound */
1584               continue;
1585            }
1586         }
1587         /* (see comments for the first pass) */
1588         if (t < 0.0) t = 0.0;
1589         /* t is a change of xN[q], on which xB[i] reaches its bound;
1590            if t <= tmax, all basic variables can violate their bounds
1591            only within relaxation tolerance delta; we can use this
1592            freedom and choose basic variable having largest influence
1593            coefficient to avoid possible numeric instability */
1594         if (t <= tmax && big < fabs(alfa))
1595            p = i, p_stat = i_stat, teta = t, big = fabs(alfa);
1596      }
1597      /* something must be chosen on the second pass */
1598      xassert(p != 0);
1599done: /* store the index and status of basic variable xB[p] chosen */
1600      csa->p = p;
1601      if (p > 0 && type[head[p]] == GLP_FX)
1602         csa->p_stat = GLP_NS;
1603      else
1604         csa->p_stat = p_stat;
1605      /* store corresponding change of non-basic variable xN[q] */
1606#ifdef GLP_DEBUG
1607      xassert(teta >= 0.0);
1608#endif
1609      csa->teta = s * teta;
1610      return;
1611}
1612
1613/***********************************************************************
1614*  eval_rho - compute pivot row of the inverse
1615*
1616*  This routine computes the pivot (p-th) row of the inverse inv(B),
1617*  which corresponds to basic variable xB[p] chosen:
1618*
1619*     rho = inv(B') * e[p],
1620*
1621*  where B' is a matrix transposed to the current basis matrix, e[p]
1622*  is unity vector. */
1623
1624static void eval_rho(struct csa *csa, double rho[])
1625{     int m = csa->m;
1626      int p = csa->p;
1627      double *e = rho;
1628      int i;
1629#ifdef GLP_DEBUG
1630      xassert(1 <= p && p <= m);
1631#endif
1632      /* construct the right-hand side vector e[p] */
1633      for (i = 1; i <= m; i++)
1634         e[i] = 0.0;
1635      e[p] = 1.0;
1636      /* solve system B'* rho = e[p] */
1637      xassert(csa->valid);
1638      bfd_btran(csa->bfd, rho);
1639      return;
1640}
1641
1642/***********************************************************************
1643*  refine_rho - refine pivot row of the inverse
1644*
1645*  This routine refines the pivot row of the inverse inv(B) assuming
1646*  that it was previously computed by the routine eval_rho. */
1647
1648static void refine_rho(struct csa *csa, double rho[])
1649{     int m = csa->m;
1650      int p = csa->p;
1651      double *e = csa->work3;
1652      int i;
1653#ifdef GLP_DEBUG
1654      xassert(1 <= p && p <= m);
1655#endif
1656      /* construct the right-hand side vector e[p] */
1657      for (i = 1; i <= m; i++)
1658         e[i] = 0.0;
1659      e[p] = 1.0;
1660      /* refine solution of B'* rho = e[p] */
1661      refine_btran(csa, e, rho);
1662      return;
1663}
1664
1665/***********************************************************************
1666*  eval_trow - compute pivot row of the simplex table
1667*
1668*  This routine computes the pivot row of the simplex table, which
1669*  corresponds to basic variable xB[p] chosen.
1670*
1671*  The pivot row is the following vector:
1672*
1673*     trow = T'* e[p] = - N'* inv(B') * e[p] = - N' * rho,
1674*
1675*  where rho is the pivot row of the inverse inv(B) previously computed
1676*  by the routine eval_rho.
1677*
1678*  Note that elements of the pivot row corresponding to fixed non-basic
1679*  variables are not computed. */
1680
1681static void eval_trow(struct csa *csa, double rho[])
1682{     int m = csa->m;
1683      int n = csa->n;
1684#ifdef GLP_DEBUG
1685      char *stat = csa->stat;
1686#endif
1687      int *N_ptr = csa->N_ptr;
1688      int *N_len = csa->N_len;
1689      int *N_ind = csa->N_ind;
1690      double *N_val = csa->N_val;
1691      int *trow_ind = csa->trow_ind;
1692      double *trow_vec = csa->trow_vec;
1693      int i, j, beg, end, ptr, nnz;
1694      double temp;
1695      /* clear the pivot row */
1696      for (j = 1; j <= n; j++)
1697         trow_vec[j] = 0.0;
1698      /* compute the pivot row as a linear combination of rows of the
1699         matrix N: trow = - rho[1] * N'[1] - ... - rho[m] * N'[m] */
1700      for (i = 1; i <= m; i++)
1701      {  temp = rho[i];
1702         if (temp == 0.0) continue;
1703         /* trow := trow - rho[i] * N'[i] */
1704         beg = N_ptr[i];
1705         end = beg + N_len[i];
1706         for (ptr = beg; ptr < end; ptr++)
1707         {
1708#ifdef GLP_DEBUG
1709            j = N_ind[ptr];
1710            xassert(1 <= j && j <= n);
1711            xassert(stat[j] != GLP_NS);
1712#endif
1713            trow_vec[N_ind[ptr]] -= temp * N_val[ptr];
1714         }
1715      }
1716      /* construct sparse pattern of the pivot row */
1717      nnz = 0;
1718      for (j = 1; j <= n; j++)
1719      {  if (trow_vec[j] != 0.0)
1720            trow_ind[++nnz] = j;
1721      }
1722      csa->trow_nnz = nnz;
1723      return;
1724}
1725
1726/***********************************************************************
1727*  update_bbar - update values of basic variables
1728*
1729*  This routine updates values of all basic variables for the adjacent
1730*  basis. */
1731
1732static void update_bbar(struct csa *csa)
1733{
1734#ifdef GLP_DEBUG
1735      int m = csa->m;
1736      int n = csa->n;
1737#endif
1738      double *bbar = csa->bbar;
1739      int q = csa->q;
1740      int tcol_nnz = csa->tcol_nnz;
1741      int *tcol_ind = csa->tcol_ind;
1742      double *tcol_vec = csa->tcol_vec;
1743      int p = csa->p;
1744      double teta = csa->teta;
1745      int i, pos;
1746#ifdef GLP_DEBUG
1747      xassert(1 <= q && q <= n);
1748      xassert(p < 0 || 1 <= p && p <= m);
1749#endif
1750      /* if xN[q] leaves the basis, compute its value in the adjacent
1751         basis, where it will replace xB[p] */
1752      if (p > 0)
1753         bbar[p] = get_xN(csa, q) + teta;
1754      /* update values of other basic variables (except xB[p], because
1755         it will be replaced by xN[q]) */
1756      if (teta == 0.0) goto done;
1757      for (pos = 1; pos <= tcol_nnz; pos++)
1758      {  i = tcol_ind[pos];
1759         /* skip xB[p] */
1760         if (i == p) continue;
1761         /* (change of xB[i]) = alfa[i,q] * (change of xN[q]) */
1762         bbar[i] += tcol_vec[i] * teta;
1763      }
1764done: return;
1765}
1766
1767/***********************************************************************
1768*  reeval_cost - recompute reduced cost of non-basic variable xN[q]
1769*
1770*  This routine recomputes reduced cost of non-basic variable xN[q] for
1771*  the current basis more accurately using its direct definition:
1772*
1773*     d[q] = cN[q] - N'[q] * pi =
1774*
1775*          = cN[q] - N'[q] * (inv(B') * cB) =
1776*
1777*          = cN[q] - (cB' * inv(B) * N[q]) =
1778*
1779*          = cN[q] + cB' * (pivot column).
1780*
1781*  It is assumed that the pivot column of the simplex table is already
1782*  computed. */
1783
1784static double reeval_cost(struct csa *csa)
1785{     int m = csa->m;
1786#ifdef GLP_DEBUG
1787      int n = csa->n;
1788#endif
1789      double *coef = csa->coef;
1790      int *head = csa->head;
1791      int q = csa->q;
1792      int tcol_nnz = csa->tcol_nnz;
1793      int *tcol_ind = csa->tcol_ind;
1794      double *tcol_vec = csa->tcol_vec;
1795      int i, pos;
1796      double dq;
1797#ifdef GLP_DEBUG
1798      xassert(1 <= q && q <= n);
1799#endif
1800      dq = coef[head[m+q]];
1801      for (pos = 1; pos <= tcol_nnz; pos++)
1802      {  i = tcol_ind[pos];
1803#ifdef GLP_DEBUG
1804         xassert(1 <= i && i <= m);
1805#endif
1806         dq += coef[head[i]] * tcol_vec[i];
1807      }
1808      return dq;
1809}
1810
1811/***********************************************************************
1812*  update_cbar - update reduced costs of non-basic variables
1813*
1814*  This routine updates reduced costs of all (except fixed) non-basic
1815*  variables for the adjacent basis. */
1816
1817static void update_cbar(struct csa *csa)
1818{
1819#ifdef GLP_DEBUG
1820      int n = csa->n;
1821#endif
1822      double *cbar = csa->cbar;
1823      int q = csa->q;
1824      int trow_nnz = csa->trow_nnz;
1825      int *trow_ind = csa->trow_ind;
1826      double *trow_vec = csa->trow_vec;
1827      int j, pos;
1828      double new_dq;
1829#ifdef GLP_DEBUG
1830      xassert(1 <= q && q <= n);
1831#endif
1832      /* compute reduced cost of xB[p] in the adjacent basis, where it
1833         will replace xN[q] */
1834#ifdef GLP_DEBUG
1835      xassert(trow_vec[q] != 0.0);
1836#endif
1837      new_dq = (cbar[q] /= trow_vec[q]);
1838      /* update reduced costs of other non-basic variables (except
1839         xN[q], because it will be replaced by xB[p]) */
1840      for (pos = 1; pos <= trow_nnz; pos++)
1841      {  j = trow_ind[pos];
1842         /* skip xN[q] */
1843         if (j == q) continue;
1844         cbar[j] -= trow_vec[j] * new_dq;
1845      }
1846      return;
1847}
1848
1849/***********************************************************************
1850*  update_gamma - update steepest edge coefficients
1851*
1852*  This routine updates steepest-edge coefficients for the adjacent
1853*  basis. */
1854
1855static void update_gamma(struct csa *csa)
1856{     int m = csa->m;
1857#ifdef GLP_DEBUG
1858      int n = csa->n;
1859#endif
1860      char *type = csa->type;
1861      int *A_ptr = csa->A_ptr;
1862      int *A_ind = csa->A_ind;
1863      double *A_val = csa->A_val;
1864      int *head = csa->head;
1865      char *refsp = csa->refsp;
1866      double *gamma = csa->gamma;
1867      int q = csa->q;
1868      int tcol_nnz = csa->tcol_nnz;
1869      int *tcol_ind = csa->tcol_ind;
1870      double *tcol_vec = csa->tcol_vec;
1871      int p = csa->p;
1872      int trow_nnz = csa->trow_nnz;
1873      int *trow_ind = csa->trow_ind;
1874      double *trow_vec = csa->trow_vec;
1875      double *u = csa->work3;
1876      int i, j, k, pos, beg, end, ptr;
1877      double gamma_q, delta_q, pivot, s, t, t1, t2;
1878#ifdef GLP_DEBUG
1879      xassert(1 <= p && p <= m);
1880      xassert(1 <= q && q <= n);
1881#endif
1882      /* the basis changes, so decrease the count */
1883      xassert(csa->refct > 0);
1884      csa->refct--;
1885      /* recompute gamma[q] for the current basis more accurately and
1886         compute auxiliary vector u */
1887      gamma_q = delta_q = (refsp[head[m+q]] ? 1.0 : 0.0);
1888      for (i = 1; i <= m; i++) u[i] = 0.0;
1889      for (pos = 1; pos <= tcol_nnz; pos++)
1890      {  i = tcol_ind[pos];
1891         if (refsp[head[i]])
1892         {  u[i] = t = tcol_vec[i];
1893            gamma_q += t * t;
1894         }
1895         else
1896            u[i] = 0.0;
1897      }
1898      xassert(csa->valid);
1899      bfd_btran(csa->bfd, u);
1900      /* update gamma[k] for other non-basic variables (except fixed
1901         variables and xN[q], because it will be replaced by xB[p]) */
1902      pivot = trow_vec[q];
1903#ifdef GLP_DEBUG
1904      xassert(pivot != 0.0);
1905#endif
1906      for (pos = 1; pos <= trow_nnz; pos++)
1907      {  j = trow_ind[pos];
1908         /* skip xN[q] */
1909         if (j == q) continue;
1910         /* compute t */
1911         t = trow_vec[j] / pivot;
1912         /* compute inner product s = N'[j] * u */
1913         k = head[m+j]; /* x[k] = xN[j] */
1914         if (k <= m)
1915            s = u[k];
1916         else
1917         {  s = 0.0;
1918            beg = A_ptr[k-m];
1919            end = A_ptr[k-m+1];
1920            for (ptr = beg; ptr < end; ptr++)
1921               s -= A_val[ptr] * u[A_ind[ptr]];
1922         }
1923         /* compute gamma[k] for the adjacent basis */
1924         t1 = gamma[j] + t * t * gamma_q + 2.0 * t * s;
1925         t2 = (refsp[k] ? 1.0 : 0.0) + delta_q * t * t;
1926         gamma[j] = (t1 >= t2 ? t1 : t2);
1927         if (gamma[j] < DBL_EPSILON) gamma[j] = DBL_EPSILON;
1928      }
1929      /* compute gamma[q] for the adjacent basis */
1930      if (type[head[p]] == GLP_FX)
1931         gamma[q] = 1.0;
1932      else
1933      {  gamma[q] = gamma_q / (pivot * pivot);
1934         if (gamma[q] < DBL_EPSILON) gamma[q] = DBL_EPSILON;
1935      }
1936      return;
1937}
1938
1939/***********************************************************************
1940*  err_in_bbar - compute maximal relative error in primal solution
1941*
1942*  This routine returns maximal relative error:
1943*
1944*     max |beta[i] - bbar[i]| / (1 + |beta[i]|),
1945*
1946*  where beta and bbar are, respectively, directly computed and the
1947*  current (updated) values of basic variables.
1948*
1949*  NOTE: The routine is intended only for debugginig purposes. */
1950
1951static double err_in_bbar(struct csa *csa)
1952{     int m = csa->m;
1953      double *bbar = csa->bbar;
1954      int i;
1955      double e, emax, *beta;
1956      beta = xcalloc(1+m, sizeof(double));
1957      eval_beta(csa, beta);
1958      emax = 0.0;
1959      for (i = 1; i <= m; i++)
1960      {  e = fabs(beta[i] - bbar[i]) / (1.0 + fabs(beta[i]));
1961         if (emax < e) emax = e;
1962      }
1963      xfree(beta);
1964      return emax;
1965}
1966
1967/***********************************************************************
1968*  err_in_cbar - compute maximal relative error in dual solution
1969*
1970*  This routine returns maximal relative error:
1971*
1972*     max |cost[j] - cbar[j]| / (1 + |cost[j]|),
1973*
1974*  where cost and cbar are, respectively, directly computed and the
1975*  current (updated) reduced costs of non-basic non-fixed variables.
1976*
1977*  NOTE: The routine is intended only for debugginig purposes. */
1978
1979static double err_in_cbar(struct csa *csa)
1980{     int m = csa->m;
1981      int n = csa->n;
1982      char *stat = csa->stat;
1983      double *cbar = csa->cbar;
1984      int j;
1985      double e, emax, cost, *pi;
1986      pi = xcalloc(1+m, sizeof(double));
1987      eval_pi(csa, pi);
1988      emax = 0.0;
1989      for (j = 1; j <= n; j++)
1990      {  if (stat[j] == GLP_NS) continue;
1991         cost = eval_cost(csa, pi, j);
1992         e = fabs(cost - cbar[j]) / (1.0 + fabs(cost));
1993         if (emax < e) emax = e;
1994      }
1995      xfree(pi);
1996      return emax;
1997}
1998
1999/***********************************************************************
2000*  err_in_gamma - compute maximal relative error in steepest edge cff.
2001*
2002*  This routine returns maximal relative error:
2003*
2004*     max |gamma'[j] - gamma[j]| / (1 + |gamma'[j]),
2005*
2006*  where gamma'[j] and gamma[j] are, respectively, directly computed
2007*  and the current (updated) steepest edge coefficients for non-basic
2008*  non-fixed variable x[j].
2009*
2010*  NOTE: The routine is intended only for debugginig purposes. */
2011
2012static double err_in_gamma(struct csa *csa)
2013{     int n = csa->n;
2014      char *stat = csa->stat;
2015      double *gamma = csa->gamma;
2016      int j;
2017      double e, emax, temp;
2018      emax = 0.0;
2019      for (j = 1; j <= n; j++)
2020      {  if (stat[j] == GLP_NS)
2021         {  xassert(gamma[j] == 1.0);
2022            continue;
2023         }
2024         temp = eval_gamma(csa, j);
2025         e = fabs(temp - gamma[j]) / (1.0 + fabs(temp));
2026         if (emax < e) emax = e;
2027      }
2028      return emax;
2029}
2030
2031/***********************************************************************
2032*  change_basis - change basis header
2033*
2034*  This routine changes the basis header to make it corresponding to
2035*  the adjacent basis. */
2036
2037static void change_basis(struct csa *csa)
2038{     int m = csa->m;
2039#ifdef GLP_DEBUG
2040      int n = csa->n;
2041      char *type = csa->type;
2042#endif
2043      int *head = csa->head;
2044      char *stat = csa->stat;
2045      int q = csa->q;
2046      int p = csa->p;
2047      int p_stat = csa->p_stat;
2048      int k;
2049#ifdef GLP_DEBUG
2050      xassert(1 <= q && q <= n);
2051#endif
2052      if (p < 0)
2053      {  /* xN[q] goes to its opposite bound */
2054#ifdef GLP_DEBUG
2055         k = head[m+q]; /* x[k] = xN[q] */
2056         xassert(1 <= k && k <= m+n);
2057         xassert(type[k] == GLP_DB);
2058#endif
2059         switch (stat[q])
2060         {  case GLP_NL:
2061               /* xN[q] increases */
2062               stat[q] = GLP_NU;
2063               break;
2064            case GLP_NU:
2065               /* xN[q] decreases */
2066               stat[q] = GLP_NL;
2067               break;
2068            default:
2069               xassert(stat != stat);
2070         }
2071      }
2072      else
2073      {  /* xB[p] leaves the basis, xN[q] enters the basis */
2074#ifdef GLP_DEBUG
2075         xassert(1 <= p && p <= m);
2076         k = head[p]; /* x[k] = xB[p] */
2077         switch (p_stat)
2078         {  case GLP_NL:
2079               /* xB[p] goes to its lower bound */
2080               xassert(type[k] == GLP_LO || type[k] == GLP_DB);
2081               break;
2082            case GLP_NU:
2083               /* xB[p] goes to its upper bound */
2084               xassert(type[k] == GLP_UP || type[k] == GLP_DB);
2085               break;
2086            case GLP_NS:
2087               /* xB[p] goes to its fixed value */
2088               xassert(type[k] == GLP_NS);
2089               break;
2090            default:
2091               xassert(p_stat != p_stat);
2092         }
2093#endif
2094         /* xB[p] <-> xN[q] */
2095         k = head[p], head[p] = head[m+q], head[m+q] = k;
2096         stat[q] = (char)p_stat;
2097      }
2098      return;
2099}
2100
2101/***********************************************************************
2102*  set_aux_obj - construct auxiliary objective function
2103*
2104*  The auxiliary objective function is a separable piecewise linear
2105*  convex function, which is the sum of primal infeasibilities:
2106*
2107*     z = t[1] + ... + t[m+n] -> minimize,
2108*
2109*  where:
2110*
2111*            / lb[k] - x[k], if x[k] < lb[k]
2112*            |
2113*     t[k] = <  0, if lb[k] <= x[k] <= ub[k]
2114*            |
2115*            \ x[k] - ub[k], if x[k] > ub[k]
2116*
2117*  This routine computes objective coefficients for the current basis
2118*  and returns the number of non-zero terms t[k]. */
2119
2120static int set_aux_obj(struct csa *csa, double tol_bnd)
2121{     int m = csa->m;
2122      int n = csa->n;
2123      char *type = csa->type;
2124      double *lb = csa->lb;
2125      double *ub = csa->ub;
2126      double *coef = csa->coef;
2127      int *head = csa->head;
2128      double *bbar = csa->bbar;
2129      int i, k, cnt = 0;
2130      double eps;
2131      /* use a bit more restrictive tolerance */
2132      tol_bnd *= 0.90;
2133      /* clear all objective coefficients */
2134      for (k = 1; k <= m+n; k++)
2135         coef[k] = 0.0;
2136      /* walk through the list of basic variables */
2137      for (i = 1; i <= m; i++)
2138      {  k = head[i]; /* x[k] = xB[i] */
2139         if (type[k] == GLP_LO || type[k] == GLP_DB ||
2140             type[k] == GLP_FX)
2141         {  /* x[k] has lower bound */
2142            eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
2143            if (bbar[i] < lb[k] - eps)
2144            {  /* and violates it */
2145               coef[k] = -1.0;
2146               cnt++;
2147            }
2148         }
2149         if (type[k] == GLP_UP || type[k] == GLP_DB ||
2150             type[k] == GLP_FX)
2151         {  /* x[k] has upper bound */
2152            eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
2153            if (bbar[i] > ub[k] + eps)
2154            {  /* and violates it */
2155               coef[k] = +1.0;
2156               cnt++;
2157            }
2158         }
2159      }
2160      return cnt;
2161}
2162
2163/***********************************************************************
2164*  set_orig_obj - restore original objective function
2165*
2166*  This routine assigns scaled original objective coefficients to the
2167*  working objective function. */
2168
2169static void set_orig_obj(struct csa *csa)
2170{     int m = csa->m;
2171      int n = csa->n;
2172      double *coef = csa->coef;
2173      double *obj = csa->obj;
2174      double zeta = csa->zeta;
2175      int i, j;
2176      for (i = 1; i <= m; i++)
2177         coef[i] = 0.0;
2178      for (j = 1; j <= n; j++)
2179         coef[m+j] = zeta * obj[j];
2180      return;
2181}
2182
2183/***********************************************************************
2184*  check_stab - check numerical stability of basic solution
2185*
2186*  If the current basic solution is primal feasible (or pseudo feasible
2187*  on phase I) within a tolerance, this routine returns zero, otherwise
2188*  it returns non-zero. */
2189
2190static int check_stab(struct csa *csa, double tol_bnd)
2191{     int m = csa->m;
2192#ifdef GLP_DEBUG
2193      int n = csa->n;
2194#endif
2195      char *type = csa->type;
2196      double *lb = csa->lb;
2197      double *ub = csa->ub;
2198      double *coef = csa->coef;
2199      int *head = csa->head;
2200      int phase = csa->phase;
2201      double *bbar = csa->bbar;
2202      int i, k;
2203      double eps;
2204      /* walk through the list of basic variables */
2205      for (i = 1; i <= m; i++)
2206      {  k = head[i]; /* x[k] = xB[i] */
2207#ifdef GLP_DEBUG
2208         xassert(1 <= k && k <= m+n);
2209#endif
2210         if (phase == 1 && coef[k] < 0.0)
2211         {  /* x[k] must not be greater than its lower bound */
2212#ifdef GLP_DEBUG
2213            xassert(type[k] == GLP_LO || type[k] == GLP_DB ||
2214                    type[k] == GLP_FX);
2215#endif
2216            eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
2217            if (bbar[i] > lb[k] + eps) return 1;
2218         }
2219         else if (phase == 1 && coef[k] > 0.0)
2220         {  /* x[k] must not be less than its upper bound */
2221#ifdef GLP_DEBUG
2222            xassert(type[k] == GLP_UP || type[k] == GLP_DB ||
2223                    type[k] == GLP_FX);
2224#endif
2225            eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
2226            if (bbar[i] < ub[k] - eps) return 1;
2227         }
2228         else
2229         {  /* either phase = 1 and coef[k] = 0, or phase = 2 */
2230            if (type[k] == GLP_LO || type[k] == GLP_DB ||
2231                type[k] == GLP_FX)
2232            {  /* x[k] must not be less than its lower bound */
2233               eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
2234               if (bbar[i] < lb[k] - eps) return 1;
2235            }
2236            if (type[k] == GLP_UP || type[k] == GLP_DB ||
2237                type[k] == GLP_FX)
2238            {  /* x[k] must not be greater then its upper bound */
2239               eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
2240               if (bbar[i] > ub[k] + eps) return 1;
2241            }
2242         }
2243      }
2244      /* basic solution is primal feasible within a tolerance */
2245      return 0;
2246}
2247
2248/***********************************************************************
2249*  check_feas - check primal feasibility of basic solution
2250*
2251*  If the current basic solution is primal feasible within a tolerance,
2252*  this routine returns zero, otherwise it returns non-zero. */
2253
2254static int check_feas(struct csa *csa, double tol_bnd)
2255{     int m = csa->m;
2256#ifdef GLP_DEBUG
2257      int n = csa->n;
2258      char *type = csa->type;
2259#endif
2260      double *lb = csa->lb;
2261      double *ub = csa->ub;
2262      double *coef = csa->coef;
2263      int *head = csa->head;
2264      double *bbar = csa->bbar;
2265      int i, k;
2266      double eps;
2267      xassert(csa->phase == 1);
2268      /* walk through the list of basic variables */
2269      for (i = 1; i <= m; i++)
2270      {  k = head[i]; /* x[k] = xB[i] */
2271#ifdef GLP_DEBUG
2272         xassert(1 <= k && k <= m+n);
2273#endif
2274         if (coef[k] < 0.0)
2275         {  /* check if x[k] still violates its lower bound */
2276#ifdef GLP_DEBUG
2277            xassert(type[k] == GLP_LO || type[k] == GLP_DB ||
2278                    type[k] == GLP_FX);
2279#endif
2280            eps = tol_bnd * (1.0 + kappa * fabs(lb[k]));
2281            if (bbar[i] < lb[k] - eps) return 1;
2282         }
2283         else if (coef[k] > 0.0)
2284         {  /* check if x[k] still violates its upper bound */
2285#ifdef GLP_DEBUG
2286            xassert(type[k] == GLP_UP || type[k] == GLP_DB ||
2287                    type[k] == GLP_FX);
2288#endif
2289            eps = tol_bnd * (1.0 + kappa * fabs(ub[k]));
2290            if (bbar[i] > ub[k] + eps) return 1;
2291         }
2292      }
2293      /* basic solution is primal feasible within a tolerance */
2294      return 0;
2295}
2296
2297/***********************************************************************
2298*  eval_obj - compute original objective function
2299*
2300*  This routine computes the current value of the original objective
2301*  function. */
2302
2303static double eval_obj(struct csa *csa)
2304{     int m = csa->m;
2305      int n = csa->n;
2306      double *obj = csa->obj;
2307      int *head = csa->head;
2308      double *bbar = csa->bbar;
2309      int i, j, k;
2310      double sum;
2311      sum = obj[0];
2312      /* walk through the list of basic variables */
2313      for (i = 1; i <= m; i++)
2314      {  k = head[i]; /* x[k] = xB[i] */
2315#ifdef GLP_DEBUG
2316         xassert(1 <= k && k <= m+n);
2317#endif
2318         if (k > m)
2319            sum += obj[k-m] * bbar[i];
2320      }
2321      /* walk through the list of non-basic variables */
2322      for (j = 1; j <= n; j++)
2323      {  k = head[m+j]; /* x[k] = xN[j] */
2324#ifdef GLP_DEBUG
2325         xassert(1 <= k && k <= m+n);
2326#endif
2327         if (k > m)
2328            sum += obj[k-m] * get_xN(csa, j);
2329      }
2330      return sum;
2331}
2332
2333/***********************************************************************
2334*  display - display the search progress
2335*
2336*  This routine displays some information about the search progress
2337*  that includes:
2338*
2339*  the search phase;
2340*
2341*  the number of simplex iterations performed by the solver;
2342*
2343*  the original objective value;
2344*
2345*  the sum of (scaled) primal infeasibilities;
2346*
2347*  the number of basic fixed variables. */
2348
2349static void display(struct csa *csa, const glp_smcp *parm, int spec)
2350{     int m = csa->m;
2351#ifdef GLP_DEBUG
2352      int n = csa->n;
2353#endif
2354      char *type = csa->type;
2355      double *lb = csa->lb;
2356      double *ub = csa->ub;
2357      int phase = csa->phase;
2358      int *head = csa->head;
2359      double *bbar = csa->bbar;
2360      int i, k, cnt;
2361      double sum;
2362      if (parm->msg_lev < GLP_MSG_ON) goto skip;
2363      if (parm->out_dly > 0 &&
2364         1000.0 * xdifftime(xtime(), csa->tm_beg) < parm->out_dly)
2365         goto skip;
2366      if (csa->it_cnt == csa->it_dpy) goto skip;
2367      if (!spec && csa->it_cnt % parm->out_frq != 0) goto skip;
2368      /* compute the sum of primal infeasibilities and determine the
2369         number of basic fixed variables */
2370      sum = 0.0, cnt = 0;
2371      for (i = 1; i <= m; i++)
2372      {  k = head[i]; /* x[k] = xB[i] */
2373#ifdef GLP_DEBUG
2374         xassert(1 <= k && k <= m+n);
2375#endif
2376         if (type[k] == GLP_LO || type[k] == GLP_DB ||
2377             type[k] == GLP_FX)
2378         {  /* x[k] has lower bound */
2379            if (bbar[i] < lb[k])
2380               sum += (lb[k] - bbar[i]);
2381         }
2382         if (type[k] == GLP_UP || type[k] == GLP_DB ||
2383             type[k] == GLP_FX)
2384         {  /* x[k] has upper bound */
2385            if (bbar[i] > ub[k])
2386               sum += (bbar[i] - ub[k]);
2387         }
2388         if (type[k] == GLP_FX) cnt++;
2389      }
2390      xprintf("%c%6d: obj = %17.9e  infeas = %10.3e (%d)\n",
2391         phase == 1 ? ' ' : '*', csa->it_cnt, eval_obj(csa), sum, cnt);
2392      csa->it_dpy = csa->it_cnt;
2393skip: return;
2394}
2395
2396/***********************************************************************
2397*  store_sol - store basic solution back to the problem object
2398*
2399*  This routine stores basic solution components back to the problem
2400*  object. */
2401
2402static void store_sol(struct csa *csa, glp_prob *lp, int p_stat,
2403      int d_stat, int ray)
2404{     int m = csa->m;
2405      int n = csa->n;
2406      double zeta = csa->zeta;
2407      int *head = csa->head;
2408      char *stat = csa->stat;
2409      double *bbar = csa->bbar;
2410      double *cbar = csa->cbar;
2411      int i, j, k;
2412#ifdef GLP_DEBUG
2413      xassert(lp->m == m);
2414      xassert(lp->n == n);
2415#endif
2416      /* basis factorization */
2417#ifdef GLP_DEBUG
2418      xassert(!lp->valid && lp->bfd == NULL);
2419      xassert(csa->valid && csa->bfd != NULL);
2420#endif
2421      lp->valid = 1, csa->valid = 0;
2422      lp->bfd = csa->bfd, csa->bfd = NULL;
2423      memcpy(&lp->head[1], &head[1], m * sizeof(int));
2424      /* basic solution status */
2425      lp->pbs_stat = p_stat;
2426      lp->dbs_stat = d_stat;
2427      /* objective function value */
2428      lp->obj_val = eval_obj(csa);
2429      /* simplex iteration count */
2430      lp->it_cnt = csa->it_cnt;
2431      /* unbounded ray */
2432      lp->some = ray;
2433      /* basic variables */
2434      for (i = 1; i <= m; i++)
2435      {  k = head[i]; /* x[k] = xB[i] */
2436#ifdef GLP_DEBUG
2437         xassert(1 <= k && k <= m+n);
2438#endif
2439         if (k <= m)
2440         {  GLPROW *row = lp->row[k];
2441            row->stat = GLP_BS;
2442            row->bind = i;
2443            row->prim = bbar[i] / row->rii;
2444            row->dual = 0.0;
2445         }
2446         else
2447         {  GLPCOL *col = lp->col[k-m];
2448            col->stat = GLP_BS;
2449            col->bind = i;
2450            col->prim = bbar[i] * col->sjj;
2451            col->dual = 0.0;
2452         }
2453      }
2454      /* non-basic variables */
2455      for (j = 1; j <= n; j++)
2456      {  k = head[m+j]; /* x[k] = xN[j] */
2457#ifdef GLP_DEBUG
2458         xassert(1 <= k && k <= m+n);
2459#endif
2460         if (k <= m)
2461         {  GLPROW *row = lp->row[k];
2462            row->stat = stat[j];
2463            row->bind = 0;
2464#if 0
2465            row->prim = get_xN(csa, j) / row->rii;
2466#else
2467            switch (stat[j])
2468            {  case GLP_NL:
2469                  row->prim = row->lb; break;
2470               case GLP_NU:
2471                  row->prim = row->ub; break;
2472               case GLP_NF:
2473                  row->prim = 0.0; break;
2474               case GLP_NS:
2475                  row->prim = row->lb; break;
2476               default:
2477                  xassert(stat != stat);
2478            }
2479#endif
2480            row->dual = (cbar[j] * row->rii) / zeta;
2481         }
2482         else
2483         {  GLPCOL *col = lp->col[k-m];
2484            col->stat = stat[j];
2485            col->bind = 0;
2486#if 0
2487            col->prim = get_xN(csa, j) * col->sjj;
2488#else
2489            switch (stat[j])
2490            {  case GLP_NL:
2491                  col->prim = col->lb; break;
2492               case GLP_NU:
2493                  col->prim = col->ub; break;
2494               case GLP_NF:
2495                  col->prim = 0.0; break;
2496               case GLP_NS:
2497                  col->prim = col->lb; break;
2498               default:
2499                  xassert(stat != stat);
2500            }
2501#endif
2502            col->dual = (cbar[j] / col->sjj) / zeta;
2503         }
2504      }
2505      return;
2506}
2507
2508/***********************************************************************
2509*  free_csa - deallocate common storage area
2510*
2511*  This routine frees all the memory allocated to arrays in the common
2512*  storage area (CSA). */
2513
2514static void free_csa(struct csa *csa)
2515{     xfree(csa->type);
2516      xfree(csa->lb);
2517      xfree(csa->ub);
2518      xfree(csa->coef);
2519      xfree(csa->obj);
2520      xfree(csa->A_ptr);
2521      xfree(csa->A_ind);
2522      xfree(csa->A_val);
2523      xfree(csa->head);
2524      xfree(csa->stat);
2525      xfree(csa->N_ptr);
2526      xfree(csa->N_len);
2527      xfree(csa->N_ind);
2528      xfree(csa->N_val);
2529      xfree(csa->bbar);
2530      xfree(csa->cbar);
2531      xfree(csa->refsp);
2532      xfree(csa->gamma);
2533      xfree(csa->tcol_ind);
2534      xfree(csa->tcol_vec);
2535      xfree(csa->trow_ind);
2536      xfree(csa->trow_vec);
2537      xfree(csa->work1);
2538      xfree(csa->work2);
2539      xfree(csa->work3);
2540      xfree(csa->work4);
2541      xfree(csa);
2542      return;
2543}
2544
2545/***********************************************************************
2546*  spx_primal - core LP solver based on the primal simplex method
2547*
2548*  SYNOPSIS
2549*
2550*  #include "glpspx.h"
2551*  int spx_primal(glp_prob *lp, const glp_smcp *parm);
2552*
2553*  DESCRIPTION
2554*
2555*  The routine spx_primal is a core LP solver based on the two-phase
2556*  primal simplex method.
2557*
2558*  RETURNS
2559*
2560*  0  LP instance has been successfully solved.
2561*
2562*  GLP_EITLIM
2563*     Iteration limit has been exhausted.
2564*
2565*  GLP_ETMLIM
2566*     Time limit has been exhausted.
2567*
2568*  GLP_EFAIL
2569*     The solver failed to solve LP instance. */
2570
2571int spx_primal(glp_prob *lp, const glp_smcp *parm)
2572{     struct csa *csa;
2573      int binv_st = 2;
2574      /* status of basis matrix factorization:
2575         0 - invalid; 1 - just computed; 2 - updated */
2576      int bbar_st = 0;
2577      /* status of primal values of basic variables:
2578         0 - invalid; 1 - just computed; 2 - updated */
2579      int cbar_st = 0;
2580      /* status of reduced costs of non-basic variables:
2581         0 - invalid; 1 - just computed; 2 - updated */
2582      int rigorous = 0;
2583      /* rigorous mode flag; this flag is used to enable iterative
2584         refinement on computing pivot rows and columns of the simplex
2585         table */
2586      int check = 0;
2587      int p_stat, d_stat, ret;
2588      /* allocate and initialize the common storage area */
2589      csa = alloc_csa(lp);
2590      init_csa(csa, lp);
2591      if (parm->msg_lev >= GLP_MSG_DBG)
2592         xprintf("Objective scale factor = %g\n", csa->zeta);
2593loop: /* main loop starts here */
2594      /* compute factorization of the basis matrix */
2595      if (binv_st == 0)
2596      {  ret = invert_B(csa);
2597         if (ret != 0)
2598         {  if (parm->msg_lev >= GLP_MSG_ERR)
2599            {  xprintf("Error: unable to factorize the basis matrix (%d"
2600                  ")\n", ret);
2601               xprintf("Sorry, basis recovery procedure not implemented"
2602                  " yet\n");
2603            }
2604            xassert(!lp->valid && lp->bfd == NULL);
2605            lp->bfd = csa->bfd, csa->bfd = NULL;
2606            lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
2607            lp->obj_val = 0.0;
2608            lp->it_cnt = csa->it_cnt;
2609            lp->some = 0;
2610            ret = GLP_EFAIL;
2611            goto done;
2612         }
2613         csa->valid = 1;
2614         binv_st = 1; /* just computed */
2615         /* invalidate basic solution components */
2616         bbar_st = cbar_st = 0;
2617      }
2618      /* compute primal values of basic variables */
2619      if (bbar_st == 0)
2620      {  eval_bbar(csa);
2621         bbar_st = 1; /* just computed */
2622         /* determine the search phase, if not determined yet */
2623         if (csa->phase == 0)
2624         {  if (set_aux_obj(csa, parm->tol_bnd) > 0)
2625            {  /* current basic solution is primal infeasible */
2626               /* start to minimize the sum of infeasibilities */
2627               csa->phase = 1;
2628            }
2629            else
2630            {  /* current basic solution is primal feasible */
2631               /* start to minimize the original objective function */
2632               set_orig_obj(csa);
2633               csa->phase = 2;
2634            }
2635            xassert(check_stab(csa, parm->tol_bnd) == 0);
2636            /* working objective coefficients have been changed, so
2637               invalidate reduced costs */
2638            cbar_st = 0;
2639            display(csa, parm, 1);
2640         }
2641         /* make sure that the current basic solution remains primal
2642            feasible (or pseudo feasible on phase I) */
2643         if (check_stab(csa, parm->tol_bnd))
2644         {  /* there are excessive bound violations due to round-off
2645               errors */
2646            if (parm->msg_lev >= GLP_MSG_ERR)
2647               xprintf("Warning: numerical instability (primal simplex,"
2648                  " phase %s)\n", csa->phase == 1 ? "I" : "II");
2649            /* restart the search */
2650            csa->phase = 0;
2651            binv_st = 0;
2652            rigorous = 5;
2653            goto loop;
2654         }
2655      }
2656      xassert(csa->phase == 1 || csa->phase == 2);
2657      /* on phase I we do not need to wait until the current basic
2658         solution becomes dual feasible; it is sufficient to make sure
2659         that no basic variable violates its bounds */
2660      if (csa->phase == 1 && !check_feas(csa, parm->tol_bnd))
2661      {  /* the current basis is primal feasible; switch to phase II */
2662         csa->phase = 2;
2663         set_orig_obj(csa);
2664         cbar_st = 0;
2665         display(csa, parm, 1);
2666      }
2667      /* compute reduced costs of non-basic variables */
2668      if (cbar_st == 0)
2669      {  eval_cbar(csa);
2670         cbar_st = 1; /* just computed */
2671      }
2672      /* redefine the reference space, if required */
2673      switch (parm->pricing)
2674      {  case GLP_PT_STD:
2675            break;
2676         case GLP_PT_PSE:
2677            if (csa->refct == 0) reset_refsp(csa);
2678            break;
2679         default:
2680            xassert(parm != parm);
2681      }
2682      /* at this point the basis factorization and all basic solution
2683         components are valid */
2684      xassert(binv_st && bbar_st && cbar_st);
2685      /* check accuracy of current basic solution components (only for
2686         debugging) */
2687      if (check)
2688      {  double e_bbar = err_in_bbar(csa);
2689         double e_cbar = err_in_cbar(csa);
2690         double e_gamma =
2691            (parm->pricing == GLP_PT_PSE ? err_in_gamma(csa) : 0.0);
2692         xprintf("e_bbar = %10.3e; e_cbar = %10.3e; e_gamma = %10.3e\n",
2693            e_bbar, e_cbar, e_gamma);
2694         xassert(e_bbar <= 1e-5 && e_cbar <= 1e-5 && e_gamma <= 1e-3);
2695      }
2696      /* check if the iteration limit has been exhausted */
2697      if (parm->it_lim < INT_MAX &&
2698          csa->it_cnt - csa->it_beg >= parm->it_lim)
2699      {  if (bbar_st != 1 || csa->phase == 2 && cbar_st != 1)
2700         {  if (bbar_st != 1) bbar_st = 0;
2701            if (csa->phase == 2 && cbar_st != 1) cbar_st = 0;
2702            goto loop;
2703         }
2704         display(csa, parm, 1);
2705         if (parm->msg_lev >= GLP_MSG_ALL)
2706            xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
2707         switch (csa->phase)
2708         {  case 1:
2709               p_stat = GLP_INFEAS;
2710               set_orig_obj(csa);
2711               eval_cbar(csa);
2712               break;
2713            case 2:
2714               p_stat = GLP_FEAS;
2715               break;
2716            default:
2717               xassert(csa != csa);
2718         }
2719         chuzc(csa, parm->tol_dj);
2720         d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
2721         store_sol(csa, lp, p_stat, d_stat, 0);
2722         ret = GLP_EITLIM;
2723         goto done;
2724      }
2725      /* check if the time limit has been exhausted */
2726      if (parm->tm_lim < INT_MAX &&
2727          1000.0 * xdifftime(xtime(), csa->tm_beg) >= parm->tm_lim)
2728      {  if (bbar_st != 1 || csa->phase == 2 && cbar_st != 1)
2729         {  if (bbar_st != 1) bbar_st = 0;
2730            if (csa->phase == 2 && cbar_st != 1) cbar_st = 0;
2731            goto loop;
2732         }
2733         display(csa, parm, 1);
2734         if (parm->msg_lev >= GLP_MSG_ALL)
2735            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
2736         switch (csa->phase)
2737         {  case 1:
2738               p_stat = GLP_INFEAS;
2739               set_orig_obj(csa);
2740               eval_cbar(csa);
2741               break;
2742            case 2:
2743               p_stat = GLP_FEAS;
2744               break;
2745            default:
2746               xassert(csa != csa);
2747         }
2748         chuzc(csa, parm->tol_dj);
2749         d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
2750         store_sol(csa, lp, p_stat, d_stat, 0);
2751         ret = GLP_ETMLIM;
2752         goto done;
2753      }
2754      /* display the search progress */
2755      display(csa, parm, 0);
2756      /* choose non-basic variable xN[q] */
2757      chuzc(csa, parm->tol_dj);
2758      if (csa->q == 0)
2759      {  if (bbar_st != 1 || cbar_st != 1)
2760         {  if (bbar_st != 1) bbar_st = 0;
2761            if (cbar_st != 1) cbar_st = 0;
2762            goto loop;
2763         }
2764         display(csa, parm, 1);
2765         switch (csa->phase)
2766         {  case 1:
2767               if (parm->msg_lev >= GLP_MSG_ALL)
2768                  xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
2769               p_stat = GLP_NOFEAS;
2770               set_orig_obj(csa);
2771               eval_cbar(csa);
2772               chuzc(csa, parm->tol_dj);
2773               d_stat = (csa->q == 0 ? GLP_FEAS : GLP_INFEAS);
2774               break;
2775            case 2:
2776               if (parm->msg_lev >= GLP_MSG_ALL)
2777                  xprintf("OPTIMAL SOLUTION FOUND\n");
2778               p_stat = d_stat = GLP_FEAS;
2779               break;
2780            default:
2781               xassert(csa != csa);
2782         }
2783         store_sol(csa, lp, p_stat, d_stat, 0);
2784         ret = 0;
2785         goto done;
2786      }
2787      /* compute pivot column of the simplex table */
2788      eval_tcol(csa);
2789      if (rigorous) refine_tcol(csa);
2790      sort_tcol(csa, parm->tol_piv);
2791      /* check accuracy of the reduced cost of xN[q] */
2792      {  double d1 = csa->cbar[csa->q]; /* less accurate */
2793         double d2 = reeval_cost(csa);  /* more accurate */
2794         xassert(d1 != 0.0);
2795         if (fabs(d1 - d2) > 1e-5 * (1.0 + fabs(d2)) ||
2796             !(d1 < 0.0 && d2 < 0.0 || d1 > 0.0 && d2 > 0.0))
2797         {  if (parm->msg_lev >= GLP_MSG_DBG)
2798               xprintf("d1 = %.12g; d2 = %.12g\n", d1, d2);
2799            if (cbar_st != 1 || !rigorous)
2800            {  if (cbar_st != 1) cbar_st = 0;
2801               rigorous = 5;
2802               goto loop;
2803            }
2804         }
2805         /* replace cbar[q] by more accurate value keeping its sign */
2806         if (d1 > 0.0)
2807            csa->cbar[csa->q] = (d2 > 0.0 ? d2 : +DBL_EPSILON);
2808         else
2809            csa->cbar[csa->q] = (d2 < 0.0 ? d2 : -DBL_EPSILON);
2810      }
2811      /* choose basic variable xB[p] */
2812      switch (parm->r_test)
2813      {  case GLP_RT_STD:
2814            chuzr(csa, 0.0);
2815            break;
2816         case GLP_RT_HAR:
2817            chuzr(csa, 0.30 * parm->tol_bnd);
2818            break;
2819         default:
2820            xassert(parm != parm);
2821      }
2822      if (csa->p == 0)
2823      {  if (bbar_st != 1 || cbar_st != 1 || !rigorous)
2824         {  if (bbar_st != 1) bbar_st = 0;
2825            if (cbar_st != 1) cbar_st = 0;
2826            rigorous = 1;
2827            goto loop;
2828         }
2829         display(csa, parm, 1);
2830         switch (csa->phase)
2831         {  case 1:
2832               if (parm->msg_lev >= GLP_MSG_ERR)
2833                  xprintf("Error: unable to choose basic variable on ph"
2834                     "ase I\n");
2835               xassert(!lp->valid && lp->bfd == NULL);
2836               lp->bfd = csa->bfd, csa->bfd = NULL;
2837               lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
2838               lp->obj_val = 0.0;
2839               lp->it_cnt = csa->it_cnt;
2840               lp->some = 0;
2841               ret = GLP_EFAIL;
2842               break;
2843            case 2:
2844               if (parm->msg_lev >= GLP_MSG_ALL)
2845                  xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
2846               store_sol(csa, lp, GLP_FEAS, GLP_NOFEAS,
2847                  csa->head[csa->m+csa->q]);
2848               ret = 0;
2849               break;
2850            default:
2851               xassert(csa != csa);
2852         }
2853         goto done;
2854      }
2855      /* check if the pivot element is acceptable */
2856      if (csa->p > 0)
2857      {  double piv = csa->tcol_vec[csa->p];
2858         double eps = 1e-5 * (1.0 + 0.01 * csa->tcol_max);
2859         if (fabs(piv) < eps)
2860         {  if (parm->msg_lev >= GLP_MSG_DBG)
2861               xprintf("piv = %.12g; eps = %g\n", piv, eps);
2862            if (!rigorous)
2863            {  rigorous = 5;
2864               goto loop;
2865            }
2866         }
2867      }
2868      /* now xN[q] and xB[p] have been chosen anyhow */
2869      /* compute pivot row of the simplex table */
2870      if (csa->p > 0)
2871      {  double *rho = csa->work4;
2872         eval_rho(csa, rho);
2873         if (rigorous) refine_rho(csa, rho);
2874         eval_trow(csa, rho);
2875      }
2876      /* accuracy check based on the pivot element */
2877      if (csa->p > 0)
2878      {  double piv1 = csa->tcol_vec[csa->p]; /* more accurate */
2879         double piv2 = csa->trow_vec[csa->q]; /* less accurate */
2880         xassert(piv1 != 0.0);
2881         if (fabs(piv1 - piv2) > 1e-8 * (1.0 + fabs(piv1)) ||
2882             !(piv1 > 0.0 && piv2 > 0.0 || piv1 < 0.0 && piv2 < 0.0))
2883         {  if (parm->msg_lev >= GLP_MSG_DBG)
2884               xprintf("piv1 = %.12g; piv2 = %.12g\n", piv1, piv2);
2885            if (binv_st != 1 || !rigorous)
2886            {  if (binv_st != 1) binv_st = 0;
2887               rigorous = 5;
2888               goto loop;
2889            }
2890            /* use more accurate version in the pivot row */
2891            if (csa->trow_vec[csa->q] == 0.0)
2892            {  csa->trow_nnz++;
2893               xassert(csa->trow_nnz <= csa->n);
2894               csa->trow_ind[csa->trow_nnz] = csa->q;
2895            }
2896            csa->trow_vec[csa->q] = piv1;
2897         }
2898      }
2899      /* update primal values of basic variables */
2900      update_bbar(csa);
2901      bbar_st = 2; /* updated */
2902      /* update reduced costs of non-basic variables */
2903      if (csa->p > 0)
2904      {  update_cbar(csa);
2905         cbar_st = 2; /* updated */
2906         /* on phase I objective coefficient of xB[p] in the adjacent
2907            basis becomes zero */
2908         if (csa->phase == 1)
2909         {  int k = csa->head[csa->p]; /* x[k] = xB[p] -> xN[q] */
2910            csa->cbar[csa->q] -= csa->coef[k];
2911            csa->coef[k] = 0.0;
2912         }
2913      }
2914      /* update steepest edge coefficients */
2915      if (csa->p > 0)
2916      {  switch (parm->pricing)
2917         {  case GLP_PT_STD:
2918               break;
2919            case GLP_PT_PSE:
2920               if (csa->refct > 0) update_gamma(csa);
2921               break;
2922            default:
2923               xassert(parm != parm);
2924         }
2925      }
2926      /* update factorization of the basis matrix */
2927      if (csa->p > 0)
2928      {  ret = update_B(csa, csa->p, csa->head[csa->m+csa->q]);
2929         if (ret == 0)
2930            binv_st = 2; /* updated */
2931         else
2932         {  csa->valid = 0;
2933            binv_st = 0; /* invalid */
2934         }
2935      }
2936      /* update matrix N */
2937      if (csa->p > 0)
2938      {  del_N_col(csa, csa->q, csa->head[csa->m+csa->q]);
2939         if (csa->type[csa->head[csa->p]] != GLP_FX)
2940            add_N_col(csa, csa->q, csa->head[csa->p]);
2941      }
2942      /* change the basis header */
2943      change_basis(csa);
2944      /* iteration complete */
2945      csa->it_cnt++;
2946      if (rigorous > 0) rigorous--;
2947      goto loop;
2948done: /* deallocate the common storage area */
2949      free_csa(csa);
2950      /* return to the calling program */
2951      return ret;
2952}
2953
2954/* eof */
Note: See TracBrowser for help on using the repository browser.