1
/* glpapi12.c (basis factorization and simplex tableau routines) */
3
/***********************************************************************
4
* This code is part of GLPK (GNU Linear Programming Kit).
6
* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7
* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
8
* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9
* E-mail: <mao@gnu.org>.
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.
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.
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
***********************************************************************/
27
/***********************************************************************
30
* glp_bf_exists - check if the basis factorization exists
34
* int glp_bf_exists(glp_prob *lp);
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)
45
ret = (lp->m == 0 || lp->valid);
49
/***********************************************************************
52
* glp_factorize - compute the basis factorization
56
* int glp_factorize(glp_prob *lp);
60
* The routine glp_factorize computes the basis factorization for the
61
* current basis associated with the specified problem object.
65
* 0 The basis factorization has been successfully computed.
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
73
* The basis matrix is singular within the working precision.
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;
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] */
87
/* build j-th column of the basic matrix, which is k-th column of
88
the scaled augmented matrix (I | -R*A*S) */
90
{ /* x[k] is auxiliary variable */
96
{ /* x[k] is structural variable */
98
for (aij = lp->col[k-m]->ptr; aij != NULL; aij = aij->c_next)
100
ind[len] = aij->row->i;
101
val[len] = - aij->row->rii * aij->val * aij->col->sjj;
107
static void copy_bfcp(glp_prob *lp);
109
int glp_factorize(glp_prob *lp)
112
GLPROW **row = lp->row;
113
GLPCOL **col = lp->col;
114
int *head = lp->head;
116
/* invalidate the basis factorization */
118
/* build the basis header */
120
for (k = 1; k <= m+n; k++)
122
{ stat = row[k]->stat;
126
{ stat = col[k-m]->stat;
132
{ /* too many basic variables */
144
{ /* too few basic variables */
148
/* try to factorize the basis matrix */
150
{ if (lp->bfd == NULL)
151
{ lp->bfd = bfd_create_it();
154
switch (bfd_factorize(lp->bfd, m, lp->head, b_col, lp))
159
/* singular matrix */
163
/* ill-conditioned matrix */
171
/* factorization successful */
173
fini: /* bring the return code to the calling program */
177
/***********************************************************************
180
* glp_bf_updated - check if the basis factorization has been updated
184
* int glp_bf_updated(glp_prob *lp);
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)
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);
199
cnt = (lp->m == 0 ? 0 : bfd_get_count(lp->bfd));
204
/***********************************************************************
207
* glp_get_bfcp - retrieve basis factorization control parameters
211
* void glp_get_bfcp(glp_prob *lp, glp_bfcp *parm);
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.
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;
225
{ parm->type = GLP_BF_FT;
227
parm->piv_tol = 0.10;
230
parm->eps_tol = 1e-15;
231
parm->max_gro = 1e+10;
233
parm->upd_tol = 1e-6;
238
memcpy(parm, bfcp, sizeof(glp_bfcp));
242
/***********************************************************************
245
* glp_set_bfcp - change basis factorization control parameters
249
* void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm);
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.
257
* New values of the control parameters should be passed in a structure
258
* glp_bfcp, which the parameter parm points to.
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;
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;
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);
291
void glp_set_bfcp(glp_prob *lp, const glp_bfcp *parm)
292
{ glp_bfcp *bfcp = lp->bfcp;
294
{ /* reset to default values */
296
xfree(bfcp), lp->bfcp = NULL;
299
{ /* set to specified values */
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",
307
if (bfcp->lu_size < 0)
308
xerror("glp_set_bfcp: lu_size = %d; invalid parameter\n",
310
if (!(0.0 < bfcp->piv_tol && bfcp->piv_tol < 1.0))
311
xerror("glp_set_bfcp: piv_tol = %g; invalid parameter\n",
313
if (bfcp->piv_lim < 1)
314
xerror("glp_set_bfcp: piv_lim = %d; invalid parameter\n",
316
if (!(bfcp->suhl == GLP_ON || bfcp->suhl == GLP_OFF))
317
xerror("glp_set_bfcp: suhl = %d; invalid parameter\n",
319
if (!(0.0 <= bfcp->eps_tol && bfcp->eps_tol <= 1e-6))
320
xerror("glp_set_bfcp: eps_tol = %g; invalid parameter\n",
322
if (bfcp->max_gro < 1.0)
323
xerror("glp_set_bfcp: max_gro = %g; invalid parameter\n",
325
if (!(1 <= bfcp->nfs_max && bfcp->nfs_max <= 32767))
326
xerror("glp_set_bfcp: nfs_max = %d; invalid parameter\n",
328
if (!(0.0 < bfcp->upd_tol && bfcp->upd_tol < 1.0))
329
xerror("glp_set_bfcp: upd_tol = %g; invalid parameter\n",
331
if (!(1 <= bfcp->nrs_max && bfcp->nrs_max <= 32767))
332
xerror("glp_set_bfcp: nrs_max = %d; invalid parameter\n",
334
if (bfcp->rs_size < 0)
335
xerror("glp_set_bfcp: rs_size = %d; invalid parameter\n",
337
if (bfcp->rs_size == 0)
338
bfcp->rs_size = 20 * bfcp->nrs_max;
340
if (lp->bfd != NULL) copy_bfcp(lp);
344
/***********************************************************************
347
* glp_get_bhead - retrieve the basis header information
351
* int glp_get_bhead(glp_prob *lp, int k);
355
* The routine glp_get_bhead returns the basis header information for
356
* the current basis associated with the specified problem object.
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);
373
/***********************************************************************
376
* glp_get_row_bind - retrieve row index in the basis header
380
* int glp_get_row_bind(glp_prob *lp, int i);
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"
394
if (!(1 <= i && i <= lp->m))
395
xerror("glp_get_row_bind: i = %d; row number out of range\n",
397
return lp->row[i]->bind;
400
/***********************************************************************
403
* glp_get_col_bind - retrieve column index in the basis header
407
* int glp_get_col_bind(glp_prob *lp, int j);
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"
421
if (!(1 <= j && j <= lp->n))
422
xerror("glp_get_col_bind: j = %d; column number out of range\n"
424
return lp->col[j]->bind;
427
/***********************************************************************
430
* glp_ftran - perform forward transformation (solve system B*x = b)
434
* void glp_ftran(glp_prob *lp, double x[]);
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.
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.
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
453
* A~" = (I | A") = (I | R*A*S) = (R*I*inv(R) | R*A*S) =
455
* = R*(I | A)*S~ = R*A~*S~,
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
460
* S~ = diag(inv(R) | S) (2)
462
* is an augmented diagonal scaling matrix.
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:
472
* A~" = (B" | N") = (R*B*SB | R*N*SN), (4)
474
* where SB and SN are parts of the augmented scaling matrix S~, which
475
* correspond to basic and non-basic variables, respectively. Therefore
479
* which is the scaled basis matrix. */
481
void glp_ftran(glp_prob *lp, double x[])
483
GLPROW **row = lp->row;
484
GLPCOL **col = lp->col;
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");
491
for (i = 1; i <= m; i++)
493
/* x" := inv(B")*b" */
494
if (m > 0) bfd_ftran(lp->bfd, x);
496
for (i = 1; i <= m; i++)
501
x[i] *= col[k-m]->sjj;
506
/***********************************************************************
509
* glp_btran - perform backward transformation (solve system B'*x = b)
513
* void glp_btran(glp_prob *lp, double x[]);
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.
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.
529
* See comments to the routine glp_ftran. */
531
void glp_btran(glp_prob *lp, double x[])
533
GLPROW **row = lp->row;
534
GLPCOL **col = lp->col;
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");
541
for (i = 1; i <= m; i++)
546
x[i] *= col[k-m]->sjj;
548
/* x" := inv[(B")']*b" */
549
if (m > 0) bfd_btran(lp->bfd, x);
551
for (i = 1; i <= m; i++)
556
/***********************************************************************
559
* glp_warm_up - "warm up" LP basis
563
* int glp_warm_up(glp_prob *P);
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).
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.
577
* 0 The operation has been successfully performed.
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
585
* The basis matrix is singular within the working precision.
588
* The basis matrix is ill-conditioned. */
590
int glp_warm_up(glp_prob *P)
595
double eps, temp, *work;
596
/* invalidate basic solution */
597
P->pbs_stat = P->dbs_stat = GLP_UNDEF;
600
for (i = 1; i <= P->m; i++)
602
row->prim = row->dual = 0.0;
604
for (j = 1; j <= P->n; j++)
606
col->prim = col->dual = 0.0;
608
/* compute the basis factorization, if necessary */
609
if (!glp_bf_exists(P))
610
{ ret = glp_factorize(P);
611
if (ret != 0) goto done;
613
/* allocate working array */
614
work = xcalloc(1+P->m, sizeof(double));
615
/* determine and store values of non-basic variables, compute
617
for (i = 1; i <= P->m; i++)
619
for (i = 1; i <= P->m; i++)
621
if (row->stat == GLP_BS)
623
else if (row->stat == GLP_NL)
625
else if (row->stat == GLP_NU)
627
else if (row->stat == GLP_NF)
629
else if (row->stat == GLP_NS)
633
/* N[j] is i-th column of matrix (I|-A) */
634
work[i] -= row->prim;
636
for (j = 1; j <= P->n; j++)
638
if (col->stat == GLP_BS)
640
else if (col->stat == GLP_NL)
642
else if (col->stat == GLP_NU)
644
else if (col->stat == GLP_NF)
646
else if (col->stat == GLP_NS)
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;
656
/* compute vector of basic variables xB = - inv(B) * N * xN */
658
/* store values of basic variables, check primal feasibility */
659
P->pbs_stat = GLP_FEAS;
660
for (i = 1; i <= P->m; i++)
662
if (row->stat != GLP_BS)
664
row->prim = work[row->bind];
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;
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;
677
for (j = 1; j <= P->n; j++)
679
if (col->stat != GLP_BS)
681
col->prim = work[col->bind];
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;
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;
694
/* compute value of the objective function */
696
for (j = 1; j <= P->n; j++)
698
P->obj_val += col->coef * col->prim;
700
/* build vector cB of objective coefficients at basic variables */
701
for (i = 1; i <= P->m; i++)
703
for (j = 1; j <= P->n; j++)
705
if (col->stat == GLP_BS)
706
work[col->bind] = col->coef;
708
/* compute vector of simplex multipliers pi = inv(B') * cB */
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++)
715
if (row->stat == GLP_BS)
719
/* N[j] is i-th column of matrix (I|-A) */
720
row->dual = - work[i];
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;
727
for (j = 1; j <= P->n; j++)
729
if (col->stat == GLP_BS)
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];
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;
743
/* free working array */
749
/***********************************************************************
752
* glp_eval_tab_row - compute row of the simplex tableau
756
* int glp_eval_tab_row(glp_prob *lp, int k, int ind[], double val[]);
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.
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.
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).
776
* The computed row shows how the specified basic variable x[k] = xB[i]
777
* depends on non-basic variables:
779
* xB[i] = alfa[i,1]*xN[1] + alfa[i,2]*xN[2] + ... + alfa[i,n]*xN[n],
781
* where alfa[i,j] are elements of the simplex table row, xN[j] are
782
* non-basic (auxiliary and structural) variables.
786
* The routine returns number of non-zero elements in the simplex table
787
* row stored in the arrays ind and val.
791
* The system of equality constraints of the LP problem is:
795
* where xR is the vector of auxliary variables, xS is the vector of
796
* structural variables, A is the matrix of constraint coefficients.
798
* The system (1) can be written in homogenous form as follows:
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.
806
* By definition for the current basis we have:
810
* where B is the basis matrix. Thus, the system (2) can be written as:
812
* B * xB + N * xN = 0. (4)
814
* From (4) it follows that:
820
* A^ = - inv(B) * N (6)
822
* is called the simplex table.
824
* It is understood that i-th row of the simplex table is:
826
* e * A^ = - e * inv(B) * N, (7)
828
* where e is a unity vector with e[i] = 1.
830
* To compute i-th row of the simplex table the routine first computes
831
* i-th row of the inverse:
833
* rho = inv(B') * e, (8)
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:
838
* alfa[i,j] = - rho * N[j] for all j, (9)
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[])
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"
851
if (!(1 <= k && k <= m+n))
852
xerror("glp_eval_tab_row: k = %d; variable number out of range"
854
/* determine xB[i] which corresponds to x[k] */
856
i = glp_get_row_bind(lp, k);
858
i = glp_get_col_bind(lp, k-m);
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;
870
/* compute i-th row of the simplex table */
872
for (k = 1; k <= m+n; k++)
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) */
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);
886
for (t = 1; t <= lll; t++) alfa += rho[iii[t]] * vvv[t];
888
/* store alfa[i,j] */
889
if (alfa != 0.0) len++, ind[len] = k, val[len] = alfa;
892
/* free working arrays */
896
/* return to the calling program */
900
/***********************************************************************
903
* glp_eval_tab_col - compute column of the simplex tableau
907
* int glp_eval_tab_col(glp_prob *lp, int k, int ind[], double val[]);
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.
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.
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).
927
* The computed column shows how basic variables depend on the specified
928
* non-basic variable x[k] = xN[j]:
930
* xB[1] = ... + alfa[1,j]*xN[j] + ...
931
* xB[2] = ... + alfa[2,j]*xN[j] + ...
933
* xB[m] = ... + alfa[m,j]*xN[j] + ...
935
* where alfa[i,j] are elements of the simplex table column, xB[i] are
936
* basic (auxiliary and structural) variables.
940
* The routine returns number of non-zero elements in the simplex table
941
* column stored in the arrays ind and val.
945
* As it was explained in comments to the routine glp_eval_tab_row (see
946
* above) the simplex table is the following matrix:
948
* A^ = - inv(B) * N. (1)
950
* Therefore j-th column of the simplex table is:
952
* A^ * e = - inv(B) * N * e = - inv(B) * N[j], (2)
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[])
963
if (!(m == 0 || lp->valid))
964
xerror("glp_eval_tab_col: basis factorization does not exist\n"
966
if (!(1 <= k && k <= m+n))
967
xerror("glp_eval_tab_col: k = %d; variable number out of range"
970
stat = glp_get_row_stat(lp, k);
972
stat = glp_get_col_stat(lp, k-m);
974
xerror("glp_eval_tab_col: k = %d; variable must be non-basic",
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;
980
{ /* x[k] is auxiliary variable, so N[k] is a unity column */
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];
989
/* compute column of the simplex table, which corresponds to the
990
specified non-basic variable x[k] */
993
for (t = 1; t <= m; t++)
996
ind[len] = glp_get_bhead(lp, t);
1001
/* return to the calling program */
1005
/***********************************************************************
1008
* glp_transform_row - transform explicitly specified row
1012
* int glp_transform_row(glp_prob *P, int len, int ind[], double val[]);
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.
1020
* The explicitly specified row may be thought as a linear form:
1022
* x = a[1]*x[m+1] + a[2]*x[m+2] + ... + a[n]*x[m+n], (1)
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.
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.
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:
1037
* x = alfa[1]*xN[1] + alfa[2]*xN[2] + ... + alfa[n]*xN[n], (2)
1039
* where xN[j] are non-basic (auxiliary or structural) variables, n is
1040
* the number of columns in the LP problem object.
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.
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.
1057
* The explicitly specified row (1) is transformed in the same way as it
1058
* were the objective function row.
1060
* From (1) it follows that:
1062
* x = aB * xB + aN * xN, (3)
1064
* where xB is the vector of basic variables, xN is the vector of
1065
* non-basic variables.
1067
* The simplex table, which corresponds to the current basis, is:
1069
* xB = [-inv(B) * N] * xN. (4)
1071
* Therefore substituting xB from (4) to (3) we have:
1073
* x = aB * [-inv(B) * N] * xN + aN * xN =
1075
* = rho * (-N) * xN + aN * xN = alfa * xN,
1079
* rho = inv(B') * aB, (6)
1083
* alfa = aN + rho * (-N) (7)
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 "
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",
1101
for (t = 1; t <= len; t++)
1103
if (!(1 <= j && j <= n))
1104
xerror("glp_transform_row: ind[%d] = %d; column index out o"
1107
xerror("glp_transform_row: val[%d] = 0; zero coefficient no"
1110
xerror("glp_transform_row: ind[%d] = %d; duplicate column i"
1111
"ndices not allowed\n", t, j);
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]);
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 */
1126
for (i = 1; i <= m; i++)
1127
{ if (glp_get_row_stat(P, i) != GLP_BS)
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)
1142
lll = glp_get_mat_col(P, j, iii, vvv);
1143
for (t = 1; t <= lll; t++) alfa += vvv[t] * rho[iii[t]];
1159
/***********************************************************************
1162
* glp_transform_col - transform explicitly specified column
1166
* int glp_transform_col(glp_prob *P, int len, int ind[], double val[]);
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.
1174
* The explicitly specified column may be thought as if it were added
1175
* to the original system of equality constraints:
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
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.
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.
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:
1196
* xB[1] = ... + alfa[1]*x
1197
* xB[2] = ... + alfa[2]*x (2)
1199
* xB[m] = ... + alfa[m]*x
1201
* where xB are basic (auxiliary and structural) variables, m is the
1202
* number of rows in the problem object.
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
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.
1220
* The explicitly specified column (1) is transformed in the same way
1221
* as any other column of the constraint matrix using the formula:
1223
* alfa = inv(B) * a, (3)
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[])
1230
if (!glp_bf_exists(P))
1231
xerror("glp_transform_col: basis factorization does not exist "
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",
1240
for (t = 1; t <= len; t++)
1242
if (!(1 <= i && i <= m))
1243
xerror("glp_transform_col: ind[%d] = %d; row index out of r"
1246
xerror("glp_transform_col: val[%d] = 0; zero coefficient no"
1249
xerror("glp_transform_col: ind[%d] = %d; duplicate row indi"
1250
"ces not allowed\n", t, i);
1253
/* solve the system B*a = alfa to compute the vector alfa */
1254
alfa = a, glp_ftran(P, alfa);
1255
/* store resultant coefficients */
1257
for (i = 1; i <= m; i++)
1258
{ if (alfa[i] != 0.0)
1260
ind[len] = glp_get_bhead(P, i);
1268
/***********************************************************************
1271
* glp_prim_rtest - perform primal ratio test
1275
* int glp_prim_rtest(glp_prob *P, int len, const int ind[],
1276
* const double val[], int dir, double eps);
1280
* The routine glp_prim_rtest performs the primal ratio test using an
1281
* explicitly specified column of the simplex table.
1283
* The current basic solution associated with the LP problem object
1284
* must be primal feasible.
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):
1290
* xB[1] = ... + alfa[1] * x + ...
1291
* xB[2] = ... + alfa[2] * x + ... (*)
1293
* xB[m] = ... + alfa[m] * x + ...
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.
1303
* The parameter dir specifies direction in which the variable x changes
1304
* on entering the basis: +1 means increasing, -1 means decreasing.
1306
* The parameter eps is an absolute tolerance (small positive number)
1307
* used by the routine to skip small alfa[j] of the row (*).
1309
* The routine determines which basic variable (among specified in
1310
* ind[1], ..., ind[len]) should leave the basis in order to keep primal
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.
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 "
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 */
1345
if (!(1 <= k && k <= m+n))
1346
xerror("glp_prim_rtest: ind[%d] = %d; variable number out o"
1348
/* determine type, bounds, status and primal value of basic
1349
variable xB[i] = x[k] in the current basic solution */
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);
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);
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
1371
alfa = (dir > 0 ? + val[t] : - val[t]);
1372
/* analyze main cases */
1374
{ /* xB[i] is free variable */
1377
else if (type == GLP_LO)
1378
lo: { /* xB[i] has an lower bound */
1379
if (alfa > - eps) continue;
1380
temp = (lb - beta) / alfa;
1382
else if (type == GLP_UP)
1383
up: { /* xB[i] has an upper bound */
1384
if (alfa < + eps) continue;
1385
temp = (ub - beta) / alfa;
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;
1391
else if (type == GLP_FX)
1392
{ /* xB[i] is fixed variable */
1393
if (- eps < alfa && alfa < + eps) continue;
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);
1408
/* return index of the pivot element chosen */
1412
/***********************************************************************
1415
* glp_dual_rtest - perform dual ratio test
1419
* int glp_dual_rtest(glp_prob *P, int len, const int ind[],
1420
* const double val[], int dir, double eps);
1424
* The routine glp_dual_rtest performs the dual ratio test using an
1425
* explicitly specified row of the simplex table.
1427
* The current basic solution associated with the LP problem object
1428
* must be dual feasible.
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:
1434
* x = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n]. (*)
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.
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.
1448
* The parameter eps is an absolute tolerance (small positive number)
1449
* used by the routine to skip small alfa[j] of the row (*).
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
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.
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")
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 */
1489
if (!(1 <= k && k <= m+n))
1490
xerror("glp_dual_rtest: ind[%d] = %d; variable number out o"
1492
/* determine status and reduced cost of non-basic variable
1493
x[k] = xN[j] in the current basic solution */
1495
{ stat = glp_get_row_stat(P, k);
1496
cost = glp_get_row_dual(P, k);
1499
{ stat = glp_get_col_stat(P, k-m);
1500
cost = glp_get_col_dual(P, k-m);
1503
xerror("glp_dual_rtest: ind[%d] = %d; basic variable not al"
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
1509
alfa = (dir > 0 ? + val[t] : - val[t]);
1510
/* analyze main cases */
1512
{ /* xN[j] is on its lower bound */
1513
if (alfa < + eps) continue;
1514
temp = (obj * cost) / alfa;
1516
else if (stat == GLP_NU)
1517
{ /* xN[j] is on its upper bound */
1518
if (alfa > - eps) continue;
1519
temp = (obj * cost) / alfa;
1521
else if (stat == GLP_NF)
1522
{ /* xN[j] is non-basic free variable */
1523
if (- eps < alfa && alfa < + eps) continue;
1526
else if (stat == GLP_NS)
1527
{ /* xN[j] is non-basic fixed variable */
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);
1542
/* return index of the pivot element chosen */
1546
/***********************************************************************
1549
* glp_analyze_row - simulate one iteration of dual simplex method
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);
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).
1566
* The current basic solution associated with the problem object passed
1567
* to the routine must be dual feasible, and its primal components must
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:
1577
* y = alfa[1] * xN[1] + alfa[2] * xN[2] + ... + alfa[n] * xN[n],
1579
* where y is an auxiliary variable of the row, alfa[j] is an influence
1580
* coefficient, xN[j] is a non-basic variable.
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.
1589
* The parameters type and rhs specify the row type and its right-hand
1592
* type = GLP_LO: y = sum alfa[j] * xN[j] >= rhs
1594
* type = GLP_UP: y = sum alfa[j] * xN[j] <= rhs
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.
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):
1604
* piv index in the array ind and val, 1 <= piv <= len, determining
1605
* the non-basic variable, which would enter the adjacent basis;
1607
* x value of the non-basic variable in the current basis;
1609
* dx difference between values of the non-basic variable in the
1610
* adjacent and current bases, dx = x.new - x.old;
1612
* y value of the row (i.e. of its auxiliary variable) in the
1615
* dy difference between values of the row in the adjacent and
1616
* current bases, dy = y.new - y.old;
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 "
1632
if (P->dbs_stat != GLP_FEAS)
1633
xerror("glp_analyze_row: basic solution is not dual feasible\n"
1635
/* compute the row value y = sum alfa[j] * xN[j] in the current
1637
if (!(0 <= len && len <= P->n))
1638
xerror("glp_analyze_row: len = %d; invalid row length\n", len);
1640
for (t = 1; t <= len; t++)
1641
{ /* determine value of x[k] = xN[j] in the current basis */
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);
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;
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;
1662
/* check if the row is primal infeasible in the current basis,
1663
i.e. the constraint is violated at the current point */
1666
{ /* the constraint is not violated */
1670
/* in the adjacent basis y goes to its lower bound */
1673
else if (type == GLP_UP)
1675
{ /* the constraint is not violated */
1679
/* in the adjacent basis y goes to its upper bound */
1683
xerror("glp_analyze_row: type = %d; invalid parameter\n",
1685
/* compute dy = y.new - y.old */
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);
1691
{ /* no dual feasible adjacent basis exists */
1695
/* non-basic variable x[k] = xN[j] should enter the basis */
1697
xassert(1 <= k && k <= P->m+P->n);
1698
/* determine its value in the current basis */
1700
x = P->row[k]->prim;
1702
x = P->col[k-P->m]->prim;
1703
/* compute dx = x.new - x.old = dy / alfa[j] */
1704
xassert(val[piv] != 0.0);
1706
/* compute dz = z.new - z.old = d[j] * dx, where d[j] is reduced
1707
cost of xN[j] in the current basis */
1709
dz = P->row[k]->dual * dx;
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;
1724
{ /* example program for the routine glp_analyze_row */
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);
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;
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
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
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;
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;
1786
glp_simplex(P, &parm);
1787
/* the current objective value is 301.861 */
1792
/***********************************************************************
1795
* glp_analyze_bound - analyze active bound of non-basic variable
1799
* void glp_analyze_bound(glp_prob *P, int k, double *limit1, int *var1,
1800
* double *limit2, int *var2);
1804
* The routine glp_analyze_bound analyzes the effect of varying the
1805
* active bound of specified non-basic variable.
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).
1811
* Note that the current basic solution must be optimal, and the basis
1812
* factorization must exist.
1814
* Results of the analysis have the following meaning.
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.
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.
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.
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)
1838
int m, n, stat, kase, p, len, piv, *ind;
1839
double x, new_x, ll, uu, xx, delta, *val;
1841
if (P == NULL || P->magic != GLP_PROB_MAGIC)
1842
xerror("glp_analyze_bound: P = %p; invalid problem object\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"
1852
/* retrieve information about the specified non-basic variable
1853
x[k] whose active bound is to be analyzed */
1860
{ col = P->col[k-m];
1865
xerror("glp_analyze_bound: k = %d; basic variable not allowed "
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);
1882
{ /* nothing limits changing the active bound of x[k] */
1884
new_x = (kase < 0 ? -DBL_MAX : +DBL_MAX);
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);
1893
ll = glp_get_row_lb(P, row->i);
1894
uu = glp_get_row_ub(P, row->i);
1899
{ col = P->col[p-m];
1900
ll = glp_get_col_lb(P, col->j);
1901
uu = glp_get_col_ub(P, col->j);
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);
1914
{ /* delta x[p] > 0, so x[p] goes toward its upper bound */
1915
xassert(uu != +DBL_MAX);
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 */
1925
{ if (value1 != NULL) *value1 = new_x;
1926
if (var1 != NULL) *var1 = p;
1929
{ if (value2 != NULL) *value2 = new_x;
1930
if (var2 != NULL) *var2 = p;
1933
/* free working arrays */
1939
/***********************************************************************
1942
* glp_analyze_coef - analyze objective coefficient at basic variable
1946
* void glp_analyze_coef(glp_prob *P, int k, double *coef1, int *var1,
1947
* double *value1, double *coef2, int *var2, double *value2);
1951
* The routine glp_analyze_coef analyzes the effect of varying the
1952
* objective coefficient at specified basic variable.
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).
1958
* Note that the current basic solution must be optimal, and the basis
1959
* factorization must exist.
1961
* Results of the analysis have the following meaning.
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.
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.
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.
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.
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.
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,
2003
double lb, ub, coef, x, lim_coef, new_x, d, delta, ll, uu, xx,
2006
if (P == NULL || P->magic != GLP_PROB_MAGIC)
2007
xerror("glp_analyze_coef: P = %p; invalid problem object\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"
2017
/* retrieve information about the specified basic variable x[k]
2018
whose objective coefficient c[k] is to be analyzed */
2029
{ col = P->col[k-m];
2038
xerror("glp_analyze_coef: k = %d; non-basic variable not allow"
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
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)
2058
else if (P->dir == GLP_MAX)
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);
2066
{ /* nothing limits changing c[k] */
2067
lim_coef = (kase < 0 ? -DBL_MAX : +DBL_MAX);
2069
/* x[k] keeps its current value */
2073
/* non-basic variable x[q] limits changing coefficient c[k];
2074
determine its status and reduced cost d[k] in the current
2076
xassert(1 <= rpiv && rpiv <= rlen);
2078
xassert(1 <= q && q <= m+n);
2085
{ col = P->col[q-m];
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 */
2111
{ /* delta d[q] > 0, so d[q] being non-positive will become
2112
positive, so x[q] will decrease */
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 */
2119
xassert(stat == GLP_NL || stat == GLP_NF);
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) */
2129
row->lb = row->ub = 0.0;
2132
{ col = P->col[k-m];
2134
col->lb = col->ub = 0.0;
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] */
2143
row->lb = lb, row->ub = ub;
2146
{ col = P->col[k-m];
2148
col->lb = lb, col->ub = ub;
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 */
2158
{ /* delta x[k] = alfa[k,q] * delta x[q] > 0 */
2163
/* some basic variable x[p] limits changing non-basic variable
2164
x[q] in the adjacent basis */
2165
xassert(1 <= cpiv && cpiv <= clen);
2167
xassert(1 <= p && p <= m+n);
2171
xassert(row->stat == GLP_BS);
2172
ll = glp_get_row_lb(P, row->i);
2173
uu = glp_get_row_ub(P, row->i);
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);
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);
2191
{ /* delta x[p] > 0, so x[p] goes toward its upper bound */
2192
xassert(uu != +DBL_MAX);
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 */
2201
{ if (coef1 != NULL) *coef1 = lim_coef;
2202
if (var1 != NULL) *var1 = q;
2203
if (value1 != NULL) *value1 = new_x;
2206
{ if (coef2 != NULL) *coef2 = lim_coef;
2207
if (var2 != NULL) *var2 = q;
2208
if (value2 != NULL) *value2 = new_x;
2211
/* free working arrays */
added
removed
optional/glpk/glpapi12.c
