1 /* spxprim.c */
2 
3 /***********************************************************************
4 *  This code is part of GLPK (GNU Linear Programming Kit).
5 *  Copyright (C) 2015-2017 Free Software Foundation, Inc.
6 *  Written by Andrew Makhorin <mao@gnu.org>.
7 *
8 *  GLPK is free software: you can redistribute it and/or modify it
9 *  under the terms of the GNU General Public License as published by
10 *  the Free Software Foundation, either version 3 of the License, or
11 *  (at your option) any later version.
12 *
13 *  GLPK is distributed in the hope that it will be useful, but WITHOUT
14 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16 *  License for more details.
17 *
18 *  You should have received a copy of the GNU General Public License
19 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
20 ***********************************************************************/
21 
22 #if 1 /* 18/VII-2017 */
23 #define SCALE_Z 1
24 #endif
25 
26 #include "env.h"
27 #include "simplex.h"
28 #include "spxat.h"
29 #include "spxnt.h"
30 #include "spxchuzc.h"
31 #include "spxchuzr.h"
32 #include "spxprob.h"
33 
34 #define CHECK_ACCURACY 0
35 /* (for debugging) */
36 
37 struct csa
38 {     /* common storage area */
39       SPXLP *lp;
40       /* LP problem data and its (current) basis; this LP has m rows
41        * and n columns */
42       int dir;
43       /* original optimization direction:
44        * +1 - minimization
45        * -1 - maximization */
46 #if SCALE_Z
47       double fz;
48       /* factor used to scale original objective */
49 #endif
50       double *orig_c; /* double orig_c[1+n]; */
51       /* copy of original objective coefficients */
52       double *orig_l; /* double orig_l[1+n]; */
53       /* copy of original lower bounds */
54       double *orig_u; /* double orig_u[1+n]; */
55       /* copy of original upper bounds */
56       SPXAT *at;
57       /* mxn-matrix A of constraint coefficients, in sparse row-wise
58        * format (NULL if not used) */
59       SPXNT *nt;
60       /* mx(n-m)-matrix N composed of non-basic columns of constraint
61        * matrix A, in sparse row-wise format (NULL if not used) */
62       int phase;
63       /* search phase:
64        * 0 - not determined yet
65        * 1 - searching for primal feasible solution
66        * 2 - searching for optimal solution */
67       double *beta; /* double beta[1+m]; */
68       /* beta[i] is a primal value of basic variable xB[i] */
69       int beta_st;
70       /* status of the vector beta:
71        * 0 - undefined
72        * 1 - just computed
73        * 2 - updated */
74       double *d; /* double d[1+n-m]; */
75       /* d[j] is a reduced cost of non-basic variable xN[j] */
76       int d_st;
77       /* status of the vector d:
78        * 0 - undefined
79        * 1 - just computed
80        * 2 - updated */
81       SPXSE *se;
82       /* projected steepest edge and Devex pricing data block (NULL if
83        * not used) */
84       int num;
85       /* number of eligible non-basic variables */
86       int *list; /* int list[1+n-m]; */
87       /* list[1], ..., list[num] are indices j of eligible non-basic
88        * variables xN[j] */
89       int q;
90       /* xN[q] is a non-basic variable chosen to enter the basis */
91 #if 0 /* 11/VI-2017 */
92       double *tcol; /* double tcol[1+m]; */
93 #else
94       FVS tcol; /* FVS tcol[1:m]; */
95 #endif
96       /* q-th (pivot) column of the simplex table */
97 #if 1 /* 23/VI-2017 */
98       SPXBP *bp; /* SPXBP bp[1+2*m+1]; */
99       /* penalty function break points */
100 #endif
101       int p;
102       /* xB[p] is a basic variable chosen to leave the basis;
103        * p = 0 means that no basic variable reaches its bound;
104        * p < 0 means that non-basic variable xN[q] reaches its opposite
105        * bound before any basic variable */
106       int p_flag;
107       /* if this flag is set, the active bound of xB[p] in the adjacent
108        * basis should be set to the upper bound */
109 #if 0 /* 11/VI-2017 */
110       double *trow; /* double trow[1+n-m]; */
111 #else
112       FVS trow; /* FVS trow[1:n-m]; */
113 #endif
114       /* p-th (pivot) row of the simplex table */
115 #if 0 /* 09/VII-2017 */
116       double *work; /* double work[1+m]; */
117       /* working array */
118 #else
119       FVS work; /* FVS work[1:m]; */
120       /* working vector */
121 #endif
122       int p_stat, d_stat;
123       /* primal and dual solution statuses */
124       /*--------------------------------------------------------------*/
125       /* control parameters (see struct glp_smcp) */
126       int msg_lev;
127       /* message level */
128 #if 0 /* 23/VI-2017 */
129       int harris;
130       /* ratio test technique:
131        * 0 - textbook ratio test
132        * 1 - Harris' two pass ratio test */
133 #else
134       int r_test;
135       /* ratio test technique:
136        * GLP_RT_STD  - textbook ratio test
137        * GLP_RT_HAR  - Harris' two pass ratio test
138        * GLP_RT_FLIP - long-step ratio test (only for phase I) */
139 #endif
140       double tol_bnd, tol_bnd1;
141       /* primal feasibility tolerances */
142       double tol_dj, tol_dj1;
143       /* dual feasibility tolerances */
144       double tol_piv;
145       /* pivot tolerance */
146       int it_lim;
147       /* iteration limit */
148       int tm_lim;
149       /* time limit, milliseconds */
150       int out_frq;
151 #if 0 /* 15/VII-2017 */
152       /* display output frequency, iterations */
153 #else
154       /* display output frequency, milliseconds */
155 #endif
156       int out_dly;
157       /* display output delay, milliseconds */
158       /*--------------------------------------------------------------*/
159       /* working parameters */
160       double tm_beg;
161       /* time value at the beginning of the search */
162       int it_beg;
163       /* simplex iteration count at the beginning of the search */
164       int it_cnt;
165       /* simplex iteration count; it increases by one every time the
166        * basis changes (including the case when a non-basic variable
167        * jumps to its opposite bound) */
168       int it_dpy;
169       /* simplex iteration count at most recent display output */
170 #if 1 /* 15/VII-2017 */
171       double tm_dpy;
172       /* time value at most recent display output */
173 #endif
174       int inv_cnt;
175       /* basis factorization count since most recent display output */
176 #if 1 /* 01/VII-2017 */
177       int degen;
178       /* count of successive degenerate iterations; this count is used
179        * to detect stalling */
180 #endif
181 #if 1 /* 23/VI-2017 */
182       int ns_cnt, ls_cnt;
183       /* normal and long-step iteration counts */
184 #endif
185 };
186 
187 /***********************************************************************
188 *  set_penalty - set penalty function coefficients
189 *
190 *  This routine sets up objective coefficients of the penalty function,
191 *  which is the sum of primal infeasibilities, as follows:
192 *
193 *     if beta[i] < l[k] - eps1, set c[k] = -1,
194 *
195 *     if beta[i] > u[k] + eps2, set c[k] = +1,
196 *
197 *     otherwise, set c[k] = 0,
198 *
199 *  where beta[i] is current value of basic variable xB[i] = x[k], l[k]
200 *  and u[k] are original bounds of x[k], and
201 *
202 *     eps1 = tol + tol1 * |l[k]|,
203 *
204 *     eps2 = tol + tol1 * |u[k]|.
205 *
206 *  The routine returns the number of non-zero objective coefficients,
207 *  which is the number of basic variables violating their bounds. Thus,
208 *  if the value returned is zero, the current basis is primal feasible
209 *  within the specified tolerances. */
210 
set_penalty(struct csa * csa,double tol,double tol1)211 static int set_penalty(struct csa *csa, double tol, double tol1)
212 {     SPXLP *lp = csa->lp;
213       int m = lp->m;
214       int n = lp->n;
215       double *c = lp->c;
216       double *l = lp->l;
217       double *u = lp->u;
218       int *head = lp->head;
219       double *beta = csa->beta;
220       int i, k, count = 0;
221       double t, eps;
222       /* reset objective coefficients */
223       for (k = 0; k <= n; k++)
224          c[k] = 0.0;
225       /* walk thru the list of basic variables */
226       for (i = 1; i <= m; i++)
227       {  k = head[i]; /* x[k] = xB[i] */
228          /* check lower bound */
229          if ((t = l[k]) != -DBL_MAX)
230          {  eps = tol + tol1 * (t >= 0.0 ? +t : -t);
231             if (beta[i] < t - eps)
232             {  /* lower bound is violated */
233                c[k] = -1.0, count++;
234             }
235          }
236          /* check upper bound */
237          if ((t = u[k]) != +DBL_MAX)
238          {  eps = tol + tol1 * (t >= 0.0 ? +t : -t);
239             if (beta[i] > t + eps)
240             {  /* upper bound is violated */
241                c[k] = +1.0, count++;
242             }
243          }
244       }
245       return count;
246 }
247 
248 /***********************************************************************
249 *  check_feas - check primal feasibility of basic solution
250 *
251 *  This routine checks if the specified values of all basic variables
252 *  beta = (beta[i]) are within their bounds.
253 *
254 *  Let l[k] and u[k] be original bounds of basic variable xB[i] = x[k].
255 *  The actual bounds of x[k] are determined as follows:
256 *
257 *  1) if phase = 1 and c[k] < 0, x[k] violates its lower bound, so its
258 *     actual bounds are artificial: -inf < x[k] <= l[k];
259 *
260 *  2) if phase = 1 and c[k] > 0, x[k] violates its upper bound, so its
261 *     actual bounds are artificial: u[k] <= x[k] < +inf;
262 *
263 *  3) in all other cases (if phase = 1 and c[k] = 0, or if phase = 2)
264 *     actual bounds are original: l[k] <= x[k] <= u[k].
265 *
266 *  The parameters tol and tol1 are bound violation tolerances. The
267 *  actual bounds l'[k] and u'[k] are considered as non-violated within
268 *  the specified tolerance if
269 *
270 *     l'[k] - eps1 <= beta[i] <= u'[k] + eps2,
271 *
272 *  where eps1 = tol + tol1 * |l'[k]|, eps2 = tol + tol1 * |u'[k]|.
273 *
274 *  The routine returns one of the following codes:
275 *
276 *  0 - solution is feasible (no actual bounds are violated);
277 *
278 *  1 - solution is infeasible, however, only artificial bounds are
279 *      violated (this is possible only if phase = 1);
280 *
281 *  2 - solution is infeasible and at least one original bound is
282 *      violated. */
283 
check_feas(struct csa * csa,int phase,double tol,double tol1)284 static int check_feas(struct csa *csa, int phase, double tol, double
285       tol1)
286 {     SPXLP *lp = csa->lp;
287       int m = lp->m;
288       double *c = lp->c;
289       double *l = lp->l;
290       double *u = lp->u;
291       int *head = lp->head;
292       double *beta = csa->beta;
293       int i, k, orig, ret = 0;
294       double lk, uk, eps;
295       xassert(phase == 1 || phase == 2);
296       /* walk thru the list of basic variables */
297       for (i = 1; i <= m; i++)
298       {  k = head[i]; /* x[k] = xB[i] */
299          /* determine actual bounds of x[k] */
300          if (phase == 1 && c[k] < 0.0)
301          {  /* -inf < x[k] <= l[k] */
302             lk = -DBL_MAX, uk = l[k];
303             orig = 0; /* artificial bounds */
304          }
305          else if (phase == 1 && c[k] > 0.0)
306          {  /* u[k] <= x[k] < +inf */
307             lk = u[k], uk = +DBL_MAX;
308             orig = 0; /* artificial bounds */
309          }
310          else
311          {  /* l[k] <= x[k] <= u[k] */
312             lk = l[k], uk = u[k];
313             orig = 1; /* original bounds */
314          }
315          /* check actual lower bound */
316          if (lk != -DBL_MAX)
317          {  eps = tol + tol1 * (lk >= 0.0 ? +lk : -lk);
318             if (beta[i] < lk - eps)
319             {  /* actual lower bound is violated */
320                if (orig)
321                {  ret = 2;
322                   break;
323                }
324                ret = 1;
325             }
326          }
327          /* check actual upper bound */
328          if (uk != +DBL_MAX)
329          {  eps = tol + tol1 * (uk >= 0.0 ? +uk : -uk);
330             if (beta[i] > uk + eps)
331             {  /* actual upper bound is violated */
332                if (orig)
333                {  ret = 2;
334                   break;
335                }
336                ret = 1;
337             }
338          }
339       }
340       return ret;
341 }
342 
343 /***********************************************************************
344 *  adjust_penalty - adjust penalty function coefficients
345 *
346 *  On searching for primal feasible solution it may happen that some
347 *  basic variable xB[i] = x[k] has non-zero objective coefficient c[k]
348 *  indicating that xB[i] violates its lower (if c[k] < 0) or upper (if
349 *  c[k] > 0) original bound, but due to primal degenarcy the violation
350 *  is close to zero.
351 *
352 *  This routine identifies such basic variables and sets objective
353 *  coefficients at these variables to zero that allows avoiding zero-
354 *  step simplex iterations.
355 *
356 *  The parameters tol and tol1 are bound violation tolerances. The
357 *  original bounds l[k] and u[k] are considered as non-violated within
358 *  the specified tolerance if
359 *
360 *     l[k] - eps1 <= beta[i] <= u[k] + eps2,
361 *
362 *  where beta[i] is value of basic variable xB[i] = x[k] in the current
363 *  basis, eps1 = tol + tol1 * |l[k]|, eps2 = tol + tol1 * |u[k]|.
364 *
365 *  The routine returns the number of objective coefficients which were
366 *  set to zero. */
367 
368 #if 0
369 static int adjust_penalty(struct csa *csa, double tol, double tol1)
370 {     SPXLP *lp = csa->lp;
371       int m = lp->m;
372       double *c = lp->c;
373       double *l = lp->l;
374       double *u = lp->u;
375       int *head = lp->head;
376       double *beta = csa->beta;
377       int i, k, count = 0;
378       double t, eps;
379       xassert(csa->phase == 1);
380       /* walk thru the list of basic variables */
381       for (i = 1; i <= m; i++)
382       {  k = head[i]; /* x[k] = xB[i] */
383          if (c[k] < 0.0)
384          {  /* x[k] violates its original lower bound l[k] */
385             xassert((t = l[k]) != -DBL_MAX);
386             eps = tol + tol1 * (t >= 0.0 ? +t : -t);
387             if (beta[i] >= t - eps)
388             {  /* however, violation is close to zero */
389                c[k] = 0.0, count++;
390             }
391          }
392          else if (c[k] > 0.0)
393          {  /* x[k] violates its original upper bound u[k] */
394             xassert((t = u[k]) != +DBL_MAX);
395             eps = tol + tol1 * (t >= 0.0 ? +t : -t);
396             if (beta[i] <= t + eps)
397             {  /* however, violation is close to zero */
398                c[k] = 0.0, count++;
399             }
400          }
401       }
402       return count;
403 }
404 #else
adjust_penalty(struct csa * csa,int num,const int ind[],double tol,double tol1)405 static int adjust_penalty(struct csa *csa, int num, const int
406       ind[/*1+num*/], double tol, double tol1)
407 {     SPXLP *lp = csa->lp;
408       int m = lp->m;
409       double *c = lp->c;
410       double *l = lp->l;
411       double *u = lp->u;
412       int *head = lp->head;
413       double *beta = csa->beta;
414       int i, k, t, cnt = 0;
415       double lk, uk, eps;
416       xassert(csa->phase == 1);
417       /* walk thru the specified list of basic variables */
418       for (t = 1; t <= num; t++)
419       {  i = ind[t];
420          xassert(1 <= i && i <= m);
421          k = head[i]; /* x[k] = xB[i] */
422          if (c[k] < 0.0)
423          {  /* x[k] violates its original lower bound */
424             lk = l[k];
425             xassert(lk != -DBL_MAX);
426             eps = tol + tol1 * (lk >= 0.0 ? +lk : -lk);
427             if (beta[i] >= lk - eps)
428             {  /* however, violation is close to zero */
429                c[k] = 0.0, cnt++;
430             }
431          }
432          else if (c[k] > 0.0)
433          {  /* x[k] violates its original upper bound */
434             uk = u[k];
435             xassert(uk != +DBL_MAX);
436             eps = tol + tol1 * (uk >= 0.0 ? +uk : -uk);
437             if (beta[i] <= uk + eps)
438             {  /* however, violation is close to zero */
439                c[k] = 0.0, cnt++;
440             }
441          }
442       }
443       return cnt;
444 }
445 #endif
446 
447 #if CHECK_ACCURACY
448 /***********************************************************************
449 *  err_in_vec - compute maximal relative error between two vectors
450 *
451 *  This routine computes and returns maximal relative error between
452 *  n-vectors x and y:
453 *
454 *     err_max = max |x[i] - y[i]| / (1 + |x[i]|).
455 *
456 *  NOTE: This routine is intended only for debugginig purposes. */
457 
err_in_vec(int n,const double x[],const double y[])458 static double err_in_vec(int n, const double x[], const double y[])
459 {     int i;
460       double err, err_max;
461       err_max = 0.0;
462       for (i = 1; i <= n; i++)
463       {  err = fabs(x[i] - y[i]) / (1.0 + fabs(x[i]));
464          if (err_max < err)
465             err_max = err;
466       }
467       return err_max;
468 }
469 #endif
470 
471 #if CHECK_ACCURACY
472 /***********************************************************************
473 *  err_in_beta - compute maximal relative error in vector beta
474 *
475 *  This routine computes and returns maximal relative error in vector
476 *  of values of basic variables beta = (beta[i]).
477 *
478 *  NOTE: This routine is intended only for debugginig purposes. */
479 
err_in_beta(struct csa * csa)480 static double err_in_beta(struct csa *csa)
481 {     SPXLP *lp = csa->lp;
482       int m = lp->m;
483       double err, *beta;
484       beta = talloc(1+m, double);
485       spx_eval_beta(lp, beta);
486       err = err_in_vec(m, beta, csa->beta);
487       tfree(beta);
488       return err;
489 }
490 #endif
491 
492 #if CHECK_ACCURACY
493 /***********************************************************************
494 *  err_in_d - compute maximal relative error in vector d
495 *
496 *  This routine computes and returns maximal relative error in vector
497 *  of reduced costs of non-basic variables d = (d[j]).
498 *
499 *  NOTE: This routine is intended only for debugginig purposes. */
500 
err_in_d(struct csa * csa)501 static double err_in_d(struct csa *csa)
502 {     SPXLP *lp = csa->lp;
503       int m = lp->m;
504       int n = lp->n;
505       int j;
506       double err, *pi, *d;
507       pi = talloc(1+m, double);
508       d = talloc(1+n-m, double);
509       spx_eval_pi(lp, pi);
510       for (j = 1; j <= n-m; j++)
511          d[j] = spx_eval_dj(lp, pi, j);
512       err = err_in_vec(n-m, d, csa->d);
513       tfree(pi);
514       tfree(d);
515       return err;
516 }
517 #endif
518 
519 #if CHECK_ACCURACY
520 /***********************************************************************
521 *  err_in_gamma - compute maximal relative error in vector gamma
522 *
523 *  This routine computes and returns maximal relative error in vector
524 *  of projected steepest edge weights gamma = (gamma[j]).
525 *
526 *  NOTE: This routine is intended only for debugginig purposes. */
527 
err_in_gamma(struct csa * csa)528 static double err_in_gamma(struct csa *csa)
529 {     SPXLP *lp = csa->lp;
530       int m = lp->m;
531       int n = lp->n;
532       SPXSE *se = csa->se;
533       int j;
534       double err, *gamma;
535       xassert(se != NULL);
536       gamma = talloc(1+n-m, double);
537       for (j = 1; j <= n-m; j++)
538          gamma[j] = spx_eval_gamma_j(lp, se, j);
539       err = err_in_vec(n-m, gamma, se->gamma);
540       tfree(gamma);
541       return err;
542 }
543 #endif
544 
545 #if CHECK_ACCURACY
546 /***********************************************************************
547 *  check_accuracy - check accuracy of basic solution components
548 *
549 *  This routine checks accuracy of current basic solution components.
550 *
551 *  NOTE: This routine is intended only for debugginig purposes. */
552 
check_accuracy(struct csa * csa)553 static void check_accuracy(struct csa *csa)
554 {     double e_beta, e_d, e_gamma;
555       e_beta = err_in_beta(csa);
556       e_d = err_in_d(csa);
557       if (csa->se == NULL)
558          e_gamma = 0.;
559       else
560          e_gamma = err_in_gamma(csa);
561       xprintf("e_beta = %10.3e; e_d = %10.3e; e_gamma = %10.3e\n",
562          e_beta, e_d, e_gamma);
563       xassert(e_beta <= 1e-5 && e_d <= 1e-5 && e_gamma <= 1e-3);
564       return;
565 }
566 #endif
567 
568 /***********************************************************************
569 *  choose_pivot - choose xN[q] and xB[p]
570 *
571 *  Given the list of eligible non-basic variables this routine first
572 *  chooses non-basic variable xN[q]. This choice is always possible,
573 *  because the list is assumed to be non-empty. Then the routine
574 *  computes q-th column T[*,q] of the simplex table T[i,j] and chooses
575 *  basic variable xB[p]. If the pivot T[p,q] is small in magnitude,
576 *  the routine attempts to choose another xN[q] and xB[p] in order to
577 *  avoid badly conditioned adjacent bases. */
578 
579 #if 1 /* 17/III-2016 */
580 #define MIN_RATIO 0.0001
581 
choose_pivot(struct csa * csa)582 static int choose_pivot(struct csa *csa)
583 {     SPXLP *lp = csa->lp;
584       int m = lp->m;
585       int n = lp->n;
586       double *beta = csa->beta;
587       double *d = csa->d;
588       SPXSE *se = csa->se;
589       int *list = csa->list;
590 #if 0 /* 09/VII-2017 */
591       double *tcol = csa->work;
592 #else
593       double *tcol = csa->work.vec;
594 #endif
595       double tol_piv = csa->tol_piv;
596       int try, nnn, /*i,*/ p, p_flag, q, t;
597       double big, /*temp,*/ best_ratio;
598 #if 1 /* 23/VI-2017 */
599       double *c = lp->c;
600       int *head = lp->head;
601       SPXBP *bp = csa->bp;
602       int nbp, t_best, ret, k;
603       double dz_best;
604 #endif
605       xassert(csa->beta_st);
606       xassert(csa->d_st);
607 more: /* initial number of eligible non-basic variables */
608       nnn = csa->num;
609       /* nothing has been chosen so far */
610       csa->q = 0;
611       best_ratio = 0.0;
612 #if 0 /* 23/VI-2017 */
613       try = 0;
614 #else
615       try = ret = 0;
616 #endif
617 try:  /* choose non-basic variable xN[q] */
618       xassert(nnn > 0);
619       try++;
620       if (se == NULL)
621       {  /* Dantzig's rule */
622          q = spx_chuzc_std(lp, d, nnn, list);
623       }
624       else
625       {  /* projected steepest edge */
626          q = spx_chuzc_pse(lp, se, d, nnn, list);
627       }
628       xassert(1 <= q && q <= n-m);
629       /* compute q-th column of the simplex table */
630       spx_eval_tcol(lp, q, tcol);
631 #if 0
632       /* big := max(1, |tcol[1]|, ..., |tcol[m]|) */
633       big = 1.0;
634       for (i = 1; i <= m; i++)
635       {  temp = tcol[i];
636          if (temp < 0.0)
637             temp = - temp;
638          if (big < temp)
639             big = temp;
640       }
641 #else
642       /* this still puzzles me */
643       big = 1.0;
644 #endif
645       /* choose basic variable xB[p] */
646 #if 1 /* 23/VI-2017 */
647       if (csa->phase == 1 && csa->r_test == GLP_RT_FLIP && try <= 2)
648       {  /* long-step ratio test */
649          int t, num, num1;
650          double slope, teta_lim;
651          /* determine penalty function break points */
652          nbp = spx_ls_eval_bp(lp, beta, q, d[q], tcol, tol_piv, bp);
653          if (nbp < 2)
654             goto skip;
655          /* set initial slope */
656          slope = - fabs(d[q]);
657          /* estimate initial teta_lim */
658          teta_lim = DBL_MAX;
659          for (t = 1; t <= nbp; t++)
660          {  if (teta_lim > bp[t].teta)
661                teta_lim = bp[t].teta;
662          }
663          xassert(teta_lim >= 0.0);
664          if (teta_lim < 1e-3)
665             teta_lim = 1e-3;
666          /* nothing has been chosen so far */
667          t_best = 0, dz_best = 0.0, num = 0;
668          /* choose appropriate break point */
669          while (num < nbp)
670          {  /* select and process a new portion of break points */
671             num1 = spx_ls_select_bp(lp, tcol, nbp, bp, num, &slope,
672                teta_lim);
673             for (t = num+1; t <= num1; t++)
674             {  int i = (bp[t].i >= 0 ? bp[t].i : -bp[t].i);
675                xassert(0 <= i && i <= m);
676                if (i == 0 || fabs(tcol[i]) / big >= MIN_RATIO)
677                {  if (dz_best > bp[t].dz)
678                      t_best = t, dz_best = bp[t].dz;
679                }
680 #if 0
681                if (i == 0)
682                {  /* do not consider further break points beyond this
683                    * point, where xN[q] reaches its opposite bound;
684                    * in principle (see spx_ls_eval_bp), this break
685                    * point should be the last one, however, due to
686                    * round-off errors there may be other break points
687                    * with the same teta beyond this one */
688                   slope = +1.0;
689                }
690 #endif
691             }
692             if (slope > 0.0)
693             {  /* penalty function starts increasing */
694                break;
695             }
696             /* penalty function continues decreasing */
697             num = num1;
698             teta_lim += teta_lim;
699          }
700          if (dz_best == 0.0)
701             goto skip;
702          /* the choice has been made */
703          xassert(1 <= t_best && t_best <= num1);
704          if (t_best == 1)
705          {  /* the very first break point was chosen; it is reasonable
706              * to use the short-step ratio test */
707             goto skip;
708          }
709          csa->q = q;
710          memcpy(&csa->tcol.vec[1], &tcol[1], m * sizeof(double));
711          fvs_gather_vec(&csa->tcol, DBL_EPSILON);
712          if (bp[t_best].i == 0)
713          {  /* xN[q] goes to its opposite bound */
714             csa->p = -1;
715             csa->p_flag = 0;
716             best_ratio = 1.0;
717          }
718          else if (bp[t_best].i > 0)
719          {  /* xB[p] leaves the basis and goes to its lower bound */
720             csa->p = + bp[t_best].i;
721             xassert(1 <= csa->p && csa->p <= m);
722             csa->p_flag = 0;
723             best_ratio = fabs(tcol[csa->p]) / big;
724          }
725          else
726          {  /* xB[p] leaves the basis and goes to its upper bound */
727             csa->p = - bp[t_best].i;
728             xassert(1 <= csa->p && csa->p <= m);
729             csa->p_flag = 1;
730             best_ratio = fabs(tcol[csa->p]) / big;
731          }
732 #if 0
733          xprintf("num1 = %d; t_best = %d; dz = %g\n", num1, t_best,
734             bp[t_best].dz);
735 #endif
736          ret = 1;
737          goto done;
738 skip:    ;
739       }
740 #endif
741 #if 0 /* 23/VI-2017 */
742       if (!csa->harris)
743 #else
744       if (csa->r_test == GLP_RT_STD)
745 #endif
746       {  /* textbook ratio test */
747          p = spx_chuzr_std(lp, csa->phase, beta, q,
748             d[q] < 0.0 ? +1. : -1., tcol, &p_flag, tol_piv,
749             .30 * csa->tol_bnd, .30 * csa->tol_bnd1);
750       }
751       else
752       {  /* Harris' two-pass ratio test */
753          p = spx_chuzr_harris(lp, csa->phase, beta, q,
754             d[q] < 0.0 ? +1. : -1., tcol, &p_flag , tol_piv,
755             .50 * csa->tol_bnd, .50 * csa->tol_bnd1);
756       }
757       if (p <= 0)
758       {  /* primal unboundedness or special case */
759          csa->q = q;
760 #if 0 /* 11/VI-2017 */
761          memcpy(&csa->tcol[1], &tcol[1], m * sizeof(double));
762 #else
763          memcpy(&csa->tcol.vec[1], &tcol[1], m * sizeof(double));
764          fvs_gather_vec(&csa->tcol, DBL_EPSILON);
765 #endif
766          csa->p = p;
767          csa->p_flag = p_flag;
768          best_ratio = 1.0;
769          goto done;
770       }
771       /* either keep previous choice or accept new choice depending on
772        * which one is better */
773       if (best_ratio < fabs(tcol[p]) / big)
774       {  csa->q = q;
775 #if 0 /* 11/VI-2017 */
776          memcpy(&csa->tcol[1], &tcol[1], m * sizeof(double));
777 #else
778          memcpy(&csa->tcol.vec[1], &tcol[1], m * sizeof(double));
779          fvs_gather_vec(&csa->tcol, DBL_EPSILON);
780 #endif
781          csa->p = p;
782          csa->p_flag = p_flag;
783          best_ratio = fabs(tcol[p]) / big;
784       }
785       /* check if the current choice is acceptable */
786       if (best_ratio >= MIN_RATIO || nnn == 1 || try == 5)
787          goto done;
788       /* try to choose other xN[q] and xB[p] */
789       /* find xN[q] in the list */
790       for (t = 1; t <= nnn; t++)
791          if (list[t] == q) break;
792       xassert(t <= nnn);
793       /* move xN[q] to the end of the list */
794       list[t] = list[nnn], list[nnn] = q;
795       /* and exclude it from consideration */
796       nnn--;
797       /* repeat the choice */
798       goto try;
799 done: /* the choice has been made */
800 #if 1 /* FIXME: currently just to avoid badly conditioned basis */
801       if (best_ratio < .001 * MIN_RATIO)
802       {  /* looks like this helps */
803          if (bfd_get_count(lp->bfd) > 0)
804             return -1;
805          /* didn't help; last chance to improve the choice */
806          if (tol_piv == csa->tol_piv)
807          {  tol_piv *= 1000.;
808             goto more;
809          }
810       }
811 #endif
812 #if 0 /* 23/VI-2017 */
813       return 0;
814 #else /* FIXME */
815       if (ret)
816       {  /* invalidate dual basic solution components */
817          csa->d_st = 0;
818          /* change penalty function coefficients at basic variables for
819           * all break points preceding the chosen one */
820          for (t = 1; t < t_best; t++)
821          {  int i = (bp[t].i >= 0 ? bp[t].i : -bp[t].i);
822             xassert(0 <= i && i <= m);
823             if (i == 0)
824             {  /* xN[q] crosses its opposite bound */
825                xassert(1 <= csa->q && csa->q <= n-m);
826                k = head[m+csa->q];
827             }
828             else
829             {  /* xB[i] crosses its (lower or upper) bound */
830                k = head[i]; /* x[k] = xB[i] */
831             }
832             c[k] += bp[t].dc;
833             xassert(c[k] == 0.0 || c[k] == +1.0 || c[k] == -1.0);
834          }
835       }
836       return ret;
837 #endif
838 }
839 #endif
840 
841 /***********************************************************************
842 *  play_bounds - play bounds of primal variables
843 *
844 *  This routine is called after the primal values of basic variables
845 *  beta[i] were updated and the basis was changed to the adjacent one.
846 *
847 *  It is assumed that before updating all the primal values beta[i]
848 *  were strongly feasible, so in the adjacent basis beta[i] remain
849 *  feasible within a tolerance, i.e. if some beta[i] violates its lower
850 *  or upper bound, the violation is insignificant.
851 *
852 *  If some beta[i] violates its lower or upper bound, this routine
853 *  changes (perturbs) the bound to remove such violation, i.e. to make
854 *  all beta[i] strongly feasible. Otherwise, if beta[i] has a feasible
855 *  value, this routine attempts to reduce (or remove) perturbation of
856 *  corresponding lower/upper bound keeping strong feasibility. */
857 
858 /* FIXME: what to do if l[k] = u[k]? */
859 
860 /* FIXME: reduce/remove perturbation if x[k] becomes non-basic? */
861 
play_bounds(struct csa * csa,int all)862 static void play_bounds(struct csa *csa, int all)
863 {     SPXLP *lp = csa->lp;
864       int m = lp->m;
865       double *c = lp->c;
866       double *l = lp->l;
867       double *u = lp->u;
868       int *head = lp->head;
869       double *orig_l = csa->orig_l;
870       double *orig_u = csa->orig_u;
871       double *beta = csa->beta;
872 #if 0 /* 11/VI-2017 */
873       const double *tcol = csa->tcol; /* was used to update beta */
874 #else
875       const double *tcol = csa->tcol.vec;
876 #endif
877       int i, k;
878       xassert(csa->phase == 1 || csa->phase == 2);
879       /* primal values beta = (beta[i]) should be valid */
880       xassert(csa->beta_st);
881       /* walk thru the list of basic variables xB = (xB[i]) */
882       for (i = 1; i <= m; i++)
883       {  if (all || tcol[i] != 0.0)
884          {  /* beta[i] has changed in the adjacent basis */
885             k = head[i]; /* x[k] = xB[i] */
886             if (csa->phase == 1 && c[k] < 0.0)
887             {  /* -inf < xB[i] <= lB[i] (artificial bounds) */
888                if (beta[i] < l[k] - 1e-9)
889                   continue;
890                /* restore actual bounds */
891                c[k] = 0.0;
892                csa->d_st = 0; /* since c[k] = cB[i] has changed */
893             }
894             if (csa->phase == 1 && c[k] > 0.0)
895             {  /* uB[i] <= xB[i] < +inf (artificial bounds) */
896                if (beta[i] > u[k] + 1e-9)
897                   continue;
898                /* restore actual bounds */
899                c[k] = 0.0;
900                csa->d_st = 0; /* since c[k] = cB[i] has changed */
901             }
902             /* lB[i] <= xB[i] <= uB[i] */
903             if (csa->phase == 1)
904                xassert(c[k] == 0.0);
905             if (l[k] != -DBL_MAX)
906             {  /* xB[i] has lower bound */
907                if (beta[i] < l[k])
908                {  /* strong feasibility means xB[i] >= lB[i] */
909 #if 0 /* 11/VI-2017 */
910                   l[k] = beta[i];
911 #else
912                   l[k] = beta[i] - 1e-9;
913 #endif
914                }
915                else if (l[k] < orig_l[k])
916                {  /* remove/reduce perturbation of lB[i] */
917                   if (beta[i] >= orig_l[k])
918                      l[k] = orig_l[k];
919                   else
920                      l[k] = beta[i];
921                }
922             }
923             if (u[k] != +DBL_MAX)
924             {  /* xB[i] has upper bound */
925                if (beta[i] > u[k])
926                {  /* strong feasibility means xB[i] <= uB[i] */
927 #if 0 /* 11/VI-2017 */
928                   u[k] = beta[i];
929 #else
930                   u[k] = beta[i] + 1e-9;
931 #endif
932                }
933                else if (u[k] > orig_u[k])
934                {  /* remove/reduce perturbation of uB[i] */
935                   if (beta[i] <= orig_u[k])
936                      u[k] = orig_u[k];
937                   else
938                      u[k] = beta[i];
939                }
940             }
941          }
942       }
943       return;
944 }
945 
remove_perturb(struct csa * csa)946 static void remove_perturb(struct csa *csa)
947 {     /* remove perturbation */
948       SPXLP *lp = csa->lp;
949       int m = lp->m;
950       int n = lp->n;
951       double *l = lp->l;
952       double *u = lp->u;
953       int *head = lp->head;
954       char *flag = lp->flag;
955       double *orig_l = csa->orig_l;
956       double *orig_u = csa->orig_u;
957       int j, k;
958       /* restore original bounds of variables */
959       memcpy(l, orig_l, (1+n) * sizeof(double));
960       memcpy(u, orig_u, (1+n) * sizeof(double));
961       /* adjust flags of fixed non-basic variables, because in the
962        * perturbed problem such variables might be changed to double-
963        * bounded type */
964       for (j = 1; j <= n-m; j++)
965       {  k = head[m+j]; /* x[k] = xN[j] */
966          if (l[k] == u[k])
967             flag[j] = 0;
968       }
969       /* removing perturbation changes primal solution components */
970       csa->phase = csa->beta_st = 0;
971 #if 1
972       if (csa->msg_lev >= GLP_MSG_ALL)
973          xprintf("Removing LP perturbation [%d]...\n",
974             csa->it_cnt);
975 #endif
976       return;
977 }
978 
979 /***********************************************************************
980 *  sum_infeas - compute sum of primal infeasibilities
981 *
982 *  This routine compute the sum of primal infeasibilities, which is the
983 *  current penalty function value. */
984 
sum_infeas(SPXLP * lp,const double beta[])985 static double sum_infeas(SPXLP *lp, const double beta[/*1+m*/])
986 {     int m = lp->m;
987       double *l = lp->l;
988       double *u = lp->u;
989       int *head = lp->head;
990       int i, k;
991       double sum = 0.0;
992       for (i = 1; i <= m; i++)
993       {  k = head[i]; /* x[k] = xB[i] */
994          if (l[k] != -DBL_MAX && beta[i] < l[k])
995             sum += l[k] - beta[i];
996          if (u[k] != +DBL_MAX && beta[i] > u[k])
997             sum += beta[i] - u[k];
998       }
999       return sum;
1000 }
1001 
1002 /***********************************************************************
1003 *  display - display search progress
1004 *
1005 *  This routine displays some information about the search progress
1006 *  that includes:
1007 *
1008 *  search phase;
1009 *
1010 *  number of simplex iterations performed by the solver;
1011 *
1012 *  original objective value;
1013 *
1014 *  sum of (scaled) primal infeasibilities;
1015 *
1016 *  number of infeasibilities (phase I) or non-optimalities (phase II);
1017 *
1018 *  number of basic factorizations since last display output. */
1019 
display(struct csa * csa,int spec)1020 static void display(struct csa *csa, int spec)
1021 {     int nnn, k;
1022       double obj, sum, *save, *save1;
1023 #if 1 /* 15/VII-2017 */
1024       double tm_cur;
1025 #endif
1026       /* check if the display output should be skipped */
1027       if (csa->msg_lev < GLP_MSG_ON) goto skip;
1028 #if 1 /* 15/VII-2017 */
1029       tm_cur = xtime();
1030 #endif
1031       if (csa->out_dly > 0 &&
1032 #if 0 /* 15/VII-2017 */
1033          1000.0 * xdifftime(xtime(), csa->tm_beg) < csa->out_dly)
1034 #else
1035          1000.0 * xdifftime(tm_cur, csa->tm_beg) < csa->out_dly)
1036 #endif
1037          goto skip;
1038       if (csa->it_cnt == csa->it_dpy) goto skip;
1039 #if 0 /* 15/VII-2017 */
1040       if (!spec && csa->it_cnt % csa->out_frq != 0) goto skip;
1041 #else
1042       if (!spec &&
1043          1000.0 * xdifftime(tm_cur, csa->tm_dpy) < csa->out_frq)
1044          goto skip;
1045 #endif
1046       /* compute original objective value */
1047       save = csa->lp->c;
1048       csa->lp->c = csa->orig_c;
1049       obj = csa->dir * spx_eval_obj(csa->lp, csa->beta);
1050       csa->lp->c = save;
1051 #if SCALE_Z
1052       obj *= csa->fz;
1053 #endif
1054       /* compute sum of (scaled) primal infeasibilities */
1055 #if 1 /* 01/VII-2017 */
1056       save = csa->lp->l;
1057       save1 = csa->lp->u;
1058       csa->lp->l = csa->orig_l;
1059       csa->lp->u = csa->orig_u;
1060 #endif
1061       sum = sum_infeas(csa->lp, csa->beta);
1062 #if 1 /* 01/VII-2017 */
1063       csa->lp->l = save;
1064       csa->lp->u = save1;
1065 #endif
1066       /* compute number of infeasibilities/non-optimalities */
1067       switch (csa->phase)
1068       {  case 1:
1069             nnn = 0;
1070             for (k = 1; k <= csa->lp->n; k++)
1071                if (csa->lp->c[k] != 0.0) nnn++;
1072             break;
1073          case 2:
1074             xassert(csa->d_st);
1075             nnn = spx_chuzc_sel(csa->lp, csa->d, csa->tol_dj,
1076                csa->tol_dj1, NULL);
1077             break;
1078          default:
1079             xassert(csa != csa);
1080       }
1081       /* display search progress */
1082       xprintf("%c%6d: obj = %17.9e inf = %11.3e (%d)",
1083          csa->phase == 2 ? '*' : ' ', csa->it_cnt, obj, sum, nnn);
1084       if (csa->inv_cnt)
1085       {  /* number of basis factorizations performed */
1086          xprintf(" %d", csa->inv_cnt);
1087          csa->inv_cnt = 0;
1088       }
1089 #if 1 /* 23/VI-2017 */
1090       if (csa->phase == 1 && csa->r_test == GLP_RT_FLIP)
1091       {  /*xprintf("   %d,%d", csa->ns_cnt, csa->ls_cnt);*/
1092          if (csa->ns_cnt + csa->ls_cnt)
1093             xprintf(" %d%%",
1094                (100 * csa->ls_cnt) / (csa->ns_cnt + csa->ls_cnt));
1095          csa->ns_cnt = csa->ls_cnt = 0;
1096       }
1097 #endif
1098       xprintf("\n");
1099       csa->it_dpy = csa->it_cnt;
1100 #if 1 /* 15/VII-2017 */
1101       csa->tm_dpy = tm_cur;
1102 #endif
1103 skip: return;
1104 }
1105 
1106 /***********************************************************************
1107 *  spx_primal - driver to the primal simplex method
1108 *
1109 *  This routine is a driver to the two-phase primal simplex method.
1110 *
1111 *  On exit this routine returns one of the following codes:
1112 *
1113 *  0  LP instance has been successfully solved.
1114 *
1115 *  GLP_EITLIM
1116 *     Iteration limit has been exhausted.
1117 *
1118 *  GLP_ETMLIM
1119 *     Time limit has been exhausted.
1120 *
1121 *  GLP_EFAIL
1122 *     The solver failed to solve LP instance. */
1123 
primal_simplex(struct csa * csa)1124 static int primal_simplex(struct csa *csa)
1125 {     /* primal simplex method main logic routine */
1126       SPXLP *lp = csa->lp;
1127       int m = lp->m;
1128       int n = lp->n;
1129       double *c = lp->c;
1130       int *head = lp->head;
1131       SPXAT *at = csa->at;
1132       SPXNT *nt = csa->nt;
1133       double *beta = csa->beta;
1134       double *d = csa->d;
1135       SPXSE *se = csa->se;
1136       int *list = csa->list;
1137 #if 0 /* 11/VI-2017 */
1138       double *tcol = csa->tcol;
1139       double *trow = csa->trow;
1140 #endif
1141 #if 0 /* 09/VII-2017 */
1142       double *pi = csa->work;
1143       double *rho = csa->work;
1144 #else
1145       double *pi = csa->work.vec;
1146       double *rho = csa->work.vec;
1147 #endif
1148       int msg_lev = csa->msg_lev;
1149       double tol_bnd = csa->tol_bnd;
1150       double tol_bnd1 = csa->tol_bnd1;
1151       double tol_dj = csa->tol_dj;
1152       double tol_dj1 = csa->tol_dj1;
1153       int perturb = -1;
1154       /* -1 = perturbation is not used, but enabled
1155        *  0 = perturbation is not used and disabled
1156        * +1 = perturbation is being used */
1157       int j, refct, ret;
1158 loop: /* main loop starts here */
1159       /* compute factorization of the basis matrix */
1160       if (!lp->valid)
1161       {  double cond;
1162          ret = spx_factorize(lp);
1163          csa->inv_cnt++;
1164          if (ret != 0)
1165          {  if (msg_lev >= GLP_MSG_ERR)
1166                xprintf("Error: unable to factorize the basis matrix (%d"
1167                   ")\n", ret);
1168             csa->p_stat = csa->d_stat = GLP_UNDEF;
1169             ret = GLP_EFAIL;
1170             goto fini;
1171          }
1172          /* check condition of the basis matrix */
1173          cond = bfd_condest(lp->bfd);
1174          if (cond > 1.0 / DBL_EPSILON)
1175          {  if (msg_lev >= GLP_MSG_ERR)
1176                xprintf("Error: basis matrix is singular to working prec"
1177                   "ision (cond = %.3g)\n", cond);
1178             csa->p_stat = csa->d_stat = GLP_UNDEF;
1179             ret = GLP_EFAIL;
1180             goto fini;
1181          }
1182          if (cond > 0.001 / DBL_EPSILON)
1183          {  if (msg_lev >= GLP_MSG_ERR)
1184                xprintf("Warning: basis matrix is ill-conditioned (cond "
1185                   "= %.3g)\n", cond);
1186          }
1187          /* invalidate basic solution components */
1188          csa->beta_st = csa->d_st = 0;
1189       }
1190       /* compute values of basic variables beta = (beta[i]) */
1191       if (!csa->beta_st)
1192       {  spx_eval_beta(lp, beta);
1193          csa->beta_st = 1; /* just computed */
1194          /* determine the search phase, if not determined yet */
1195          if (!csa->phase)
1196          {  if (set_penalty(csa, 0.97 * tol_bnd, 0.97 * tol_bnd1))
1197             {  /* current basic solution is primal infeasible */
1198                /* start to minimize the sum of infeasibilities */
1199                csa->phase = 1;
1200             }
1201             else
1202             {  /* current basic solution is primal feasible */
1203                /* start to minimize the original objective function */
1204                csa->phase = 2;
1205                memcpy(c, csa->orig_c, (1+n) * sizeof(double));
1206             }
1207             /* working objective coefficients have been changed, so
1208              * invalidate reduced costs */
1209             csa->d_st = 0;
1210          }
1211          /* make sure that the current basic solution remains primal
1212           * feasible (or pseudo-feasible on phase I) */
1213          if (perturb <= 0)
1214          {  if (check_feas(csa, csa->phase, tol_bnd, tol_bnd1))
1215             {  /* excessive bound violations due to round-off errors */
1216 #if 1 /* 01/VII-2017 */
1217                if (perturb < 0)
1218                {  if (msg_lev >= GLP_MSG_ALL)
1219                      xprintf("Perturbing LP to avoid instability [%d].."
1220                         ".\n", csa->it_cnt);
1221                   perturb = 1;
1222                   goto loop;
1223                }
1224 #endif
1225                if (msg_lev >= GLP_MSG_ERR)
1226                   xprintf("Warning: numerical instability (primal simpl"
1227                      "ex, phase %s)\n", csa->phase == 1 ? "I" : "II");
1228                /* restart the search */
1229                lp->valid = 0;
1230                csa->phase = 0;
1231                goto loop;
1232             }
1233             if (csa->phase == 1)
1234             {  int i, cnt;
1235                for (i = 1; i <= m; i++)
1236                   csa->tcol.ind[i] = i;
1237                cnt = adjust_penalty(csa, m, csa->tcol.ind,
1238                   0.99 * tol_bnd, 0.99 * tol_bnd1);
1239                if (cnt)
1240                {  /*xprintf("*** cnt = %d\n", cnt);*/
1241                   csa->d_st = 0;
1242                }
1243             }
1244          }
1245          else
1246          {  /* FIXME */
1247             play_bounds(csa, 1);
1248          }
1249       }
1250       /* at this point the search phase is determined */
1251       xassert(csa->phase == 1 || csa->phase == 2);
1252       /* compute reduced costs of non-basic variables d = (d[j]) */
1253       if (!csa->d_st)
1254       {  spx_eval_pi(lp, pi);
1255          for (j = 1; j <= n-m; j++)
1256             d[j] = spx_eval_dj(lp, pi, j);
1257          csa->d_st = 1; /* just computed */
1258       }
1259       /* reset the reference space, if necessary */
1260       if (se != NULL && !se->valid)
1261          spx_reset_refsp(lp, se), refct = 1000;
1262       /* at this point the basis factorization and all basic solution
1263        * components are valid */
1264       xassert(lp->valid && csa->beta_st && csa->d_st);
1265 #if CHECK_ACCURACY
1266       /* check accuracy of current basic solution components (only for
1267        * debugging) */
1268       check_accuracy(csa);
1269 #endif
1270       /* check if the iteration limit has been exhausted */
1271       if (csa->it_cnt - csa->it_beg >= csa->it_lim)
1272       {  if (perturb > 0)
1273          {  /* remove perturbation */
1274             remove_perturb(csa);
1275             perturb = 0;
1276          }
1277          if (csa->beta_st != 1)
1278             csa->beta_st = 0;
1279          if (csa->d_st != 1)
1280             csa->d_st = 0;
1281          if (!(csa->beta_st && csa->d_st))
1282             goto loop;
1283          display(csa, 1);
1284          if (msg_lev >= GLP_MSG_ALL)
1285             xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
1286          csa->p_stat = (csa->phase == 2 ? GLP_FEAS : GLP_INFEAS);
1287          csa->d_stat = GLP_UNDEF; /* will be set below */
1288          ret = GLP_EITLIM;
1289          goto fini;
1290       }
1291       /* check if the time limit has been exhausted */
1292       if (1000.0 * xdifftime(xtime(), csa->tm_beg) >= csa->tm_lim)
1293       {  if (perturb > 0)
1294          {  /* remove perturbation */
1295             remove_perturb(csa);
1296             perturb = 0;
1297          }
1298          if (csa->beta_st != 1)
1299             csa->beta_st = 0;
1300          if (csa->d_st != 1)
1301             csa->d_st = 0;
1302          if (!(csa->beta_st && csa->d_st))
1303             goto loop;
1304          display(csa, 1);
1305          if (msg_lev >= GLP_MSG_ALL)
1306             xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
1307          csa->p_stat = (csa->phase == 2 ? GLP_FEAS : GLP_INFEAS);
1308          csa->d_stat = GLP_UNDEF; /* will be set below */
1309          ret = GLP_ETMLIM;
1310          goto fini;
1311       }
1312       /* display the search progress */
1313       display(csa, 0);
1314       /* select eligible non-basic variables */
1315       switch (csa->phase)
1316       {  case 1:
1317             csa->num = spx_chuzc_sel(lp, d, 1e-8, 0.0, list);
1318             break;
1319          case 2:
1320             csa->num = spx_chuzc_sel(lp, d, tol_dj, tol_dj1, list);
1321             break;
1322          default:
1323             xassert(csa != csa);
1324       }
1325       /* check for optimality */
1326       if (csa->num == 0)
1327       {  if (perturb > 0 && csa->phase == 2)
1328          {  /* remove perturbation */
1329             remove_perturb(csa);
1330             perturb = 0;
1331          }
1332          if (csa->beta_st != 1)
1333             csa->beta_st = 0;
1334          if (csa->d_st != 1)
1335             csa->d_st = 0;
1336          if (!(csa->beta_st && csa->d_st))
1337             goto loop;
1338          /* current basis is optimal */
1339          display(csa, 1);
1340          switch (csa->phase)
1341          {  case 1:
1342                /* check for primal feasibility */
1343                if (!check_feas(csa, 2, tol_bnd, tol_bnd1))
1344                {  /* feasible solution found; switch to phase II */
1345                   memcpy(c, csa->orig_c, (1+n) * sizeof(double));
1346                   csa->phase = 2;
1347                   csa->d_st = 0;
1348                   goto loop;
1349                }
1350                /* no feasible solution exists */
1351 #if 1 /* 09/VII-2017 */
1352                /* FIXME: remove perturbation */
1353 #endif
1354                if (msg_lev >= GLP_MSG_ALL)
1355                   xprintf("LP HAS NO PRIMAL FEASIBLE SOLUTION\n");
1356                csa->p_stat = GLP_NOFEAS;
1357                csa->d_stat = GLP_UNDEF; /* will be set below */
1358                ret = 0;
1359                goto fini;
1360             case 2:
1361                /* optimal solution found */
1362                if (msg_lev >= GLP_MSG_ALL)
1363                   xprintf("OPTIMAL LP SOLUTION FOUND\n");
1364                csa->p_stat = csa->d_stat = GLP_FEAS;
1365                ret = 0;
1366                goto fini;
1367             default:
1368                xassert(csa != csa);
1369          }
1370       }
1371       /* choose xN[q] and xB[p] */
1372 #if 0 /* 23/VI-2017 */
1373 #if 0 /* 17/III-2016 */
1374       choose_pivot(csa);
1375 #else
1376       if (choose_pivot(csa) < 0)
1377       {  lp->valid = 0;
1378          goto loop;
1379       }
1380 #endif
1381 #else
1382       ret = choose_pivot(csa);
1383       if (ret < 0)
1384       {  lp->valid = 0;
1385          goto loop;
1386       }
1387       if (ret == 0)
1388          csa->ns_cnt++;
1389       else
1390          csa->ls_cnt++;
1391 #endif
1392       /* check for unboundedness */
1393       if (csa->p == 0)
1394       {  if (perturb > 0)
1395          {  /* remove perturbation */
1396             remove_perturb(csa);
1397             perturb = 0;
1398          }
1399          if (csa->beta_st != 1)
1400             csa->beta_st = 0;
1401          if (csa->d_st != 1)
1402             csa->d_st = 0;
1403          if (!(csa->beta_st && csa->d_st))
1404             goto loop;
1405          display(csa, 1);
1406          switch (csa->phase)
1407          {  case 1:
1408                /* this should never happen */
1409                if (msg_lev >= GLP_MSG_ERR)
1410                   xprintf("Error: primal simplex failed\n");
1411                csa->p_stat = csa->d_stat = GLP_UNDEF;
1412                ret = GLP_EFAIL;
1413                goto fini;
1414             case 2:
1415                /* primal unboundedness detected */
1416                if (msg_lev >= GLP_MSG_ALL)
1417                   xprintf("LP HAS UNBOUNDED PRIMAL SOLUTION\n");
1418                csa->p_stat = GLP_FEAS;
1419                csa->d_stat = GLP_NOFEAS;
1420                ret = 0;
1421                goto fini;
1422             default:
1423                xassert(csa != csa);
1424          }
1425       }
1426 #if 1 /* 01/VII-2017 */
1427       /* check for stalling */
1428       if (csa->p > 0)
1429       {  int k;
1430          xassert(1 <= csa->p && csa->p <= m);
1431          k = head[csa->p]; /* x[k] = xB[p] */
1432          if (lp->l[k] != lp->u[k])
1433          {  if (csa->p_flag)
1434             {  /* xB[p] goes to its upper bound */
1435                xassert(lp->u[k] != +DBL_MAX);
1436                if (fabs(beta[csa->p] - lp->u[k]) >= 1e-6)
1437                {  csa->degen = 0;
1438                   goto skip1;
1439                }
1440             }
1441             else if (lp->l[k] == -DBL_MAX)
1442             {  /* unusual case */
1443                goto skip1;
1444             }
1445             else
1446             {  /* xB[p] goes to its lower bound */
1447                xassert(lp->l[k] != -DBL_MAX);
1448                if (fabs(beta[csa->p] - lp->l[k]) >= 1e-6)
1449                {  csa->degen = 0;
1450                   goto skip1;
1451                }
1452             }
1453             /* degenerate iteration has been detected */
1454             csa->degen++;
1455             if (perturb < 0 && csa->degen >= 200)
1456             {  if (msg_lev >= GLP_MSG_ALL)
1457                   xprintf("Perturbing LP to avoid stalling [%d]...\n",
1458                      csa->it_cnt);
1459                perturb = 1;
1460             }
1461 skip1:      ;
1462          }
1463       }
1464 #endif
1465       /* update values of basic variables for adjacent basis */
1466 #if 0 /* 11/VI-2017 */
1467       spx_update_beta(lp, beta, csa->p, csa->p_flag, csa->q, tcol);
1468 #else
1469       spx_update_beta_s(lp, beta, csa->p, csa->p_flag, csa->q,
1470          &csa->tcol);
1471 #endif
1472       csa->beta_st = 2;
1473       /* p < 0 means that xN[q] jumps to its opposite bound */
1474       if (csa->p < 0)
1475          goto skip;
1476       /* xN[q] enters and xB[p] leaves the basis */
1477       /* compute p-th row of inv(B) */
1478       spx_eval_rho(lp, csa->p, rho);
1479       /* compute p-th (pivot) row of the simplex table */
1480 #if 0 /* 11/VI-2017 */
1481       if (at != NULL)
1482          spx_eval_trow1(lp, at, rho, trow);
1483       else
1484          spx_nt_prod(lp, nt, trow, 1, -1.0, rho);
1485 #else
1486       if (at != NULL)
1487          spx_eval_trow1(lp, at, rho, csa->trow.vec);
1488       else
1489          spx_nt_prod(lp, nt, csa->trow.vec, 1, -1.0, rho);
1490       fvs_gather_vec(&csa->trow, DBL_EPSILON);
1491 #endif
1492       /* FIXME: tcol[p] and trow[q] should be close to each other */
1493 #if 0 /* 26/V-2017 by cmatraki */
1494       xassert(trow[csa->q] != 0.0);
1495 #else
1496       if (csa->trow.vec[csa->q] == 0.0)
1497       {  if (msg_lev >= GLP_MSG_ERR)
1498             xprintf("Error: trow[q] = 0.0\n");
1499          csa->p_stat = csa->d_stat = GLP_UNDEF;
1500          ret = GLP_EFAIL;
1501          goto fini;
1502       }
1503 #endif
1504       /* update reduced costs of non-basic variables for adjacent
1505        * basis */
1506 #if 1 /* 23/VI-2017 */
1507       /* dual solution may be invalidated due to long step */
1508       if (csa->d_st)
1509 #endif
1510 #if 0 /* 11/VI-2017 */
1511       if (spx_update_d(lp, d, csa->p, csa->q, trow, tcol) <= 1e-9)
1512 #else
1513       if (spx_update_d_s(lp, d, csa->p, csa->q, &csa->trow, &csa->tcol)
1514          <= 1e-9)
1515 #endif
1516       {  /* successful updating */
1517          csa->d_st = 2;
1518          if (csa->phase == 1)
1519          {  /* adjust reduced cost of xN[q] in adjacent basis, since
1520              * its penalty coefficient changes (see below) */
1521             d[csa->q] -= c[head[csa->p]];
1522          }
1523       }
1524       else
1525       {  /* new reduced costs are inaccurate */
1526          csa->d_st = 0;
1527       }
1528       if (csa->phase == 1)
1529       {  /* xB[p] leaves the basis replacing xN[q], so set its penalty
1530           * coefficient to zero */
1531          c[head[csa->p]] = 0.0;
1532       }
1533       /* update steepest edge weights for adjacent basis, if used */
1534       if (se != NULL)
1535       {  if (refct > 0)
1536 #if 0 /* 11/VI-2017 */
1537          {  if (spx_update_gamma(lp, se, csa->p, csa->q, trow, tcol)
1538                <= 1e-3)
1539 #else /* FIXME: spx_update_gamma_s */
1540          {  if (spx_update_gamma(lp, se, csa->p, csa->q, csa->trow.vec,
1541                csa->tcol.vec) <= 1e-3)
1542 #endif
1543             {  /* successful updating */
1544                refct--;
1545             }
1546             else
1547             {  /* new weights are inaccurate; reset reference space */
1548                se->valid = 0;
1549             }
1550          }
1551          else
1552          {  /* too many updates; reset reference space */
1553             se->valid = 0;
1554          }
1555       }
1556       /* update matrix N for adjacent basis, if used */
1557       if (nt != NULL)
1558          spx_update_nt(lp, nt, csa->p, csa->q);
1559 skip: /* change current basis header to adjacent one */
1560       spx_change_basis(lp, csa->p, csa->p_flag, csa->q);
1561       /* and update factorization of the basis matrix */
1562       if (csa->p > 0)
1563          spx_update_invb(lp, csa->p, head[csa->p]);
1564 #if 1
1565       if (perturb <= 0)
1566       {  if (csa->phase == 1)
1567          {  int cnt;
1568             /* adjust penalty function coefficients */
1569             cnt = adjust_penalty(csa, csa->tcol.nnz, csa->tcol.ind,
1570                0.99 * tol_bnd, 0.99 * tol_bnd1);
1571             if (cnt)
1572             {  /* some coefficients were changed, so invalidate reduced
1573                 * costs of non-basic variables */
1574                /*xprintf("... cnt = %d\n", cnt);*/
1575                csa->d_st = 0;
1576             }
1577          }
1578       }
1579       else
1580       {  /* FIXME */
1581          play_bounds(csa, 0);
1582       }
1583 #endif
1584       /* simplex iteration complete */
1585       csa->it_cnt++;
1586       goto loop;
1587 fini: /* restore original objective function */
1588       memcpy(c, csa->orig_c, (1+n) * sizeof(double));
1589       /* compute reduced costs of non-basic variables and determine
1590        * solution dual status, if necessary */
1591       if (csa->p_stat != GLP_UNDEF && csa->d_stat == GLP_UNDEF)
1592       {  xassert(ret != GLP_EFAIL);
1593          spx_eval_pi(lp, pi);
1594          for (j = 1; j <= n-m; j++)
1595             d[j] = spx_eval_dj(lp, pi, j);
1596          csa->num = spx_chuzc_sel(lp, d, tol_dj, tol_dj1, NULL);
1597          csa->d_stat = (csa->num == 0 ? GLP_FEAS : GLP_INFEAS);
1598       }
1599       return ret;
1600 }
1601 
1602 int spx_primal(glp_prob *P, const glp_smcp *parm)
1603 {     /* driver to the primal simplex method */
1604       struct csa csa_, *csa = &csa_;
1605       SPXLP lp;
1606       SPXAT at;
1607       SPXNT nt;
1608       SPXSE se;
1609       int ret, *map, *daeh;
1610 #if SCALE_Z
1611       int i, j, k;
1612 #endif
1613       /* build working LP and its initial basis */
1614       memset(csa, 0, sizeof(struct csa));
1615       csa->lp = &lp;
1616       spx_init_lp(csa->lp, P, parm->excl);
1617       spx_alloc_lp(csa->lp);
1618       map = talloc(1+P->m+P->n, int);
1619       spx_build_lp(csa->lp, P, parm->excl, parm->shift, map);
1620       spx_build_basis(csa->lp, P, map);
1621       switch (P->dir)
1622       {  case GLP_MIN:
1623             csa->dir = +1;
1624             break;
1625          case GLP_MAX:
1626             csa->dir = -1;
1627             break;
1628          default:
1629             xassert(P != P);
1630       }
1631 #if SCALE_Z
1632       csa->fz = 0.0;
1633       for (k = 1; k <= csa->lp->n; k++)
1634       {  double t = fabs(csa->lp->c[k]);
1635          if (csa->fz < t)
1636             csa->fz = t;
1637       }
1638       if (csa->fz <= 1000.0)
1639          csa->fz = 1.0;
1640       else
1641          csa->fz /= 1000.0;
1642       /*xprintf("csa->fz = %g\n", csa->fz);*/
1643       for (k = 0; k <= csa->lp->n; k++)
1644          csa->lp->c[k] /= csa->fz;
1645 #endif
1646       csa->orig_c = talloc(1+csa->lp->n, double);
1647       memcpy(csa->orig_c, csa->lp->c, (1+csa->lp->n) * sizeof(double));
1648 #if 1 /*PERTURB*/
1649       csa->orig_l = talloc(1+csa->lp->n, double);
1650       memcpy(csa->orig_l, csa->lp->l, (1+csa->lp->n) * sizeof(double));
1651       csa->orig_u = talloc(1+csa->lp->n, double);
1652       memcpy(csa->orig_u, csa->lp->u, (1+csa->lp->n) * sizeof(double));
1653 #else
1654       csa->orig_l = csa->orig_u = NULL;
1655 #endif
1656       switch (parm->aorn)
1657       {  case GLP_USE_AT:
1658             /* build matrix A in row-wise format */
1659             csa->at = &at;
1660             csa->nt = NULL;
1661             spx_alloc_at(csa->lp, csa->at);
1662             spx_build_at(csa->lp, csa->at);
1663             break;
1664          case GLP_USE_NT:
1665             /* build matrix N in row-wise format for initial basis */
1666             csa->at = NULL;
1667             csa->nt = &nt;
1668             spx_alloc_nt(csa->lp, csa->nt);
1669             spx_init_nt(csa->lp, csa->nt);
1670             spx_build_nt(csa->lp, csa->nt);
1671             break;
1672          default:
1673             xassert(parm != parm);
1674       }
1675       /* allocate and initialize working components */
1676       csa->phase = 0;
1677       csa->beta = talloc(1+csa->lp->m, double);
1678       csa->beta_st = 0;
1679       csa->d = talloc(1+csa->lp->n-csa->lp->m, double);
1680       csa->d_st = 0;
1681       switch (parm->pricing)
1682       {  case GLP_PT_STD:
1683             csa->se = NULL;
1684             break;
1685          case GLP_PT_PSE:
1686             csa->se = &se;
1687             spx_alloc_se(csa->lp, csa->se);
1688             break;
1689          default:
1690             xassert(parm != parm);
1691       }
1692       csa->list = talloc(1+csa->lp->n-csa->lp->m, int);
1693 #if 0 /* 11/VI-2017 */
1694       csa->tcol = talloc(1+csa->lp->m, double);
1695       csa->trow = talloc(1+csa->lp->n-csa->lp->m, double);
1696 #else
1697       fvs_alloc_vec(&csa->tcol, csa->lp->m);
1698       fvs_alloc_vec(&csa->trow, csa->lp->n-csa->lp->m);
1699 #endif
1700 #if 1 /* 23/VI-2017 */
1701       csa->bp = NULL;
1702 #endif
1703 #if 0 /* 09/VII-2017 */
1704       csa->work = talloc(1+csa->lp->m, double);
1705 #else
1706       fvs_alloc_vec(&csa->work, csa->lp->m);
1707 #endif
1708       /* initialize control parameters */
1709       csa->msg_lev = parm->msg_lev;
1710 #if 0 /* 23/VI-2017 */
1711       switch (parm->r_test)
1712       {  case GLP_RT_STD:
1713             csa->harris = 0;
1714             break;
1715          case GLP_RT_HAR:
1716 #if 1 /* 16/III-2016 */
1717          case GLP_RT_FLIP:
1718             /* FIXME */
1719             /* currently for primal simplex GLP_RT_FLIP is equivalent
1720              * to GLP_RT_HAR */
1721 #endif
1722             csa->harris = 1;
1723             break;
1724          default:
1725             xassert(parm != parm);
1726       }
1727 #else
1728       switch (parm->r_test)
1729       {  case GLP_RT_STD:
1730          case GLP_RT_HAR:
1731             break;
1732          case GLP_RT_FLIP:
1733             csa->bp = talloc(1+2*csa->lp->m+1, SPXBP);
1734             break;
1735          default:
1736             xassert(parm != parm);
1737       }
1738       csa->r_test = parm->r_test;
1739 #endif
1740       csa->tol_bnd = parm->tol_bnd;
1741       csa->tol_bnd1 = .001 * parm->tol_bnd;
1742       csa->tol_dj = parm->tol_dj;
1743       csa->tol_dj1 = .001 * parm->tol_dj;
1744       csa->tol_piv = parm->tol_piv;
1745       csa->it_lim = parm->it_lim;
1746       csa->tm_lim = parm->tm_lim;
1747       csa->out_frq = parm->out_frq;
1748       csa->out_dly = parm->out_dly;
1749       /* initialize working parameters */
1750       csa->tm_beg = xtime();
1751       csa->it_beg = csa->it_cnt = P->it_cnt;
1752       csa->it_dpy = -1;
1753 #if 1 /* 15/VII-2017 */
1754       csa->tm_dpy = 0.0;
1755 #endif
1756       csa->inv_cnt = 0;
1757 #if 1 /* 01/VII-2017 */
1758       csa->degen = 0;
1759 #endif
1760 #if 1 /* 23/VI-2017 */
1761       csa->ns_cnt = csa->ls_cnt = 0;
1762 #endif
1763       /* try to solve working LP */
1764       ret = primal_simplex(csa);
1765       /* return basis factorization back to problem object */
1766       P->valid = csa->lp->valid;
1767       P->bfd = csa->lp->bfd;
1768       /* set solution status */
1769       P->pbs_stat = csa->p_stat;
1770       P->dbs_stat = csa->d_stat;
1771       /* if the solver failed, do not store basis header and basic
1772        * solution components to problem object */
1773       if (ret == GLP_EFAIL)
1774          goto skip;
1775       /* convert working LP basis to original LP basis and store it to
1776        * problem object */
1777       daeh = talloc(1+csa->lp->n, int);
1778       spx_store_basis(csa->lp, P, map, daeh);
1779       /* compute simplex multipliers for final basic solution found by
1780        * the solver */
1781 #if 0 /* 09/VII-2017 */
1782       spx_eval_pi(csa->lp, csa->work);
1783 #else
1784       spx_eval_pi(csa->lp, csa->work.vec);
1785 #endif
1786       /* convert working LP solution to original LP solution and store
1787        * it into the problem object */
1788 #if SCALE_Z
1789       for (i = 1; i <= csa->lp->m; i++)
1790          csa->work.vec[i] *= csa->fz;
1791       for (j = 1; j <= csa->lp->n-csa->lp->m; j++)
1792          csa->d[j] *= csa->fz;
1793 #endif
1794 #if 0 /* 09/VII-2017 */
1795       spx_store_sol(csa->lp, P, SHIFT, map, daeh, csa->beta, csa->work,
1796          csa->d);
1797 #else
1798       spx_store_sol(csa->lp, P, parm->shift, map, daeh, csa->beta,
1799          csa->work.vec, csa->d);
1800 #endif
1801       tfree(daeh);
1802       /* save simplex iteration count */
1803       P->it_cnt = csa->it_cnt;
1804       /* report auxiliary/structural variable causing unboundedness */
1805       P->some = 0;
1806       if (csa->p_stat == GLP_FEAS && csa->d_stat == GLP_NOFEAS)
1807       {  int k, kk;
1808          /* xN[q] = x[k] causes unboundedness */
1809          xassert(1 <= csa->q && csa->q <= csa->lp->n - csa->lp->m);
1810          k = csa->lp->head[csa->lp->m + csa->q];
1811          xassert(1 <= k && k <= csa->lp->n);
1812          /* convert to number of original variable */
1813          for (kk = 1; kk <= P->m + P->n; kk++)
1814          {  if (abs(map[kk]) == k)
1815             {  P->some = kk;
1816                break;
1817             }
1818          }
1819          xassert(P->some != 0);
1820       }
1821 skip: /* deallocate working objects and arrays */
1822       spx_free_lp(csa->lp);
1823       tfree(map);
1824       tfree(csa->orig_c);
1825 #if 1 /*PERTURB*/
1826       tfree(csa->orig_l);
1827       tfree(csa->orig_u);
1828 #endif
1829       if (csa->at != NULL)
1830          spx_free_at(csa->lp, csa->at);
1831       if (csa->nt != NULL)
1832          spx_free_nt(csa->lp, csa->nt);
1833       tfree(csa->beta);
1834       tfree(csa->d);
1835       if (csa->se != NULL)
1836          spx_free_se(csa->lp, csa->se);
1837       tfree(csa->list);
1838 #if 0 /* 11/VI-2017 */
1839       tfree(csa->tcol);
1840       tfree(csa->trow);
1841 #else
1842       fvs_free_vec(&csa->tcol);
1843       fvs_free_vec(&csa->trow);
1844 #endif
1845 #if 1 /* 23/VI-2017 */
1846       if (csa->bp != NULL)
1847          tfree(csa->bp);
1848 #endif
1849 #if 0 /* 09/VII-2017 */
1850       tfree(csa->work);
1851 #else
1852       fvs_free_vec(&csa->work);
1853 #endif
1854       /* return to calling program */
1855       return ret;
1856 }
1857 
1858 /* eof */
1859