1 /* glpapi12.c (basis factorization and simplex tableau routines) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 * Copyright (C) 2000-2013 Free Software Foundation, Inc.
6 * Written by Andrew Makhorin <mao@gnu.org>.
7 *
8 * GLPK is free software: you can redistribute it and/or modify it
9 * under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * GLPK is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
20 ***********************************************************************/
21
22 #include "draft.h"
23 #include "env.h"
24 #include "prob.h"
25
26 /***********************************************************************
27 * NAME
28 *
29 * glp_bf_exists - check if the basis factorization exists
30 *
31 * SYNOPSIS
32 *
33 * int glp_bf_exists(glp_prob *lp);
34 *
35 * RETURNS
36 *
37 * If the basis factorization for the current basis associated with
38 * the specified problem object exists and therefore is available for
39 * computations, the routine glp_bf_exists returns non-zero. Otherwise
40 * the routine returns zero. */
41
glp_bf_exists(glp_prob * lp)42 int glp_bf_exists(glp_prob *lp)
43 { int ret;
44 ret = (lp->m == 0 || lp->valid);
45 return ret;
46 }
47
48 /***********************************************************************
49 * NAME
50 *
51 * glp_factorize - compute the basis factorization
52 *
53 * SYNOPSIS
54 *
55 * int glp_factorize(glp_prob *lp);
56 *
57 * DESCRIPTION
58 *
59 * The routine glp_factorize computes the basis factorization for the
60 * current basis associated with the specified problem object.
61 *
62 * RETURNS
63 *
64 * 0 The basis factorization has been successfully computed.
65 *
66 * GLP_EBADB
67 * The basis matrix is invalid, i.e. the number of basic (auxiliary
68 * and structural) variables differs from the number of rows in the
69 * problem object.
70 *
71 * GLP_ESING
72 * The basis matrix is singular within the working precision.
73 *
74 * GLP_ECOND
75 * The basis matrix is ill-conditioned. */
76
b_col(void * info,int j,int ind[],double val[])77 static int b_col(void *info, int j, int ind[], double val[])
78 { glp_prob *lp = info;
79 int m = lp->m;
80 GLPAIJ *aij;
81 int k, len;
82 xassert(1 <= j && j <= m);
83 /* determine the ordinal number of basic auxiliary or structural
84 variable x[k] corresponding to basic variable xB[j] */
85 k = lp->head[j];
86 /* build j-th column of the basic matrix, which is k-th column of
87 the scaled augmented matrix (I | -R*A*S) */
88 if (k <= m)
89 { /* x[k] is auxiliary variable */
90 len = 1;
91 ind[1] = k;
92 val[1] = 1.0;
93 }
94 else
95 { /* x[k] is structural variable */
96 len = 0;
97 for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next)
98 { len++;
99 ind[len] = aij->row->i;
100 val[len] = - aij->row->rii * aij->val * aij->col->sjj;
101 }
102 }
103 return len;
104 }
105
glp_factorize(glp_prob * lp)106 int glp_factorize(glp_prob *lp)
107 { int m = lp->m;
108 int n = lp->n;
109 GLPROW **row = lp->row;
110 GLPCOL **col = lp->col;
111 int *head = lp->head;
112 int j, k, stat, ret;
113 /* invalidate the basis factorization */
114 lp->valid = 0;
115 /* build the basis header */
116 j = 0;
117 for (k = 1; k <= m+n; k++)
118 { if (k <= m)
119 { stat = row[k]->stat;
120 row[k]->bind = 0;
121 }
122 else
123 { stat = col[k-m]->stat;
124 col[k-m]->bind = 0;
125 }
126 if (stat == GLP_BS)
127 { j++;
128 if (j > m)
129 { /* too many basic variables */
130 ret = GLP_EBADB;
131 goto fini;
132 }
133 head[j] = k;
134 if (k <= m)
135 row[k]->bind = j;
136 else
137 col[k-m]->bind = j;
138 }
139 }
140 if (j < m)
141 { /* too few basic variables */
142 ret = GLP_EBADB;
143 goto fini;
144 }
145 /* try to factorize the basis matrix */
146 if (m > 0)
147 { if (lp->bfd == NULL)
148 { lp->bfd = bfd_create_it();
149 #if 0 /* 08/III-2014 */
150 copy_bfcp(lp);
151 #endif
152 }
153 switch (bfd_factorize(lp->bfd, m, /*lp->head,*/ b_col, lp))
154 { case 0:
155 /* ok */
156 break;
157 case BFD_ESING:
158 /* singular matrix */
159 ret = GLP_ESING;
160 goto fini;
161 case BFD_ECOND:
162 /* ill-conditioned matrix */
163 ret = GLP_ECOND;
164 goto fini;
165 default:
166 xassert(lp != lp);
167 }
168 lp->valid = 1;
169 }
170 /* factorization successful */
171 ret = 0;
172 fini: /* bring the return code to the calling program */
173 return ret;
174 }
175
176 /***********************************************************************
177 * NAME
178 *
179 * glp_bf_updated - check if the basis factorization has been updated
180 *
181 * SYNOPSIS
182 *
183 * int glp_bf_updated(glp_prob *lp);
184 *
185 * RETURNS
186 *
187 * If the basis factorization has been just computed from scratch, the
188 * routine glp_bf_updated returns zero. Otherwise, if the factorization
189 * has been updated one or more times, the routine returns non-zero. */
190
glp_bf_updated(glp_prob * lp)191 int glp_bf_updated(glp_prob *lp)
192 { int cnt;
193 if (!(lp->m == 0 || lp->valid))
194 xerror("glp_bf_update: basis factorization does not exist\n");
195 #if 0 /* 15/XI-2009 */
196 cnt = (lp->m == 0 ? 0 : lp->bfd->upd_cnt);
197 #else
198 cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd));
199 #endif
200 return cnt;
201 }
202
203 /***********************************************************************
204 * NAME
205 *
206 * glp_get_bfcp - retrieve basis factorization control parameters
207 *
208 * SYNOPSIS
209 *
210 * void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
211 *
212 * DESCRIPTION
213 *
214 * The routine glp_get_bfcp retrieves control parameters, which are
215 * used on computing and updating the basis factorization associated
216 * with the specified problem object.
217 *
218 * Current values of control parameters are stored by the routine in
219 * a glp_bfcp structure, which the parameter parm points to. */
220
221 #if 1 /* 08/III-2014 */
glp_get_bfcp(glp_prob * P,glp_bfcp * parm)222 void glp_get_bfcp(glp_prob *P, glp_bfcp *parm)
223 { if (P->bfd == NULL)
224 P->bfd = bfd_create_it();
225 bfd_get_bfcp(P->bfd, parm);
226 return;
227 }
228 #endif
229
230 /***********************************************************************
231 * NAME
232 *
233 * glp_set_bfcp - change basis factorization control parameters
234 *
235 * SYNOPSIS
236 *
237 * void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
238 *
239 * DESCRIPTION
240 *
241 * The routine glp_set_bfcp changes control parameters, which are used
242 * by internal GLPK routines in computing and updating the basis
243 * factorization associated with the specified problem object.
244 *
245 * New values of the control parameters should be passed in a structure
246 * glp_bfcp, which the parameter parm points to.
247 *
248 * The parameter parm can be specified as NULL, in which case all
249 * control parameters are reset to their default values. */
250
251 #if 1 /* 08/III-2014 */
glp_set_bfcp(glp_prob * P,const glp_bfcp * parm)252 void glp_set_bfcp(glp_prob *P, const glp_bfcp *parm)
253 { if (P->bfd == NULL)
254 P->bfd = bfd_create_it();
255 if (parm != NULL)
256 { if (!(parm->type == GLP_BF_LUF + GLP_BF_FT ||
257 parm->type == GLP_BF_LUF + GLP_BF_BG ||
258 parm->type == GLP_BF_LUF + GLP_BF_GR ||
259 parm->type == GLP_BF_BTF + GLP_BF_BG ||
260 parm->type == GLP_BF_BTF + GLP_BF_GR))
261 xerror("glp_set_bfcp: type = 0x%02X; invalid parameter\n",
262 parm->type);
263 if (!(0.0 < parm->piv_tol && parm->piv_tol < 1.0))
264 xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n",
265 parm->piv_tol);
266 if (parm->piv_lim < 1)
267 xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n",
268 parm->piv_lim);
269 if (!(parm->suhl == GLP_ON || parm->suhl == GLP_OFF))
270 xerror("glp_set_bfcp: suhl = %d; invalid parameter\n",
271 parm->suhl);
272 if (!(0.0 <= parm->eps_tol && parm->eps_tol <= 1e-6))
273 xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n",
274 parm->eps_tol);
275 if (!(1 <= parm->nfs_max && parm->nfs_max <= 32767))
276 xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n",
277 parm->nfs_max);
278 if (!(1 <= parm->nrs_max && parm->nrs_max <= 32767))
279 xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n",
280 parm->nrs_max);
281 }
282 bfd_set_bfcp(P->bfd, parm);
283 return;
284 }
285 #endif
286
287 /***********************************************************************
288 * NAME
289 *
290 * glp_get_bhead - retrieve the basis header information
291 *
292 * SYNOPSIS
293 *
294 * int glp_get_bhead(glp_prob *lp, int k);
295 *
296 * DESCRIPTION
297 *
298 * The routine glp_get_bhead returns the basis header information for
299 * the current basis associated with the specified problem object.
300 *
301 * RETURNS
302 *
303 * If xB[k], 1 <= k <= m, is i-th auxiliary variable (1 <= i <= m), the
304 * routine returns i. Otherwise, if xB[k] is j-th structural variable
305 * (1 <= j <= n), the routine returns m+j. Here m is the number of rows
306 * and n is the number of columns in the problem object. */
307
glp_get_bhead(glp_prob * lp,int k)308 int glp_get_bhead(glp_prob *lp, int k)
309 { if (!(lp->m == 0 || lp->valid))
310 xerror("glp_get_bhead: basis factorization does not exist\n");
311 if (!(1 <= k && k <= lp->m))
312 xerror("glp_get_bhead: k = %d; index out of range\n", k);
313 return lp->head[k];
314 }
315
316 /***********************************************************************
317 * NAME
318 *
319 * glp_get_row_bind - retrieve row index in the basis header
320 *
321 * SYNOPSIS
322 *
323 * int glp_get_row_bind(glp_prob *lp, int i);
324 *
325 * RETURNS
326 *
327 * The routine glp_get_row_bind returns the index k of basic variable
328 * xB[k], 1 <= k <= m, which is i-th auxiliary variable, 1 <= i <= m,
329 * in the current basis associated with the specified problem object,
330 * where m is the number of rows. However, if i-th auxiliary variable
331 * is non-basic, the routine returns zero. */
332
glp_get_row_bind(glp_prob * lp,int i)333 int glp_get_row_bind(glp_prob *lp, int i)
334 { if (!(lp->m == 0 || lp->valid))
335 xerror("glp_get_row_bind: basis factorization does not exist\n"
336 );
337 if (!(1 <= i && i <= lp->m))
338 xerror("glp_get_row_bind: i = %d; row number out of range\n",
339 i);
340 return lp->row[i]->bind;
341 }
342
343 /***********************************************************************
344 * NAME
345 *
346 * glp_get_col_bind - retrieve column index in the basis header
347 *
348 * SYNOPSIS
349 *
350 * int glp_get_col_bind(glp_prob *lp, int j);
351 *
352 * RETURNS
353 *
354 * The routine glp_get_col_bind returns the index k of basic variable
355 * xB[k], 1 <= k <= m, which is j-th structural variable, 1 <= j <= n,
356 * in the current basis associated with the specified problem object,
357 * where m is the number of rows, n is the number of columns. However,
358 * if j-th structural variable is non-basic, the routine returns zero.*/
359
glp_get_col_bind(glp_prob * lp,int j)360 int glp_get_col_bind(glp_prob *lp, int j)
361 { if (!(lp->m == 0 || lp->valid))
362 xerror("glp_get_col_bind: basis factorization does not exist\n"
363 );
364 if (!(1 <= j && j <= lp->n))
365 xerror("glp_get_col_bind: j = %d; column number out of range\n"
366 , j);
367 return lp->col[j]->bind;
368 }
369
370 /***********************************************************************
371 * NAME
372 *
373 * glp_ftran - perform forward transformation (solve system B*x = b)
374 *
375 * SYNOPSIS
376 *
377 * void glp_ftran(glp_prob *lp, double x[]);
378 *
379 * DESCRIPTION
380 *
381 * The routine glp_ftran performs forward transformation, i.e. solves
382 * the system B*x = b, where B is the basis matrix corresponding to the
383 * current basis for the specified problem object, x is the vector of
384 * unknowns to be computed, b is the vector of right-hand sides.
385 *
386 * On entry elements of the vector b should be stored in dense format
387 * in locations x[1], ..., x[m], where m is the number of rows. On exit
388 * the routine stores elements of the vector x in the same locations.
389 *
390 * SCALING/UNSCALING
391 *
392 * Let A~ = (I | -A) is the augmented constraint matrix of the original
393 * (unscaled) problem. In the scaled LP problem instead the matrix A the
394 * scaled matrix A" = R*A*S is actually used, so
395 *
396 * A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) =
397 * (1)
398 * = R*(I | A)*S~ = R*A~*S~,
399 *
400 * is the scaled augmented constraint matrix, where R and S are diagonal
401 * scaling matrices used to scale rows and columns of the matrix A, and
402 *
403 * S~ = diag(inv(R) | S) (2)
404 *
405 * is an augmented diagonal scaling matrix.
406 *
407 * By definition:
408 *
409 * A~ = (B | N), (3)
410 *
411 * where B is the basic matrix, which consists of basic columns of the
412 * augmented constraint matrix A~, and N is a matrix, which consists of
413 * non-basic columns of A~. From (1) it follows that:
414 *
415 * A~" = (B" | N") = (R*B*SB | R*N*SN), (4)
416 *
417 * where SB and SN are parts of the augmented scaling matrix S~, which
418 * correspond to basic and non-basic variables, respectively. Therefore
419 *
420 * B" = R*B*SB, (5)
421 *
422 * which is the scaled basis matrix. */
423
glp_ftran(glp_prob * lp,double x[])424 void glp_ftran(glp_prob *lp, double x[])
425 { int m = lp->m;
426 GLPROW **row = lp->row;
427 GLPCOL **col = lp->col;
428 int i, k;
429 /* B*x = b ===> (R*B*SB)*(inv(SB)*x) = R*b ===>
430 B"*x" = b", where b" = R*b, x = SB*x" */
431 if (!(m == 0 || lp->valid))
432 xerror("glp_ftran: basis factorization does not exist\n");
433 /* b" := R*b */
434 for (i = 1; i <= m; i++)
435 x[i] *= row[i]->rii;
436 /* x" := inv(B")*b" */
437 if (m > 0) bfd_ftran(lp->bfd, x);
438 /* x := SB*x" */
439 for (i = 1; i <= m; i++)
440 { k = lp->head[i];
441 if (k <= m)
442 x[i] /= row[k]->rii;
443 else
444 x[i] *= col[k-m]->sjj;
445 }
446 return;
447 }
448
449 /***********************************************************************
450 * NAME
451 *
452 * glp_btran - perform backward transformation (solve system B'*x = b)
453 *
454 * SYNOPSIS
455 *
456 * void glp_btran(glp_prob *lp, double x[]);
457 *
458 * DESCRIPTION
459 *
460 * The routine glp_btran performs backward transformation, i.e. solves
461 * the system B'*x = b, where B' is a matrix transposed to the basis
462 * matrix corresponding to the current basis for the specified problem
463 * problem object, x is the vector of unknowns to be computed, b is the
464 * vector of right-hand sides.
465 *
466 * On entry elements of the vector b should be stored in dense format
467 * in locations x[1], ..., x[m], where m is the number of rows. On exit
468 * the routine stores elements of the vector x in the same locations.
469 *
470 * SCALING/UNSCALING
471 *
472 * See comments to the routine glp_ftran. */
473
glp_btran(glp_prob * lp,double x[])474 void glp_btran(glp_prob *lp, double x[])
475 { int m = lp->m;
476 GLPROW **row = lp->row;
477 GLPCOL **col = lp->col;
478 int i, k;
479 /* B'*x = b ===> (SB*B'*R)*(inv(R)*x) = SB*b ===>
480 (B")'*x" = b", where b" = SB*b, x = R*x" */
481 if (!(m == 0 || lp->valid))
482 xerror("glp_btran: basis factorization does not exist\n");
483 /* b" := SB*b */
484 for (i = 1; i <= m; i++)
485 { k = lp->head[i];
486 if (k <= m)
487 x[i] /= row[k]->rii;
488 else
489 x[i] *= col[k-m]->sjj;
490 }
491 /* x" := inv[(B")']*b" */
492 if (m > 0) bfd_btran(lp->bfd, x);
493 /* x := R*x" */
494 for (i = 1; i <= m; i++)
495 x[i] *= row[i]->rii;
496 return;
497 }
498
499 /***********************************************************************
500 * NAME
501 *
502 * glp_warm_up - "warm up" LP basis
503 *
504 * SYNOPSIS
505 *
506 * int glp_warm_up(glp_prob *P);
507 *
508 * DESCRIPTION
509 *
510 * The routine glp_warm_up "warms up" the LP basis for the specified
511 * problem object using current statuses assigned to rows and columns
512 * (that is, to auxiliary and structural variables).
513 *
514 * This operation includes computing factorization of the basis matrix
515 * (if it does not exist), computing primal and dual components of basic
516 * solution, and determining the solution status.
517 *
518 * RETURNS
519 *
520 * 0 The operation has been successfully performed.
521 *
522 * GLP_EBADB
523 * The basis matrix is invalid, i.e. the number of basic (auxiliary
524 * and structural) variables differs from the number of rows in the
525 * problem object.
526 *
527 * GLP_ESING
528 * The basis matrix is singular within the working precision.
529 *
530 * GLP_ECOND
531 * The basis matrix is ill-conditioned. */
532
glp_warm_up(glp_prob * P)533 int glp_warm_up(glp_prob *P)
534 { GLPROW *row;
535 GLPCOL *col;
536 GLPAIJ *aij;
537 int i, j, type, stat, ret;
538 double eps, temp, *work;
539 /* invalidate basic solution */
540 P->pbs_stat = P->dbs_stat = GLP_UNDEF;
541 P->obj_val = 0.0;
542 P->some = 0;
543 for (i = 1; i <= P->m; i++)
544 { row = P->row[i];
545 row->prim = row->dual = 0.0;
546 }
547 for (j = 1; j <= P->n; j++)
548 { col = P->col[j];
549 col->prim = col->dual = 0.0;
550 }
551 /* compute the basis factorization, if necessary */
552 if (!glp_bf_exists(P))
553 { ret = glp_factorize(P);
554 if (ret != 0) goto done;
555 }
556 /* allocate working array */
557 work = xcalloc(1+P->m, sizeof(double));
558 /* determine and store values of non-basic variables, compute
559 vector (- N * xN) */
560 for (i = 1; i <= P->m; i++)
561 work[i] = 0.0;
562 for (i = 1; i <= P->m; i++)
563 { row = P->row[i];
564 if (row->stat == GLP_BS)
565 continue;
566 else if (row->stat == GLP_NL)
567 row->prim = row->lb;
568 else if (row->stat == GLP_NU)
569 row->prim = row->ub;
570 else if (row->stat == GLP_NF)
571 row->prim = 0.0;
572 else if (row->stat == GLP_NS)
573 row->prim = row->lb;
574 else
575 xassert(row != row);
576 /* N[j] is i-th column of matrix (I|-A) */
577 work[i] -= row->prim;
578 }
579 for (j = 1; j <= P->n; j++)
580 { col = P->col[j];
581 if (col->stat == GLP_BS)
582 continue;
583 else if (col->stat == GLP_NL)
584 col->prim = col->lb;
585 else if (col->stat == GLP_NU)
586 col->prim = col->ub;
587 else if (col->stat == GLP_NF)
588 col->prim = 0.0;
589 else if (col->stat == GLP_NS)
590 col->prim = col->lb;
591 else
592 xassert(col != col);
593 /* N[j] is (m+j)-th column of matrix (I|-A) */
594 if (col->prim != 0.0)
595 { for (aij = col->ptr; aij != NULL; aij = aij->c_next)
596 work[aij->row->i] += aij->val * col->prim;
597 }
598 }
599 /* compute vector of basic variables xB = - inv(B) * N * xN */
600 glp_ftran(P, work);
601 /* store values of basic variables, check primal feasibility */
602 P->pbs_stat = GLP_FEAS;
603 for (i = 1; i <= P->m; i++)
604 { row = P->row[i];
605 if (row->stat != GLP_BS)
606 continue;
607 row->prim = work[row->bind];
608 type = row->type;
609 if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
610 { eps = 1e-6 + 1e-9 * fabs(row->lb);
611 if (row->prim < row->lb - eps)
612 P->pbs_stat = GLP_INFEAS;
613 }
614 if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
615 { eps = 1e-6 + 1e-9 * fabs(row->ub);
616 if (row->prim > row->ub + eps)
617 P->pbs_stat = GLP_INFEAS;
618 }
619 }
620 for (j = 1; j <= P->n; j++)
621 { col = P->col[j];
622 if (col->stat != GLP_BS)
623 continue;
624 col->prim = work[col->bind];
625 type = col->type;
626 if (type == GLP_LO || type == GLP_DB || type == GLP_FX)
627 { eps = 1e-6 + 1e-9 * fabs(col->lb);
628 if (col->prim < col->lb - eps)
629 P->pbs_stat = GLP_INFEAS;
630 }
631 if (type == GLP_UP || type == GLP_DB || type == GLP_FX)
632 { eps = 1e-6 + 1e-9 * fabs(col->ub);
633 if (col->prim > col->ub + eps)
634 P->pbs_stat = GLP_INFEAS;
635 }
636 }
637 /* compute value of the objective function */
638 P->obj_val = P->c0;
639 for (j = 1; j <= P->n; j++)
640 { col = P->col[j];
641 P->obj_val += col->coef * col->prim;
642 }
643 /* build vector cB of objective coefficients at basic variables */
644 for (i = 1; i <= P->m; i++)
645 work[i] = 0.0;
646 for (j = 1; j <= P->n; j++)
647 { col = P->col[j];
648 if (col->stat == GLP_BS)
649 work[col->bind] = col->coef;
650 }
651 /* compute vector of simplex multipliers pi = inv(B') * cB */
652 glp_btran(P, work);
653 /* compute and store reduced costs of non-basic variables d[j] =
654 c[j] - N'[j] * pi, check dual feasibility */
655 P->dbs_stat = GLP_FEAS;
656 for (i = 1; i <= P->m; i++)
657 { row = P->row[i];
658 if (row->stat == GLP_BS)
659 { row->dual = 0.0;
660 continue;
661 }
662 /* N[j] is i-th column of matrix (I|-A) */
663 row->dual = - work[i];
664 #if 0 /* 07/III-2013 */
665 type = row->type;
666 temp = (P->dir == GLP_MIN ? + row->dual : - row->dual);
667 if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
668 (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
669 P->dbs_stat = GLP_INFEAS;
670 #else
671 stat = row->stat;
672 temp = (P->dir == GLP_MIN ? + row->dual : - row->dual);
673 if ((stat == GLP_NF || stat == GLP_NL) && temp < -1e-5 ||
674 (stat == GLP_NF || stat == GLP_NU) && temp > +1e-5)
675 P->dbs_stat = GLP_INFEAS;
676 #endif
677 }
678 for (j = 1; j <= P->n; j++)
679 { col = P->col[j];
680 if (col->stat == GLP_BS)
681 { col->dual = 0.0;
682 continue;
683 }
684 /* N[j] is (m+j)-th column of matrix (I|-A) */
685 col->dual = col->coef;
686 for (aij = col->ptr; aij != NULL; aij = aij->c_next)
687 col->dual += aij->val * work[aij->row->i];
688 #if 0 /* 07/III-2013 */
689 type = col->type;
690 temp = (P->dir == GLP_MIN ? + col->dual : - col->dual);
691 if ((type == GLP_FR || type == GLP_LO) && temp < -1e-5 ||
692 (type == GLP_FR || type == GLP_UP) && temp > +1e-5)
693 P->dbs_stat = GLP_INFEAS;
694 #else
695 stat = col->stat;
696 temp = (P->dir == GLP_MIN ? + col->dual : - col->dual);
697 if ((stat == GLP_NF || stat == GLP_NL) && temp < -1e-5 ||
698 (stat == GLP_NF || stat == GLP_NU) && temp > +1e-5)
699 P->dbs_stat = GLP_INFEAS;
700 #endif
701 }
702 /* free working array */
703 xfree(work);
704 ret = 0;
705 done: return ret;
706 }
707
708 /***********************************************************************
709 * NAME
710 *
711 * glp_eval_tab_row - compute row of the simplex tableau
712 *
713 * SYNOPSIS
714 *
715 * int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]);
716 *
717 * DESCRIPTION
718 *
719 * The routine glp_eval_tab_row computes a row of the current simplex
720 * tableau for the basic variable, which is specified by the number k:
721 * if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
722 * x[k] is (k-m)-th structural variable, where m is number of rows, and
723 * n is number of columns. The current basis must be available.
724 *
725 * The routine stores column indices and numerical values of non-zero
726 * elements of the computed row using sparse format to the locations
727 * ind[1], ..., ind[len] and val[1], ..., val[len], respectively, where
728 * 0 <= len <= n is number of non-zeros returned on exit.
729 *
730 * Element indices stored in the array ind have the same sense as the
731 * index k, i.e. indices 1 to m denote auxiliary variables and indices
732 * m+1 to m+n denote structural ones (all these variables are obviously
733 * non-basic by definition).
734 *
735 * The computed row shows how the specified basic variable x[k] = xB[i]
736 * depends on non-basic variables:
737 *
738 * xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n],
739 *
740 * where alfa[i,j] are elements of the simplex table row, xN[j] are
741 * non-basic (auxiliary and structural) variables.
742 *
743 * RETURNS
744 *
745 * The routine returns number of non-zero elements in the simplex table
746 * row stored in the arrays ind and val.
747 *
748 * BACKGROUND
749 *
750 * The system of equality constraints of the LP problem is:
751 *
752 * xR = A * xS, (1)
753 *
754 * where xR is the vector of auxliary variables, xS is the vector of
755 * structural variables, A is the matrix of constraint coefficients.
756 *
757 * The system (1) can be written in homogenous form as follows:
758 *
759 * A~ * x = 0, (2)
760 *
761 * where A~ = (I | -A) is the augmented constraint matrix (has m rows
762 * and m+n columns), x = (xR | xS) is the vector of all (auxiliary and
763 * structural) variables.
764 *
765 * By definition for the current basis we have:
766 *
767 * A~ = (B | N), (3)
768 *
769 * where B is the basis matrix. Thus, the system (2) can be written as:
770 *
771 * B * xB + N * xN = 0. (4)
772 *
773 * From (4) it follows that:
774 *
775 * xB = A^ * xN, (5)
776 *
777 * where the matrix
778 *
779 * A^ = - inv(B) * N (6)
780 *
781 * is called the simplex table.
782 *
783 * It is understood that i-th row of the simplex table is:
784 *
785 * e * A^ = - e * inv(B) * N, (7)
786 *
787 * where e is a unity vector with e[i] = 1.
788 *
789 * To compute i-th row of the simplex table the routine first computes
790 * i-th row of the inverse:
791 *
792 * rho = inv(B') * e, (8)
793 *
794 * where B' is a matrix transposed to B, and then computes elements of
795 * i-th row of the simplex table as scalar products:
796 *
797 * alfa[i,j] = - rho * N[j] for all j, (9)
798 *
799 * where N[j] is a column of the augmented constraint matrix A~, which
800 * corresponds to some non-basic auxiliary or structural variable. */
801
glp_eval_tab_row(glp_prob * lp,int k,int ind[],double val[])802 int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[])
803 { int m = lp->m;
804 int n = lp->n;
805 int i, t, len, lll, *iii;
806 double alfa, *rho, *vvv;
807 if (!(m == 0 || lp->valid))
808 xerror("glp_eval_tab_row: basis factorization does not exist\n"
809 );
810 if (!(1 <= k && k <= m+n))
811 xerror("glp_eval_tab_row: k = %d; variable number out of range"
812 , k);
813 /* determine xB[i] which corresponds to x[k] */
814 if (k <= m)
815 i = glp_get_row_bind(lp, k);
816 else
817 i = glp_get_col_bind(lp, k-m);
818 if (i == 0)
819 xerror("glp_eval_tab_row: k = %d; variable must be basic", k);
820 xassert(1 <= i && i <= m);
821 /* allocate working arrays */
822 rho = xcalloc(1+m, sizeof(double));
823 iii = xcalloc(1+m, sizeof(int));
824 vvv = xcalloc(1+m, sizeof(double));
825 /* compute i-th row of the inverse; see (8) */
826 for (t = 1; t <= m; t++) rho[t] = 0.0;
827 rho[i] = 1.0;
828 glp_btran(lp, rho);
829 /* compute i-th row of the simplex table */
830 len = 0;
831 for (k = 1; k <= m+n; k++)
832 { if (k <= m)
833 { /* x[k] is auxiliary variable, so N[k] is a unity column */
834 if (glp_get_row_stat(lp, k) == GLP_BS) continue;
835 /* compute alfa[i,j]; see (9) */
836 alfa = - rho[k];
837 }
838 else
839 { /* x[k] is structural variable, so N[k] is a column of the
840 original constraint matrix A with negative sign */
841 if (glp_get_col_stat(lp, k-m) == GLP_BS) continue;
842 /* compute alfa[i,j]; see (9) */
843 lll = glp_get_mat_col(lp, k-m, iii, vvv);
844 alfa = 0.0;
845 for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t];
846 }
847 /* store alfa[i,j] */
848 if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa;
849 }
850 xassert(len <= n);
851 /* free working arrays */
852 xfree(rho);
853 xfree(iii);
854 xfree(vvv);
855 /* return to the calling program */
856 return len;
857 }
858
859 /***********************************************************************
860 * NAME
861 *
862 * glp_eval_tab_col - compute column of the simplex tableau
863 *
864 * SYNOPSIS
865 *
866 * int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]);
867 *
868 * DESCRIPTION
869 *
870 * The routine glp_eval_tab_col computes a column of the current simplex
871 * table for the non-basic variable, which is specified by the number k:
872 * if 1 <= k <= m, x[k] is k-th auxiliary variable; if m+1 <= k <= m+n,
873 * x[k] is (k-m)-th structural variable, where m is number of rows, and
874 * n is number of columns. The current basis must be available.
875 *
876 * The routine stores row indices and numerical values of non-zero
877 * elements of the computed column using sparse format to the locations
878 * ind[1], ..., ind[len] and val[1], ..., val[len] respectively, where
879 * 0 <= len <= m is number of non-zeros returned on exit.
880 *
881 * Element indices stored in the array ind have the same sense as the
882 * index k, i.e. indices 1 to m denote auxiliary variables and indices
883 * m+1 to m+n denote structural ones (all these variables are obviously
884 * basic by the definition).
885 *
886 * The computed column shows how basic variables depend on the specified
887 * non-basic variable x[k] = xN[j]:
888 *
889 * xB[1] = ... + alfa[1,j]*xN[j] + ...
890 * xB[2] = ... + alfa[2,j]*xN[j] + ...
891 * . . . . . .
892 * xB[m] = ... + alfa[m,j]*xN[j] + ...
893 *
894 * where alfa[i,j] are elements of the simplex table column, xB[i] are
895 * basic (auxiliary and structural) variables.
896 *
897 * RETURNS
898 *
899 * The routine returns number of non-zero elements in the simplex table
900 * column stored in the arrays ind and val.
901 *
902 * BACKGROUND
903 *
904 * As it was explained in comments to the routine glp_eval_tab_row (see
905 * above) the simplex table is the following matrix:
906 *
907 * A^ = - inv(B) * N. (1)
908 *
909 * Therefore j-th column of the simplex table is:
910 *
911 * A^ * e = - inv(B) * N * e = - inv(B) * N[j], (2)
912 *
913 * where e is a unity vector with e[j] = 1, B is the basis matrix, N[j]
914 * is a column of the augmented constraint matrix A~, which corresponds
915 * to the given non-basic auxiliary or structural variable. */
916
glp_eval_tab_col(glp_prob * lp,int k,int ind[],double val[])917 int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[])
918 { int m = lp->m;
919 int n = lp->n;
920 int t, len, stat;
921 double *col;
922 if (!(m == 0 || lp->valid))
923 xerror("glp_eval_tab_col: basis factorization does not exist\n"
924 );
925 if (!(1 <= k && k <= m+n))
926 xerror("glp_eval_tab_col: k = %d; variable number out of range"
927 , k);
928 if (k <= m)
929 stat = glp_get_row_stat(lp, k);
930 else
931 stat = glp_get_col_stat(lp, k-m);
932 if (stat == GLP_BS)
933 xerror("glp_eval_tab_col: k = %d; variable must be non-basic",
934 k);
935 /* obtain column N[k] with negative sign */
936 col = xcalloc(1+m, sizeof(double));
937 for (t = 1; t <= m; t++) col[t] = 0.0;
938 if (k <= m)
939 { /* x[k] is auxiliary variable, so N[k] is a unity column */
940 col[k] = -1.0;
941 }
942 else
943 { /* x[k] is structural variable, so N[k] is a column of the
944 original constraint matrix A with negative sign */
945 len = glp_get_mat_col(lp, k-m, ind, val);
946 for (t = 1; t <= len; t++) col[ind[t]] = val[t];
947 }
948 /* compute column of the simplex table, which corresponds to the
949 specified non-basic variable x[k] */
950 glp_ftran(lp, col);
951 len = 0;
952 for (t = 1; t <= m; t++)
953 { if (col[t] != 0.0)
954 { len++;
955 ind[len] = glp_get_bhead(lp, t);
956 val[len] = col[t];
957 }
958 }
959 xfree(col);
960 /* return to the calling program */
961 return len;
962 }
963
964 /***********************************************************************
965 * NAME
966 *
967 * glp_transform_row - transform explicitly specified row
968 *
969 * SYNOPSIS
970 *
971 * int glp_transform_row(glp_prob *P, int len, int ind[], double val[]);
972 *
973 * DESCRIPTION
974 *
975 * The routine glp_transform_row performs the same operation as the
976 * routine glp_eval_tab_row with exception that the row to be
977 * transformed is specified explicitly as a sparse vector.
978 *
979 * The explicitly specified row may be thought as a linear form:
980 *
981 * x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n], (1)
982 *
983 * where x is an auxiliary variable for this row, a[j] are coefficients
984 * of the linear form, x[m+j] are structural variables.
985 *
986 * On entry column indices and numerical values of non-zero elements of
987 * the row should be stored in locations ind[1], ..., ind[len] and
988 * val[1], ..., val[len], where len is the number of non-zero elements.
989 *
990 * This routine uses the system of equality constraints and the current
991 * basis in order to express the auxiliary variable x in (1) through the
992 * current non-basic variables (as if the transformed row were added to
993 * the problem object and its auxiliary variable were basic), i.e. the
994 * resultant row has the form:
995 *
996 * x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n], (2)
997 *
998 * where xN[j] are non-basic (auxiliary or structural) variables, n is
999 * the number of columns in the LP problem object.
1000 *
1001 * On exit the routine stores indices and numerical values of non-zero
1002 * elements of the resultant row (2) in locations ind[1], ..., ind[len']
1003 * and val[1], ..., val[len'], where 0 <= len' <= n is the number of
1004 * non-zero elements in the resultant row returned by the routine. Note
1005 * that indices (numbers) of non-basic variables stored in the array ind
1006 * correspond to original ordinal numbers of variables: indices 1 to m
1007 * mean auxiliary variables and indices m+1 to m+n mean structural ones.
1008 *
1009 * RETURNS
1010 *
1011 * The routine returns len', which is the number of non-zero elements in
1012 * the resultant row stored in the arrays ind and val.
1013 *
1014 * BACKGROUND
1015 *
1016 * The explicitly specified row (1) is transformed in the same way as it
1017 * were the objective function row.
1018 *
1019 * From (1) it follows that:
1020 *
1021 * x = aB * xB + aN * xN, (3)
1022 *
1023 * where xB is the vector of basic variables, xN is the vector of
1024 * non-basic variables.
1025 *
1026 * The simplex table, which corresponds to the current basis, is:
1027 *
1028 * xB = [-inv(B) * N] * xN. (4)
1029 *
1030 * Therefore substituting xB from (4) to (3) we have:
1031 *
1032 * x = aB * [-inv(B) * N] * xN + aN * xN =
1033 * (5)
1034 * = rho * (-N) * xN + aN * xN = alfa * xN,
1035 *
1036 * where:
1037 *
1038 * rho = inv(B') * aB, (6)
1039 *
1040 * and
1041 *
1042 * alfa = aN + rho * (-N) (7)
1043 *
1044 * is the resultant row computed by the routine. */
1045
glp_transform_row(glp_prob * P,int len,int ind[],double val[])1046 int glp_transform_row(glp_prob *P, int len, int ind[], double val[])
1047 { int i, j, k, m, n, t, lll, *iii;
1048 double alfa, *a, *aB, *rho, *vvv;
1049 if (!glp_bf_exists(P))
1050 xerror("glp_transform_row: basis factorization does not exist "
1051 "\n");
1052 m = glp_get_num_rows(P);
1053 n = glp_get_num_cols(P);
1054 /* unpack the row to be transformed to the array a */
1055 a = xcalloc(1+n, sizeof(double));
1056 for (j = 1; j <= n; j++) a[j] = 0.0;
1057 if (!(0 <= len && len <= n))
1058 xerror("glp_transform_row: len = %d; invalid row length\n",
1059 len);
1060 for (t = 1; t <= len; t++)
1061 { j = ind[t];
1062 if (!(1 <= j && j <= n))
1063 xerror("glp_transform_row: ind[%d] = %d; column index out o"
1064 "f range\n", t, j);
1065 if (val[t] == 0.0)
1066 xerror("glp_transform_row: val[%d] = 0; zero coefficient no"
1067 "t allowed\n", t);
1068 if (a[j] != 0.0)
1069 xerror("glp_transform_row: ind[%d] = %d; duplicate column i"
1070 "ndices not allowed\n", t, j);
1071 a[j] = val[t];
1072 }
1073 /* construct the vector aB */
1074 aB = xcalloc(1+m, sizeof(double));
1075 for (i = 1; i <= m; i++)
1076 { k = glp_get_bhead(P, i);
1077 /* xB[i] is k-th original variable */
1078 xassert(1 <= k && k <= m+n);
1079 aB[i] = (k <= m ? 0.0 : a[k-m]);
1080 }
1081 /* solve the system B'*rho = aB to compute the vector rho */
1082 rho = aB, glp_btran(P, rho);
1083 /* compute coefficients at non-basic auxiliary variables */
1084 len = 0;
1085 for (i = 1; i <= m; i++)
1086 { if (glp_get_row_stat(P, i) != GLP_BS)
1087 { alfa = - rho[i];
1088 if (alfa != 0.0)
1089 { len++;
1090 ind[len] = i;
1091 val[len] = alfa;
1092 }
1093 }
1094 }
1095 /* compute coefficients at non-basic structural variables */
1096 iii = xcalloc(1+m, sizeof(int));
1097 vvv = xcalloc(1+m, sizeof(double));
1098 for (j = 1; j <= n; j++)
1099 { if (glp_get_col_stat(P, j) != GLP_BS)
1100 { alfa = a[j];
1101 lll = glp_get_mat_col(P, j, iii, vvv);
1102 for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
1103 if (alfa != 0.0)
1104 { len++;
1105 ind[len] = m+j;
1106 val[len] = alfa;
1107 }
1108 }
1109 }
1110 xassert(len <= n);
1111 xfree(iii);
1112 xfree(vvv);
1113 xfree(aB);
1114 xfree(a);
1115 return len;
1116 }
1117
1118 /***********************************************************************
1119 * NAME
1120 *
1121 * glp_transform_col - transform explicitly specified column
1122 *
1123 * SYNOPSIS
1124 *
1125 * int glp_transform_col(glp_prob *P, int len, int ind[], double val[]);
1126 *
1127 * DESCRIPTION
1128 *
1129 * The routine glp_transform_col performs the same operation as the
1130 * routine glp_eval_tab_col with exception that the column to be
1131 * transformed is specified explicitly as a sparse vector.
1132 *
1133 * The explicitly specified column may be thought as if it were added
1134 * to the original system of equality constraints:
1135 *
1136 * x[1] = a[1,1]*x[m+1] + ... + a[1,n]*x[m+n] + a[1]*x
1137 * x[2] = a[2,1]*x[m+1] + ... + a[2,n]*x[m+n] + a[2]*x (1)
1138 * . . . . . . . . . . . . . . .
1139 * x[m] = a[m,1]*x[m+1] + ... + a[m,n]*x[m+n] + a[m]*x
1140 *
1141 * where x[i] are auxiliary variables, x[m+j] are structural variables,
1142 * x is a structural variable for the explicitly specified column, a[i]
1143 * are constraint coefficients for x.
1144 *
1145 * On entry row indices and numerical values of non-zero elements of
1146 * the column should be stored in locations ind[1], ..., ind[len] and
1147 * val[1], ..., val[len], where len is the number of non-zero elements.
1148 *
1149 * This routine uses the system of equality constraints and the current
1150 * basis in order to express the current basic variables through the
1151 * structural variable x in (1) (as if the transformed column were added
1152 * to the problem object and the variable x were non-basic), i.e. the
1153 * resultant column has the form:
1154 *
1155 * xB[1] = ... + alfa[1]*x
1156 * xB[2] = ... + alfa[2]*x (2)
1157 * . . . . . .
1158 * xB[m] = ... + alfa[m]*x
1159 *
1160 * where xB are basic (auxiliary and structural) variables, m is the
1161 * number of rows in the problem object.
1162 *
1163 * On exit the routine stores indices and numerical values of non-zero
1164 * elements of the resultant column (2) in locations ind[1], ...,
1165 * ind[len'] and val[1], ..., val[len'], where 0 <= len' <= m is the
1166 * number of non-zero element in the resultant column returned by the
1167 * routine. Note that indices (numbers) of basic variables stored in
1168 * the array ind correspond to original ordinal numbers of variables:
1169 * indices 1 to m mean auxiliary variables and indices m+1 to m+n mean
1170 * structural ones.
1171 *
1172 * RETURNS
1173 *
1174 * The routine returns len', which is the number of non-zero elements
1175 * in the resultant column stored in the arrays ind and val.
1176 *
1177 * BACKGROUND
1178 *
1179 * The explicitly specified column (1) is transformed in the same way
1180 * as any other column of the constraint matrix using the formula:
1181 *
1182 * alfa = inv(B) * a, (3)
1183 *
1184 * where alfa is the resultant column computed by the routine. */
1185
glp_transform_col(glp_prob * P,int len,int ind[],double val[])1186 int glp_transform_col(glp_prob *P, int len, int ind[], double val[])
1187 { int i, m, t;
1188 double *a, *alfa;
1189 if (!glp_bf_exists(P))
1190 xerror("glp_transform_col: basis factorization does not exist "
1191 "\n");
1192 m = glp_get_num_rows(P);
1193 /* unpack the column to be transformed to the array a */
1194 a = xcalloc(1+m, sizeof(double));
1195 for (i = 1; i <= m; i++) a[i] = 0.0;
1196 if (!(0 <= len && len <= m))
1197 xerror("glp_transform_col: len = %d; invalid column length\n",
1198 len);
1199 for (t = 1; t <= len; t++)
1200 { i = ind[t];
1201 if (!(1 <= i && i <= m))
1202 xerror("glp_transform_col: ind[%d] = %d; row index out of r"
1203 "ange\n", t, i);
1204 if (val[t] == 0.0)
1205 xerror("glp_transform_col: val[%d] = 0; zero coefficient no"
1206 "t allowed\n", t);
1207 if (a[i] != 0.0)
1208 xerror("glp_transform_col: ind[%d] = %d; duplicate row indi"
1209 "ces not allowed\n", t, i);
1210 a[i] = val[t];
1211 }
1212 /* solve the system B*a = alfa to compute the vector alfa */
1213 alfa = a, glp_ftran(P, alfa);
1214 /* store resultant coefficients */
1215 len = 0;
1216 for (i = 1; i <= m; i++)
1217 { if (alfa[i] != 0.0)
1218 { len++;
1219 ind[len] = glp_get_bhead(P, i);
1220 val[len] = alfa[i];
1221 }
1222 }
1223 xfree(a);
1224 return len;
1225 }
1226
1227 /***********************************************************************
1228 * NAME
1229 *
1230 * glp_prim_rtest - perform primal ratio test
1231 *
1232 * SYNOPSIS
1233 *
1234 * int glp_prim_rtest(glp_prob *P, int len, const int ind[],
1235 * const double val[], int dir, double eps);
1236 *
1237 * DESCRIPTION
1238 *
1239 * The routine glp_prim_rtest performs the primal ratio test using an
1240 * explicitly specified column of the simplex table.
1241 *
1242 * The current basic solution associated with the LP problem object
1243 * must be primal feasible.
1244 *
1245 * The explicitly specified column of the simplex table shows how the
1246 * basic variables xB depend on some non-basic variable x (which is not
1247 * necessarily presented in the problem object):
1248 *
1249 * xB[1] = ... + alfa[1] * x + ...
1250 * xB[2] = ... + alfa[2] * x + ... (*)
1251 * . . . . . . . .
1252 * xB[m] = ... + alfa[m] * x + ...
1253 *
1254 * The column (*) is specifed on entry to the routine using the sparse
1255 * format. Ordinal numbers of basic variables xB[i] should be placed in
1256 * locations ind[1], ..., ind[len], where ordinal number 1 to m denote
1257 * auxiliary variables, and ordinal numbers m+1 to m+n denote structural
1258 * variables. The corresponding non-zero coefficients alfa[i] should be
1259 * placed in locations val[1], ..., val[len]. The arrays ind and val are
1260 * not changed on exit.
1261 *
1262 * The parameter dir specifies direction in which the variable x changes
1263 * on entering the basis: +1 means increasing, -1 means decreasing.
1264 *
1265 * The parameter eps is an absolute tolerance (small positive number)
1266 * used by the routine to skip small alfa[j] of the row (*).
1267 *
1268 * The routine determines which basic variable (among specified in
1269 * ind[1], ..., ind[len]) should leave the basis in order to keep primal
1270 * feasibility.
1271 *
1272 * RETURNS
1273 *
1274 * The routine glp_prim_rtest returns the index piv in the arrays ind
1275 * and val corresponding to the pivot element chosen, 1 <= piv <= len.
1276 * If the adjacent basic solution is primal unbounded and therefore the
1277 * choice cannot be made, the routine returns zero.
1278 *
1279 * COMMENTS
1280 *
1281 * If the non-basic variable x is presented in the LP problem object,
1282 * the column (*) can be computed with the routine glp_eval_tab_col;
1283 * otherwise it can be computed with the routine glp_transform_col. */
1284
glp_prim_rtest(glp_prob * P,int len,const int ind[],const double val[],int dir,double eps)1285 int glp_prim_rtest(glp_prob *P, int len, const int ind[],
1286 const double val[], int dir, double eps)
1287 { int k, m, n, piv, t, type, stat;
1288 double alfa, big, beta, lb, ub, temp, teta;
1289 if (glp_get_prim_stat(P) != GLP_FEAS)
1290 xerror("glp_prim_rtest: basic solution is not primal feasible "
1291 "\n");
1292 if (!(dir == +1 || dir == -1))
1293 xerror("glp_prim_rtest: dir = %d; invalid parameter\n", dir);
1294 if (!(0.0 < eps && eps < 1.0))
1295 xerror("glp_prim_rtest: eps = %g; invalid parameter\n", eps);
1296 m = glp_get_num_rows(P);
1297 n = glp_get_num_cols(P);
1298 /* initial settings */
1299 piv = 0, teta = DBL_MAX, big = 0.0;
1300 /* walk through the entries of the specified column */
1301 for (t = 1; t <= len; t++)
1302 { /* get the ordinal number of basic variable */
1303 k = ind[t];
1304 if (!(1 <= k && k <= m+n))
1305 xerror("glp_prim_rtest: ind[%d] = %d; variable number out o"
1306 "f range\n", t, k);
1307 /* determine type, bounds, status and primal value of basic
1308 variable xB[i] = x[k] in the current basic solution */
1309 if (k <= m)
1310 { type = glp_get_row_type(P, k);
1311 lb = glp_get_row_lb(P, k);
1312 ub = glp_get_row_ub(P, k);
1313 stat = glp_get_row_stat(P, k);
1314 beta = glp_get_row_prim(P, k);
1315 }
1316 else
1317 { type = glp_get_col_type(P, k-m);
1318 lb = glp_get_col_lb(P, k-m);
1319 ub = glp_get_col_ub(P, k-m);
1320 stat = glp_get_col_stat(P, k-m);
1321 beta = glp_get_col_prim(P, k-m);
1322 }
1323 if (stat != GLP_BS)
1324 xerror("glp_prim_rtest: ind[%d] = %d; non-basic variable no"
1325 "t allowed\n", t, k);
1326 /* determine influence coefficient at basic variable xB[i]
1327 in the explicitly specified column and turn to the case of
1328 increasing the variable x in order to simplify the program
1329 logic */
1330 alfa = (dir > 0 ? + val[t] : - val[t]);
1331 /* analyze main cases */
1332 if (type == GLP_FR)
1333 { /* xB[i] is free variable */
1334 continue;
1335 }
1336 else if (type == GLP_LO)
1337 lo: { /* xB[i] has an lower bound */
1338 if (alfa > - eps) continue;
1339 temp = (lb - beta) / alfa;
1340 }
1341 else if (type == GLP_UP)
1342 up: { /* xB[i] has an upper bound */
1343 if (alfa < + eps) continue;
1344 temp = (ub - beta) / alfa;
1345 }
1346 else if (type == GLP_DB)
1347 { /* xB[i] has both lower and upper bounds */
1348 if (alfa < 0.0) goto lo; else goto up;
1349 }
1350 else if (type == GLP_FX)
1351 { /* xB[i] is fixed variable */
1352 if (- eps < alfa && alfa < + eps) continue;
1353 temp = 0.0;
1354 }
1355 else
1356 xassert(type != type);
1357 /* if the value of the variable xB[i] violates its lower or
1358 upper bound (slightly, because the current basis is assumed
1359 to be primal feasible), temp is negative; we can think this
1360 happens due to round-off errors and the value is exactly on
1361 the bound; this allows replacing temp by zero */
1362 if (temp < 0.0) temp = 0.0;
1363 /* apply the minimal ratio test */
1364 if (teta > temp || teta == temp && big < fabs(alfa))
1365 piv = t, teta = temp, big = fabs(alfa);
1366 }
1367 /* return index of the pivot element chosen */
1368 return piv;
1369 }
1370
1371 /***********************************************************************
1372 * NAME
1373 *
1374 * glp_dual_rtest - perform dual ratio test
1375 *
1376 * SYNOPSIS
1377 *
1378 * int glp_dual_rtest(glp_prob *P, int len, const int ind[],
1379 * const double val[], int dir, double eps);
1380 *
1381 * DESCRIPTION
1382 *
1383 * The routine glp_dual_rtest performs the dual ratio test using an
1384 * explicitly specified row of the simplex table.
1385 *
1386 * The current basic solution associated with the LP problem object
1387 * must be dual feasible.
1388 *
1389 * The explicitly specified row of the simplex table is a linear form
1390 * that shows how some basic variable x (which is not necessarily
1391 * presented in the problem object) depends on non-basic variables xN:
1392 *
1393 * x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*)
1394 *
1395 * The row (*) is specified on entry to the routine using the sparse
1396 * format. Ordinal numbers of non-basic variables xN[j] should be placed
1397 * in locations ind[1], ..., ind[len], where ordinal numbers 1 to m
1398 * denote auxiliary variables, and ordinal numbers m+1 to m+n denote
1399 * structural variables. The corresponding non-zero coefficients alfa[j]
1400 * should be placed in locations val[1], ..., val[len]. The arrays ind
1401 * and val are not changed on exit.
1402 *
1403 * The parameter dir specifies direction in which the variable x changes
1404 * on leaving the basis: +1 means that x goes to its lower bound, and -1
1405 * means that x goes to its upper bound.
1406 *
1407 * The parameter eps is an absolute tolerance (small positive number)
1408 * used by the routine to skip small alfa[j] of the row (*).
1409 *
1410 * The routine determines which non-basic variable (among specified in
1411 * ind[1], ..., ind[len]) should enter the basis in order to keep dual
1412 * feasibility.
1413 *
1414 * RETURNS
1415 *
1416 * The routine glp_dual_rtest returns the index piv in the arrays ind
1417 * and val corresponding to the pivot element chosen, 1 <= piv <= len.
1418 * If the adjacent basic solution is dual unbounded and therefore the
1419 * choice cannot be made, the routine returns zero.
1420 *
1421 * COMMENTS
1422 *
1423 * If the basic variable x is presented in the LP problem object, the
1424 * row (*) can be computed with the routine glp_eval_tab_row; otherwise
1425 * it can be computed with the routine glp_transform_row. */
1426
glp_dual_rtest(glp_prob * P,int len,const int ind[],const double val[],int dir,double eps)1427 int glp_dual_rtest(glp_prob *P, int len, const int ind[],
1428 const double val[], int dir, double eps)
1429 { int k, m, n, piv, t, stat;
1430 double alfa, big, cost, obj, temp, teta;
1431 if (glp_get_dual_stat(P) != GLP_FEAS)
1432 xerror("glp_dual_rtest: basic solution is not dual feasible\n")
1433 ;
1434 if (!(dir == +1 || dir == -1))
1435 xerror("glp_dual_rtest: dir = %d; invalid parameter\n", dir);
1436 if (!(0.0 < eps && eps < 1.0))
1437 xerror("glp_dual_rtest: eps = %g; invalid parameter\n", eps);
1438 m = glp_get_num_rows(P);
1439 n = glp_get_num_cols(P);
1440 /* take into account optimization direction */
1441 obj = (glp_get_obj_dir(P) == GLP_MIN ? +1.0 : -1.0);
1442 /* initial settings */
1443 piv = 0, teta = DBL_MAX, big = 0.0;
1444 /* walk through the entries of the specified row */
1445 for (t = 1; t <= len; t++)
1446 { /* get ordinal number of non-basic variable */
1447 k = ind[t];
1448 if (!(1 <= k && k <= m+n))
1449 xerror("glp_dual_rtest: ind[%d] = %d; variable number out o"
1450 "f range\n", t, k);
1451 /* determine status and reduced cost of non-basic variable
1452 x[k] = xN[j] in the current basic solution */
1453 if (k <= m)
1454 { stat = glp_get_row_stat(P, k);
1455 cost = glp_get_row_dual(P, k);
1456 }
1457 else
1458 { stat = glp_get_col_stat(P, k-m);
1459 cost = glp_get_col_dual(P, k-m);
1460 }
1461 if (stat == GLP_BS)
1462 xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al"
1463 "lowed\n", t, k);
1464 /* determine influence coefficient at non-basic variable xN[j]
1465 in the explicitly specified row and turn to the case of
1466 increasing the variable x in order to simplify the program
1467 logic */
1468 alfa = (dir > 0 ? + val[t] : - val[t]);
1469 /* analyze main cases */
1470 if (stat == GLP_NL)
1471 { /* xN[j] is on its lower bound */
1472 if (alfa < + eps) continue;
1473 temp = (obj * cost) / alfa;
1474 }
1475 else if (stat == GLP_NU)
1476 { /* xN[j] is on its upper bound */
1477 if (alfa > - eps) continue;
1478 temp = (obj * cost) / alfa;
1479 }
1480 else if (stat == GLP_NF)
1481 { /* xN[j] is non-basic free variable */
1482 if (- eps < alfa && alfa < + eps) continue;
1483 temp = 0.0;
1484 }
1485 else if (stat == GLP_NS)
1486 { /* xN[j] is non-basic fixed variable */
1487 continue;
1488 }
1489 else
1490 xassert(stat != stat);
1491 /* if the reduced cost of the variable xN[j] violates its zero
1492 bound (slightly, because the current basis is assumed to be
1493 dual feasible), temp is negative; we can think this happens
1494 due to round-off errors and the reduced cost is exact zero;
1495 this allows replacing temp by zero */
1496 if (temp < 0.0) temp = 0.0;
1497 /* apply the minimal ratio test */
1498 if (teta > temp || teta == temp && big < fabs(alfa))
1499 piv = t, teta = temp, big = fabs(alfa);
1500 }
1501 /* return index of the pivot element chosen */
1502 return piv;
1503 }
1504
1505 /***********************************************************************
1506 * NAME
1507 *
1508 * glp_analyze_row - simulate one iteration of dual simplex method
1509 *
1510 * SYNOPSIS
1511 *
1512 * int glp_analyze_row(glp_prob *P, int len, const int ind[],
1513 * const double val[], int type, double rhs, double eps, int *piv,
1514 * double *x, double *dx, double *y, double *dy, double *dz);
1515 *
1516 * DESCRIPTION
1517 *
1518 * Let the current basis be optimal or dual feasible, and there be
1519 * specified a row (constraint), which is violated by the current basic
1520 * solution. The routine glp_analyze_row simulates one iteration of the
1521 * dual simplex method to determine some information on the adjacent
1522 * basis (see below), where the specified row becomes active constraint
1523 * (i.e. its auxiliary variable becomes non-basic).
1524 *
1525 * The current basic solution associated with the problem object passed
1526 * to the routine must be dual feasible, and its primal components must
1527 * be defined.
1528 *
1529 * The row to be analyzed must be previously transformed either with
1530 * the routine glp_eval_tab_row (if the row is in the problem object)
1531 * or with the routine glp_transform_row (if the row is external, i.e.
1532 * not in the problem object). This is needed to express the row only
1533 * through (auxiliary and structural) variables, which are non-basic in
1534 * the current basis:
1535 *
1536 * y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n],
1537 *
1538 * where y is an auxiliary variable of the row, alfa[j] is an influence
1539 * coefficient, xN[j] is a non-basic variable.
1540 *
1541 * The row is passed to the routine in sparse format. Ordinal numbers
1542 * of non-basic variables are stored in locations ind[1], ..., ind[len],
1543 * where numbers 1 to m denote auxiliary variables while numbers m+1 to
1544 * m+n denote structural variables. Corresponding non-zero coefficients
1545 * alfa[j] are stored in locations val[1], ..., val[len]. The arrays
1546 * ind and val are ot changed on exit.
1547 *
1548 * The parameters type and rhs specify the row type and its right-hand
1549 * side as follows:
1550 *
1551 * type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs
1552 *
1553 * type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs
1554 *
1555 * The parameter eps is an absolute tolerance (small positive number)
1556 * used by the routine to skip small coefficients alfa[j] on performing
1557 * the dual ratio test.
1558 *
1559 * If the operation was successful, the routine stores the following
1560 * information to corresponding location (if some parameter is NULL,
1561 * its value is not stored):
1562 *
1563 * piv index in the array ind and val, 1 <= piv <= len, determining
1564 * the non-basic variable, which would enter the adjacent basis;
1565 *
1566 * x value of the non-basic variable in the current basis;
1567 *
1568 * dx difference between values of the non-basic variable in the
1569 * adjacent and current bases, dx = x.new - x.old;
1570 *
1571 * y value of the row (i.e. of its auxiliary variable) in the
1572 * current basis;
1573 *
1574 * dy difference between values of the row in the adjacent and
1575 * current bases, dy = y.new - y.old;
1576 *
1577 * dz difference between values of the objective function in the
1578 * adjacent and current bases, dz = z.new - z.old. Note that in
1579 * case of minimization dz >= 0, and in case of maximization
1580 * dz <= 0, i.e. in the adjacent basis the objective function
1581 * always gets worse (degrades). */
1582
_glp_analyze_row(glp_prob * P,int len,const int ind[],const double val[],int type,double rhs,double eps,int * _piv,double * _x,double * _dx,double * _y,double * _dy,double * _dz)1583 int _glp_analyze_row(glp_prob *P, int len, const int ind[],
1584 const double val[], int type, double rhs, double eps, int *_piv,
1585 double *_x, double *_dx, double *_y, double *_dy, double *_dz)
1586 { int t, k, dir, piv, ret = 0;
1587 double x, dx, y, dy, dz;
1588 if (P->pbs_stat == GLP_UNDEF)
1589 xerror("glp_analyze_row: primal basic solution components are "
1590 "undefined\n");
1591 if (P->dbs_stat != GLP_FEAS)
1592 xerror("glp_analyze_row: basic solution is not dual feasible\n"
1593 );
1594 /* compute the row value y = sum alfa[j] * xN[j] in the current
1595 basis */
1596 if (!(0 <= len && len <= P->n))
1597 xerror("glp_analyze_row: len = %d; invalid row length\n", len);
1598 y = 0.0;
1599 for (t = 1; t <= len; t++)
1600 { /* determine value of x[k] = xN[j] in the current basis */
1601 k = ind[t];
1602 if (!(1 <= k && k <= P->m+P->n))
1603 xerror("glp_analyze_row: ind[%d] = %d; row/column index out"
1604 " of range\n", t, k);
1605 if (k <= P->m)
1606 { /* x[k] is auxiliary variable */
1607 if (P->row[k]->stat == GLP_BS)
1608 xerror("glp_analyze_row: ind[%d] = %d; basic auxiliary v"
1609 "ariable is not allowed\n", t, k);
1610 x = P->row[k]->prim;
1611 }
1612 else
1613 { /* x[k] is structural variable */
1614 if (P->col[k-P->m]->stat == GLP_BS)
1615 xerror("glp_analyze_row: ind[%d] = %d; basic structural "
1616 "variable is not allowed\n", t, k);
1617 x = P->col[k-P->m]->prim;
1618 }
1619 y += val[t] * x;
1620 }
1621 /* check if the row is primal infeasible in the current basis,
1622 i.e. the constraint is violated at the current point */
1623 if (type == GLP_LO)
1624 { if (y >= rhs)
1625 { /* the constraint is not violated */
1626 ret = 1;
1627 goto done;
1628 }
1629 /* in the adjacent basis y goes to its lower bound */
1630 dir = +1;
1631 }
1632 else if (type == GLP_UP)
1633 { if (y <= rhs)
1634 { /* the constraint is not violated */
1635 ret = 1;
1636 goto done;
1637 }
1638 /* in the adjacent basis y goes to its upper bound */
1639 dir = -1;
1640 }
1641 else
1642 xerror("glp_analyze_row: type = %d; invalid parameter\n",
1643 type);
1644 /* compute dy = y.new - y.old */
1645 dy = rhs - y;
1646 /* perform dual ratio test to determine which non-basic variable
1647 should enter the adjacent basis to keep it dual feasible */
1648 piv = glp_dual_rtest(P, len, ind, val, dir, eps);
1649 if (piv == 0)
1650 { /* no dual feasible adjacent basis exists */
1651 ret = 2;
1652 goto done;
1653 }
1654 /* non-basic variable x[k] = xN[j] should enter the basis */
1655 k = ind[piv];
1656 xassert(1 <= k && k <= P->m+P->n);
1657 /* determine its value in the current basis */
1658 if (k <= P->m)
1659 x = P->row[k]->prim;
1660 else
1661 x = P->col[k-P->m]->prim;
1662 /* compute dx = x.new - x.old = dy / alfa[j] */
1663 xassert(val[piv] != 0.0);
1664 dx = dy / val[piv];
1665 /* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced
1666 cost of xN[j] in the current basis */
1667 if (k <= P->m)
1668 dz = P->row[k]->dual * dx;
1669 else
1670 dz = P->col[k-P->m]->dual * dx;
1671 /* store the analysis results */
1672 if (_piv != NULL) *_piv = piv;
1673 if (_x != NULL) *_x = x;
1674 if (_dx != NULL) *_dx = dx;
1675 if (_y != NULL) *_y = y;
1676 if (_dy != NULL) *_dy = dy;
1677 if (_dz != NULL) *_dz = dz;
1678 done: return ret;
1679 }
1680
1681 #if 0
1682 int main(void)
1683 { /* example program for the routine glp_analyze_row */
1684 glp_prob *P;
1685 glp_smcp parm;
1686 int i, k, len, piv, ret, ind[1+100];
1687 double rhs, x, dx, y, dy, dz, val[1+100];
1688 P = glp_create_prob();
1689 /* read plan.mps (see glpk/examples) */
1690 ret = glp_read_mps(P, GLP_MPS_DECK, NULL, "plan.mps");
1691 glp_assert(ret == 0);
1692 /* and solve it to optimality */
1693 ret = glp_simplex(P, NULL);
1694 glp_assert(ret == 0);
1695 glp_assert(glp_get_status(P) == GLP_OPT);
1696 /* the optimal objective value is 296.217 */
1697 /* we would like to know what happens if we would add a new row
1698 (constraint) to plan.mps:
1699 .01 * bin1 + .01 * bin2 + .02 * bin4 + .02 * bin5 <= 12 */
1700 /* first, we specify this new row */
1701 glp_create_index(P);
1702 len = 0;
1703 ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
1704 ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
1705 ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
1706 ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
1707 rhs = 12;
1708 /* then we can compute value of the row (i.e. of its auxiliary
1709 variable) in the current basis to see if the constraint is
1710 violated */
1711 y = 0.0;
1712 for (k = 1; k <= len; k++)
1713 y += val[k] * glp_get_col_prim(P, ind[k]);
1714 glp_printf("y = %g\n", y);
1715 /* this prints y = 15.1372, so the constraint is violated, since
1716 we require that y <= rhs = 12 */
1717 /* now we transform the row to express it only through non-basic
1718 (auxiliary and artificial) variables */
1719 len = glp_transform_row(P, len, ind, val);
1720 /* finally, we simulate one step of the dual simplex method to
1721 obtain necessary information for the adjacent basis */
1722 ret = _glp_analyze_row(P, len, ind, val, GLP_UP, rhs, 1e-9, &piv,
1723 &x, &dx, &y, &dy, &dz);
1724 glp_assert(ret == 0);
1725 glp_printf("k = %d, x = %g; dx = %g; y = %g; dy = %g; dz = %g\n",
1726 ind[piv], x, dx, y, dy, dz);
1727 /* this prints dz = 5.64418 and means that in the adjacent basis
1728 the objective function would be 296.217 + 5.64418 = 301.861 */
1729 /* now we actually include the row into the problem object; note
1730 that the arrays ind and val are clobbered, so we need to build
1731 them once again */
1732 len = 0;
1733 ind[++len] = glp_find_col(P, "BIN1"), val[len] = .01;
1734 ind[++len] = glp_find_col(P, "BIN2"), val[len] = .01;
1735 ind[++len] = glp_find_col(P, "BIN4"), val[len] = .02;
1736 ind[++len] = glp_find_col(P, "BIN5"), val[len] = .02;
1737 rhs = 12;
1738 i = glp_add_rows(P, 1);
1739 glp_set_row_bnds(P, i, GLP_UP, 0, rhs);
1740 glp_set_mat_row(P, i, len, ind, val);
1741 /* and perform one dual simplex iteration */
1742 glp_init_smcp(&parm);
1743 parm.meth = GLP_DUAL;
1744 parm.it_lim = 1;
1745 glp_simplex(P, &parm);
1746 /* the current objective value is 301.861 */
1747 return 0;
1748 }
1749 #endif
1750
1751 /***********************************************************************
1752 * NAME
1753 *
1754 * glp_analyze_bound - analyze active bound of non-basic variable
1755 *
1756 * SYNOPSIS
1757 *
1758 * void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1,
1759 * double *limit2, int *var2);
1760 *
1761 * DESCRIPTION
1762 *
1763 * The routine glp_analyze_bound analyzes the effect of varying the
1764 * active bound of specified non-basic variable.
1765 *
1766 * The non-basic variable is specified by the parameter k, where
1767 * 1 <= k <= m means auxiliary variable of corresponding row while
1768 * m+1 <= k <= m+n means structural variable (column).
1769 *
1770 * Note that the current basic solution must be optimal, and the basis
1771 * factorization must exist.
1772 *
1773 * Results of the analysis have the following meaning.
1774 *
1775 * value1 is the minimal value of the active bound, at which the basis
1776 * still remains primal feasible and thus optimal. -DBL_MAX means that
1777 * the active bound has no lower limit.
1778 *
1779 * var1 is the ordinal number of an auxiliary (1 to m) or structural
1780 * (m+1 to n) basic variable, which reaches its bound first and thereby
1781 * limits further decreasing the active bound being analyzed.
1782 * if value1 = -DBL_MAX, var1 is set to 0.
1783 *
1784 * value2 is the maximal value of the active bound, at which the basis
1785 * still remains primal feasible and thus optimal. +DBL_MAX means that
1786 * the active bound has no upper limit.
1787 *
1788 * var2 is the ordinal number of an auxiliary (1 to m) or structural
1789 * (m+1 to n) basic variable, which reaches its bound first and thereby
1790 * limits further increasing the active bound being analyzed.
1791 * if value2 = +DBL_MAX, var2 is set to 0. */
1792
glp_analyze_bound(glp_prob * P,int k,double * value1,int * var1,double * value2,int * var2)1793 void glp_analyze_bound(glp_prob *P, int k, double *value1, int *var1,
1794 double *value2, int *var2)
1795 { GLPROW *row;
1796 GLPCOL *col;
1797 int m, n, stat, kase, p, len, piv, *ind;
1798 double x, new_x, ll, uu, xx, delta, *val;
1799 #if 0 /* 04/IV-2016 */
1800 /* sanity checks */
1801 if (P == NULL || P->magic != GLP_PROB_MAGIC)
1802 xerror("glp_analyze_bound: P = %p; invalid problem object\n",
1803 P);
1804 #endif
1805 m = P->m, n = P->n;
1806 if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
1807 xerror("glp_analyze_bound: optimal basic solution required\n");
1808 if (!(m == 0 || P->valid))
1809 xerror("glp_analyze_bound: basis factorization required\n");
1810 if (!(1 <= k && k <= m+n))
1811 xerror("glp_analyze_bound: k = %d; variable number out of rang"
1812 "e\n", k);
1813 /* retrieve information about the specified non-basic variable
1814 x[k] whose active bound is to be analyzed */
1815 if (k <= m)
1816 { row = P->row[k];
1817 stat = row->stat;
1818 x = row->prim;
1819 }
1820 else
1821 { col = P->col[k-m];
1822 stat = col->stat;
1823 x = col->prim;
1824 }
1825 if (stat == GLP_BS)
1826 xerror("glp_analyze_bound: k = %d; basic variable not allowed "
1827 "\n", k);
1828 /* allocate working arrays */
1829 ind = xcalloc(1+m, sizeof(int));
1830 val = xcalloc(1+m, sizeof(double));
1831 /* compute column of the simplex table corresponding to the
1832 non-basic variable x[k] */
1833 len = glp_eval_tab_col(P, k, ind, val);
1834 xassert(0 <= len && len <= m);
1835 /* perform analysis */
1836 for (kase = -1; kase <= +1; kase += 2)
1837 { /* kase < 0 means active bound of x[k] is decreasing;
1838 kase > 0 means active bound of x[k] is increasing */
1839 /* use the primal ratio test to determine some basic variable
1840 x[p] which reaches its bound first */
1841 piv = glp_prim_rtest(P, len, ind, val, kase, 1e-9);
1842 if (piv == 0)
1843 { /* nothing limits changing the active bound of x[k] */
1844 p = 0;
1845 new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX);
1846 goto store;
1847 }
1848 /* basic variable x[p] limits changing the active bound of
1849 x[k]; determine its value in the current basis */
1850 xassert(1 <= piv && piv <= len);
1851 p = ind[piv];
1852 if (p <= m)
1853 { row = P->row[p];
1854 ll = glp_get_row_lb(P, row->i);
1855 uu = glp_get_row_ub(P, row->i);
1856 stat = row->stat;
1857 xx = row->prim;
1858 }
1859 else
1860 { col = P->col[p-m];
1861 ll = glp_get_col_lb(P, col->j);
1862 uu = glp_get_col_ub(P, col->j);
1863 stat = col->stat;
1864 xx = col->prim;
1865 }
1866 xassert(stat == GLP_BS);
1867 /* determine delta x[p] = bound of x[p] - value of x[p] */
1868 if (kase < 0 && val[piv] > 0.0 ||
1869 kase > 0 && val[piv] < 0.0)
1870 { /* delta x[p] < 0, so x[p] goes toward its lower bound */
1871 xassert(ll != -DBL_MAX);
1872 delta = ll - xx;
1873 }
1874 else
1875 { /* delta x[p] > 0, so x[p] goes toward its upper bound */
1876 xassert(uu != +DBL_MAX);
1877 delta = uu - xx;
1878 }
1879 /* delta x[p] = alfa[p,k] * delta x[k], so new x[k] = x[k] +
1880 delta x[k] = x[k] + delta x[p] / alfa[p,k] is the value of
1881 x[k] in the adjacent basis */
1882 xassert(val[piv] != 0.0);
1883 new_x = x + delta / val[piv];
1884 store: /* store analysis results */
1885 if (kase < 0)
1886 { if (value1 != NULL) *value1 = new_x;
1887 if (var1 != NULL) *var1 = p;
1888 }
1889 else
1890 { if (value2 != NULL) *value2 = new_x;
1891 if (var2 != NULL) *var2 = p;
1892 }
1893 }
1894 /* free working arrays */
1895 xfree(ind);
1896 xfree(val);
1897 return;
1898 }
1899
1900 /***********************************************************************
1901 * NAME
1902 *
1903 * glp_analyze_coef - analyze objective coefficient at basic variable
1904 *
1905 * SYNOPSIS
1906 *
1907 * void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
1908 * double *value1, double *coef2, int *var2, double *value2);
1909 *
1910 * DESCRIPTION
1911 *
1912 * The routine glp_analyze_coef analyzes the effect of varying the
1913 * objective coefficient at specified basic variable.
1914 *
1915 * The basic variable is specified by the parameter k, where
1916 * 1 <= k <= m means auxiliary variable of corresponding row while
1917 * m+1 <= k <= m+n means structural variable (column).
1918 *
1919 * Note that the current basic solution must be optimal, and the basis
1920 * factorization must exist.
1921 *
1922 * Results of the analysis have the following meaning.
1923 *
1924 * coef1 is the minimal value of the objective coefficient, at which
1925 * the basis still remains dual feasible and thus optimal. -DBL_MAX
1926 * means that the objective coefficient has no lower limit.
1927 *
1928 * var1 is the ordinal number of an auxiliary (1 to m) or structural
1929 * (m+1 to n) non-basic variable, whose reduced cost reaches its zero
1930 * bound first and thereby limits further decreasing the objective
1931 * coefficient being analyzed. If coef1 = -DBL_MAX, var1 is set to 0.
1932 *
1933 * value1 is value of the basic variable being analyzed in an adjacent
1934 * basis, which is defined as follows. Let the objective coefficient
1935 * reaches its minimal value (coef1) and continues decreasing. Then the
1936 * reduced cost of the limiting non-basic variable (var1) becomes dual
1937 * infeasible and the current basis becomes non-optimal that forces the
1938 * limiting non-basic variable to enter the basis replacing there some
1939 * basic variable that leaves the basis to keep primal feasibility.
1940 * Should note that on determining the adjacent basis current bounds
1941 * of the basic variable being analyzed are ignored as if it were free
1942 * (unbounded) variable, so it cannot leave the basis. It may happen
1943 * that no dual feasible adjacent basis exists, in which case value1 is
1944 * set to -DBL_MAX or +DBL_MAX.
1945 *
1946 * coef2 is the maximal value of the objective coefficient, at which
1947 * the basis still remains dual feasible and thus optimal. +DBL_MAX
1948 * means that the objective coefficient has no upper limit.
1949 *
1950 * var2 is the ordinal number of an auxiliary (1 to m) or structural
1951 * (m+1 to n) non-basic variable, whose reduced cost reaches its zero
1952 * bound first and thereby limits further increasing the objective
1953 * coefficient being analyzed. If coef2 = +DBL_MAX, var2 is set to 0.
1954 *
1955 * value2 is value of the basic variable being analyzed in an adjacent
1956 * basis, which is defined exactly in the same way as value1 above with
1957 * exception that now the objective coefficient is increasing. */
1958
glp_analyze_coef(glp_prob * P,int k,double * coef1,int * var1,double * value1,double * coef2,int * var2,double * value2)1959 void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
1960 double *value1, double *coef2, int *var2, double *value2)
1961 { GLPROW *row; GLPCOL *col;
1962 int m, n, type, stat, kase, p, q, dir, clen, cpiv, rlen, rpiv,
1963 *cind, *rind;
1964 double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx,
1965 *rval, *cval;
1966 #if 0 /* 04/IV-2016 */
1967 /* sanity checks */
1968 if (P == NULL || P->magic != GLP_PROB_MAGIC)
1969 xerror("glp_analyze_coef: P = %p; invalid problem object\n",
1970 P);
1971 #endif
1972 m = P->m, n = P->n;
1973 if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
1974 xerror("glp_analyze_coef: optimal basic solution required\n");
1975 if (!(m == 0 || P->valid))
1976 xerror("glp_analyze_coef: basis factorization required\n");
1977 if (!(1 <= k && k <= m+n))
1978 xerror("glp_analyze_coef: k = %d; variable number out of range"
1979 "\n", k);
1980 /* retrieve information about the specified basic variable x[k]
1981 whose objective coefficient c[k] is to be analyzed */
1982 if (k <= m)
1983 { row = P->row[k];
1984 type = row->type;
1985 lb = row->lb;
1986 ub = row->ub;
1987 coef = 0.0;
1988 stat = row->stat;
1989 x = row->prim;
1990 }
1991 else
1992 { col = P->col[k-m];
1993 type = col->type;
1994 lb = col->lb;
1995 ub = col->ub;
1996 coef = col->coef;
1997 stat = col->stat;
1998 x = col->prim;
1999 }
2000 if (stat != GLP_BS)
2001 xerror("glp_analyze_coef: k = %d; non-basic variable not allow"
2002 "ed\n", k);
2003 /* allocate working arrays */
2004 cind = xcalloc(1+m, sizeof(int));
2005 cval = xcalloc(1+m, sizeof(double));
2006 rind = xcalloc(1+n, sizeof(int));
2007 rval = xcalloc(1+n, sizeof(double));
2008 /* compute row of the simplex table corresponding to the basic
2009 variable x[k] */
2010 rlen = glp_eval_tab_row(P, k, rind, rval);
2011 xassert(0 <= rlen && rlen <= n);
2012 /* perform analysis */
2013 for (kase = -1; kase <= +1; kase += 2)
2014 { /* kase < 0 means objective coefficient c[k] is decreasing;
2015 kase > 0 means objective coefficient c[k] is increasing */
2016 /* note that decreasing c[k] is equivalent to increasing dual
2017 variable lambda[k] and vice versa; we need to correctly set
2018 the dir flag as required by the routine glp_dual_rtest */
2019 if (P->dir == GLP_MIN)
2020 dir = - kase;
2021 else if (P->dir == GLP_MAX)
2022 dir = + kase;
2023 else
2024 xassert(P != P);
2025 /* use the dual ratio test to determine non-basic variable
2026 x[q] whose reduced cost d[q] reaches zero bound first */
2027 rpiv = glp_dual_rtest(P, rlen, rind, rval, dir, 1e-9);
2028 if (rpiv == 0)
2029 { /* nothing limits changing c[k] */
2030 lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX);
2031 q = 0;
2032 /* x[k] keeps its current value */
2033 new_x = x;
2034 goto store;
2035 }
2036 /* non-basic variable x[q] limits changing coefficient c[k];
2037 determine its status and reduced cost d[k] in the current
2038 basis */
2039 xassert(1 <= rpiv && rpiv <= rlen);
2040 q = rind[rpiv];
2041 xassert(1 <= q && q <= m+n);
2042 if (q <= m)
2043 { row = P->row[q];
2044 stat = row->stat;
2045 d = row->dual;
2046 }
2047 else
2048 { col = P->col[q-m];
2049 stat = col->stat;
2050 d = col->dual;
2051 }
2052 /* note that delta d[q] = new d[q] - d[q] = - d[q], because
2053 new d[q] = 0; delta d[q] = alfa[k,q] * delta c[k], so
2054 delta c[k] = delta d[q] / alfa[k,q] = - d[q] / alfa[k,q] */
2055 xassert(rval[rpiv] != 0.0);
2056 delta = - d / rval[rpiv];
2057 /* compute new c[k] = c[k] + delta c[k], which is the limiting
2058 value of the objective coefficient c[k] */
2059 lim_coef = coef + delta;
2060 /* let c[k] continue decreasing/increasing that makes d[q]
2061 dual infeasible and forces x[q] to enter the basis;
2062 to perform the primal ratio test we need to know in which
2063 direction x[q] changes on entering the basis; we determine
2064 that analyzing the sign of delta d[q] (see above), since
2065 d[q] may be close to zero having wrong sign */
2066 /* let, for simplicity, the problem is minimization */
2067 if (kase < 0 && rval[rpiv] > 0.0 ||
2068 kase > 0 && rval[rpiv] < 0.0)
2069 { /* delta d[q] < 0, so d[q] being non-negative will become
2070 negative, so x[q] will increase */
2071 dir = +1;
2072 }
2073 else
2074 { /* delta d[q] > 0, so d[q] being non-positive will become
2075 positive, so x[q] will decrease */
2076 dir = -1;
2077 }
2078 /* if the problem is maximization, correct the direction */
2079 if (P->dir == GLP_MAX) dir = - dir;
2080 /* check that we didn't make a silly mistake */
2081 if (dir > 0)
2082 xassert(stat == GLP_NL || stat == GLP_NF);
2083 else
2084 xassert(stat == GLP_NU || stat == GLP_NF);
2085 /* compute column of the simplex table corresponding to the
2086 non-basic variable x[q] */
2087 clen = glp_eval_tab_col(P, q, cind, cval);
2088 /* make x[k] temporarily free (unbounded) */
2089 if (k <= m)
2090 { row = P->row[k];
2091 row->type = GLP_FR;
2092 row->lb = row->ub = 0.0;
2093 }
2094 else
2095 { col = P->col[k-m];
2096 col->type = GLP_FR;
2097 col->lb = col->ub = 0.0;
2098 }
2099 /* use the primal ratio test to determine some basic variable
2100 which leaves the basis */
2101 cpiv = glp_prim_rtest(P, clen, cind, cval, dir, 1e-9);
2102 /* restore original bounds of the basic variable x[k] */
2103 if (k <= m)
2104 { row = P->row[k];
2105 row->type = type;
2106 row->lb = lb, row->ub = ub;
2107 }
2108 else
2109 { col = P->col[k-m];
2110 col->type = type;
2111 col->lb = lb, col->ub = ub;
2112 }
2113 if (cpiv == 0)
2114 { /* non-basic variable x[q] can change unlimitedly */
2115 if (dir < 0 && rval[rpiv] > 0.0 ||
2116 dir > 0 && rval[rpiv] < 0.0)
2117 { /* delta x[k] = alfa[k,q] * delta x[q] < 0 */
2118 new_x = -DBL_MAX;
2119 }
2120 else
2121 { /* delta x[k] = alfa[k,q] * delta x[q] > 0 */
2122 new_x = +DBL_MAX;
2123 }
2124 goto store;
2125 }
2126 /* some basic variable x[p] limits changing non-basic variable
2127 x[q] in the adjacent basis */
2128 xassert(1 <= cpiv && cpiv <= clen);
2129 p = cind[cpiv];
2130 xassert(1 <= p && p <= m+n);
2131 xassert(p != k);
2132 if (p <= m)
2133 { row = P->row[p];
2134 xassert(row->stat == GLP_BS);
2135 ll = glp_get_row_lb(P, row->i);
2136 uu = glp_get_row_ub(P, row->i);
2137 xx = row->prim;
2138 }
2139 else
2140 { col = P->col[p-m];
2141 xassert(col->stat == GLP_BS);
2142 ll = glp_get_col_lb(P, col->j);
2143 uu = glp_get_col_ub(P, col->j);
2144 xx = col->prim;
2145 }
2146 /* determine delta x[p] = new x[p] - x[p] */
2147 if (dir < 0 && cval[cpiv] > 0.0 ||
2148 dir > 0 && cval[cpiv] < 0.0)
2149 { /* delta x[p] < 0, so x[p] goes toward its lower bound */
2150 xassert(ll != -DBL_MAX);
2151 delta = ll - xx;
2152 }
2153 else
2154 { /* delta x[p] > 0, so x[p] goes toward its upper bound */
2155 xassert(uu != +DBL_MAX);
2156 delta = uu - xx;
2157 }
2158 /* compute new x[k] = x[k] + alfa[k,q] * delta x[q], where
2159 delta x[q] = delta x[p] / alfa[p,q] */
2160 xassert(cval[cpiv] != 0.0);
2161 new_x = x + (rval[rpiv] / cval[cpiv]) * delta;
2162 store: /* store analysis results */
2163 if (kase < 0)
2164 { if (coef1 != NULL) *coef1 = lim_coef;
2165 if (var1 != NULL) *var1 = q;
2166 if (value1 != NULL) *value1 = new_x;
2167 }
2168 else
2169 { if (coef2 != NULL) *coef2 = lim_coef;
2170 if (var2 != NULL) *var2 = q;
2171 if (value2 != NULL) *value2 = new_x;
2172 }
2173 }
2174 /* free working arrays */
2175 xfree(cind);
2176 xfree(cval);
2177 xfree(rind);
2178 xfree(rval);
2179 return;
2180 }
2181
2182 /* eof */
2183