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