lemon-project-template-glpk

comparison deps/glpk/src/glpspx01.c @ 11:4fc6ad2fb8a6

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