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