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