COIN-OR::LEMON - Graph Library

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