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