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