1 | /* glpbfd.c (LP basis factorization driver) */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
9 | * E-mail: <mao@gnu.org>. |
---|
10 | * |
---|
11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
12 | * under the terms of the GNU General Public License as published by |
---|
13 | * the Free Software Foundation, either version 3 of the License, or |
---|
14 | * (at your option) any later version. |
---|
15 | * |
---|
16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
19 | * License for more details. |
---|
20 | * |
---|
21 | * You should have received a copy of the GNU General Public License |
---|
22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
23 | ***********************************************************************/ |
---|
24 | |
---|
25 | typedef struct BFD BFD; |
---|
26 | |
---|
27 | #define GLPBFD_PRIVATE |
---|
28 | #include "glpapi.h" |
---|
29 | #include "glpfhv.h" |
---|
30 | #include "glplpf.h" |
---|
31 | |
---|
32 | /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ |
---|
33 | |
---|
34 | #define M_MAX 100000000 /* = 100*10^6 */ |
---|
35 | /* maximal order of the basis matrix */ |
---|
36 | |
---|
37 | struct BFD |
---|
38 | { /* LP basis factorization */ |
---|
39 | int valid; |
---|
40 | /* factorization is valid only if this flag is set */ |
---|
41 | int type; |
---|
42 | /* factorization type: |
---|
43 | GLP_BF_FT - LUF + Forrest-Tomlin |
---|
44 | GLP_BF_BG - LUF + Schur compl. + Bartels-Golub |
---|
45 | GLP_BF_GR - LUF + Schur compl. + Givens rotation */ |
---|
46 | FHV *fhv; |
---|
47 | /* LP basis factorization (GLP_BF_FT) */ |
---|
48 | LPF *lpf; |
---|
49 | /* LP basis factorization (GLP_BF_BG, GLP_BF_GR) */ |
---|
50 | int lu_size; /* luf.sv_size */ |
---|
51 | double piv_tol; /* luf.piv_tol */ |
---|
52 | int piv_lim; /* luf.piv_lim */ |
---|
53 | int suhl; /* luf.suhl */ |
---|
54 | double eps_tol; /* luf.eps_tol */ |
---|
55 | double max_gro; /* luf.max_gro */ |
---|
56 | int nfs_max; /* fhv.hh_max */ |
---|
57 | double upd_tol; /* fhv.upd_tol */ |
---|
58 | int nrs_max; /* lpf.n_max */ |
---|
59 | int rs_size; /* lpf.v_size */ |
---|
60 | /* internal control parameters */ |
---|
61 | int upd_lim; |
---|
62 | /* the factorization update limit */ |
---|
63 | int upd_cnt; |
---|
64 | /* the factorization update count */ |
---|
65 | }; |
---|
66 | |
---|
67 | /*********************************************************************** |
---|
68 | * NAME |
---|
69 | * |
---|
70 | * bfd_create_it - create LP basis factorization |
---|
71 | * |
---|
72 | * SYNOPSIS |
---|
73 | * |
---|
74 | * #include "glpbfd.h" |
---|
75 | * BFD *bfd_create_it(void); |
---|
76 | * |
---|
77 | * DESCRIPTION |
---|
78 | * |
---|
79 | * The routine bfd_create_it creates a program object, which represents |
---|
80 | * a factorization of LP basis. |
---|
81 | * |
---|
82 | * RETURNS |
---|
83 | * |
---|
84 | * The routine bfd_create_it returns a pointer to the object created. */ |
---|
85 | |
---|
86 | BFD *bfd_create_it(void) |
---|
87 | { BFD *bfd; |
---|
88 | bfd = xmalloc(sizeof(BFD)); |
---|
89 | bfd->valid = 0; |
---|
90 | bfd->type = GLP_BF_FT; |
---|
91 | bfd->fhv = NULL; |
---|
92 | bfd->lpf = NULL; |
---|
93 | bfd->lu_size = 0; |
---|
94 | bfd->piv_tol = 0.10; |
---|
95 | bfd->piv_lim = 4; |
---|
96 | bfd->suhl = 1; |
---|
97 | bfd->eps_tol = 1e-15; |
---|
98 | bfd->max_gro = 1e+10; |
---|
99 | bfd->nfs_max = 100; |
---|
100 | bfd->upd_tol = 1e-6; |
---|
101 | bfd->nrs_max = 100; |
---|
102 | bfd->rs_size = 1000; |
---|
103 | bfd->upd_lim = -1; |
---|
104 | bfd->upd_cnt = 0; |
---|
105 | return bfd; |
---|
106 | } |
---|
107 | |
---|
108 | /**********************************************************************/ |
---|
109 | |
---|
110 | void bfd_set_parm(BFD *bfd, const void *_parm) |
---|
111 | { /* change LP basis factorization control parameters */ |
---|
112 | const glp_bfcp *parm = _parm; |
---|
113 | xassert(bfd != NULL); |
---|
114 | bfd->type = parm->type; |
---|
115 | bfd->lu_size = parm->lu_size; |
---|
116 | bfd->piv_tol = parm->piv_tol; |
---|
117 | bfd->piv_lim = parm->piv_lim; |
---|
118 | bfd->suhl = parm->suhl; |
---|
119 | bfd->eps_tol = parm->eps_tol; |
---|
120 | bfd->max_gro = parm->max_gro; |
---|
121 | bfd->nfs_max = parm->nfs_max; |
---|
122 | bfd->upd_tol = parm->upd_tol; |
---|
123 | bfd->nrs_max = parm->nrs_max; |
---|
124 | bfd->rs_size = parm->rs_size; |
---|
125 | return; |
---|
126 | } |
---|
127 | |
---|
128 | /*********************************************************************** |
---|
129 | * NAME |
---|
130 | * |
---|
131 | * bfd_factorize - compute LP basis factorization |
---|
132 | * |
---|
133 | * SYNOPSIS |
---|
134 | * |
---|
135 | * #include "glpbfd.h" |
---|
136 | * int bfd_factorize(BFD *bfd, int m, int bh[], int (*col)(void *info, |
---|
137 | * int j, int ind[], double val[]), void *info); |
---|
138 | * |
---|
139 | * DESCRIPTION |
---|
140 | * |
---|
141 | * The routine bfd_factorize computes the factorization of the basis |
---|
142 | * matrix B specified by the routine col. |
---|
143 | * |
---|
144 | * The parameter bfd specified the basis factorization data structure |
---|
145 | * created with the routine bfd_create_it. |
---|
146 | * |
---|
147 | * The parameter m specifies the order of B, m > 0. |
---|
148 | * |
---|
149 | * The array bh specifies the basis header: bh[j], 1 <= j <= m, is the |
---|
150 | * number of j-th column of B in some original matrix. The array bh is |
---|
151 | * optional and can be specified as NULL. |
---|
152 | * |
---|
153 | * The formal routine col specifies the matrix B to be factorized. To |
---|
154 | * obtain j-th column of A the routine bfd_factorize calls the routine |
---|
155 | * col with the parameter j (1 <= j <= n). In response the routine col |
---|
156 | * should store row indices and numerical values of non-zero elements |
---|
157 | * of j-th column of B to locations ind[1,...,len] and val[1,...,len], |
---|
158 | * respectively, where len is the number of non-zeros in j-th column |
---|
159 | * returned on exit. Neither zero nor duplicate elements are allowed. |
---|
160 | * |
---|
161 | * The parameter info is a transit pointer passed to the routine col. |
---|
162 | * |
---|
163 | * RETURNS |
---|
164 | * |
---|
165 | * 0 The factorization has been successfully computed. |
---|
166 | * |
---|
167 | * BFD_ESING |
---|
168 | * The specified matrix is singular within the working precision. |
---|
169 | * |
---|
170 | * BFD_ECOND |
---|
171 | * The specified matrix is ill-conditioned. |
---|
172 | * |
---|
173 | * For more details see comments to the routine luf_factorize. */ |
---|
174 | |
---|
175 | int bfd_factorize(BFD *bfd, int m, const int bh[], int (*col) |
---|
176 | (void *info, int j, int ind[], double val[]), void *info) |
---|
177 | { LUF *luf; |
---|
178 | int nov, ret; |
---|
179 | xassert(bfd != NULL); |
---|
180 | xassert(1 <= m && m <= M_MAX); |
---|
181 | /* invalidate the factorization */ |
---|
182 | bfd->valid = 0; |
---|
183 | /* create the factorization, if necessary */ |
---|
184 | nov = 0; |
---|
185 | switch (bfd->type) |
---|
186 | { case GLP_BF_FT: |
---|
187 | if (bfd->lpf != NULL) |
---|
188 | lpf_delete_it(bfd->lpf), bfd->lpf = NULL; |
---|
189 | if (bfd->fhv == NULL) |
---|
190 | bfd->fhv = fhv_create_it(), nov = 1; |
---|
191 | break; |
---|
192 | case GLP_BF_BG: |
---|
193 | case GLP_BF_GR: |
---|
194 | if (bfd->fhv != NULL) |
---|
195 | fhv_delete_it(bfd->fhv), bfd->fhv = NULL; |
---|
196 | if (bfd->lpf == NULL) |
---|
197 | bfd->lpf = lpf_create_it(), nov = 1; |
---|
198 | break; |
---|
199 | default: |
---|
200 | xassert(bfd != bfd); |
---|
201 | } |
---|
202 | /* set control parameters specific to LUF */ |
---|
203 | if (bfd->fhv != NULL) |
---|
204 | luf = bfd->fhv->luf; |
---|
205 | else if (bfd->lpf != NULL) |
---|
206 | luf = bfd->lpf->luf; |
---|
207 | else |
---|
208 | xassert(bfd != bfd); |
---|
209 | if (nov) luf->new_sva = bfd->lu_size; |
---|
210 | luf->piv_tol = bfd->piv_tol; |
---|
211 | luf->piv_lim = bfd->piv_lim; |
---|
212 | luf->suhl = bfd->suhl; |
---|
213 | luf->eps_tol = bfd->eps_tol; |
---|
214 | luf->max_gro = bfd->max_gro; |
---|
215 | /* set control parameters specific to FHV */ |
---|
216 | if (bfd->fhv != NULL) |
---|
217 | { if (nov) bfd->fhv->hh_max = bfd->nfs_max; |
---|
218 | bfd->fhv->upd_tol = bfd->upd_tol; |
---|
219 | } |
---|
220 | /* set control parameters specific to LPF */ |
---|
221 | if (bfd->lpf != NULL) |
---|
222 | { if (nov) bfd->lpf->n_max = bfd->nrs_max; |
---|
223 | if (nov) bfd->lpf->v_size = bfd->rs_size; |
---|
224 | } |
---|
225 | /* try to factorize the basis matrix */ |
---|
226 | if (bfd->fhv != NULL) |
---|
227 | { switch (fhv_factorize(bfd->fhv, m, col, info)) |
---|
228 | { case 0: |
---|
229 | break; |
---|
230 | case FHV_ESING: |
---|
231 | ret = BFD_ESING; |
---|
232 | goto done; |
---|
233 | case FHV_ECOND: |
---|
234 | ret = BFD_ECOND; |
---|
235 | goto done; |
---|
236 | default: |
---|
237 | xassert(bfd != bfd); |
---|
238 | } |
---|
239 | } |
---|
240 | else if (bfd->lpf != NULL) |
---|
241 | { switch (lpf_factorize(bfd->lpf, m, bh, col, info)) |
---|
242 | { case 0: |
---|
243 | /* set the Schur complement update type */ |
---|
244 | switch (bfd->type) |
---|
245 | { case GLP_BF_BG: |
---|
246 | /* Bartels-Golub update */ |
---|
247 | bfd->lpf->scf->t_opt = SCF_TBG; |
---|
248 | break; |
---|
249 | case GLP_BF_GR: |
---|
250 | /* Givens rotation update */ |
---|
251 | bfd->lpf->scf->t_opt = SCF_TGR; |
---|
252 | break; |
---|
253 | default: |
---|
254 | xassert(bfd != bfd); |
---|
255 | } |
---|
256 | break; |
---|
257 | case LPF_ESING: |
---|
258 | ret = BFD_ESING; |
---|
259 | goto done; |
---|
260 | case LPF_ECOND: |
---|
261 | ret = BFD_ECOND; |
---|
262 | goto done; |
---|
263 | default: |
---|
264 | xassert(bfd != bfd); |
---|
265 | } |
---|
266 | } |
---|
267 | else |
---|
268 | xassert(bfd != bfd); |
---|
269 | /* the basis matrix has been successfully factorized */ |
---|
270 | bfd->valid = 1; |
---|
271 | bfd->upd_cnt = 0; |
---|
272 | ret = 0; |
---|
273 | done: /* return to the calling program */ |
---|
274 | return ret; |
---|
275 | } |
---|
276 | |
---|
277 | /*********************************************************************** |
---|
278 | * NAME |
---|
279 | * |
---|
280 | * bfd_ftran - perform forward transformation (solve system B*x = b) |
---|
281 | * |
---|
282 | * SYNOPSIS |
---|
283 | * |
---|
284 | * #include "glpbfd.h" |
---|
285 | * void bfd_ftran(BFD *bfd, double x[]); |
---|
286 | * |
---|
287 | * DESCRIPTION |
---|
288 | * |
---|
289 | * The routine bfd_ftran performs forward transformation, i.e. solves |
---|
290 | * the system B*x = b, where B is the basis matrix, x is the vector of |
---|
291 | * unknowns to be computed, b is the vector of right-hand sides. |
---|
292 | * |
---|
293 | * On entry elements of the vector b should be stored in dense format |
---|
294 | * in locations x[1], ..., x[m], where m is the number of rows. On exit |
---|
295 | * the routine stores elements of the vector x in the same locations. */ |
---|
296 | |
---|
297 | void bfd_ftran(BFD *bfd, double x[]) |
---|
298 | { xassert(bfd != NULL); |
---|
299 | xassert(bfd->valid); |
---|
300 | if (bfd->fhv != NULL) |
---|
301 | fhv_ftran(bfd->fhv, x); |
---|
302 | else if (bfd->lpf != NULL) |
---|
303 | lpf_ftran(bfd->lpf, x); |
---|
304 | else |
---|
305 | xassert(bfd != bfd); |
---|
306 | return; |
---|
307 | } |
---|
308 | |
---|
309 | /*********************************************************************** |
---|
310 | * NAME |
---|
311 | * |
---|
312 | * bfd_btran - perform backward transformation (solve system B'*x = b) |
---|
313 | * |
---|
314 | * SYNOPSIS |
---|
315 | * |
---|
316 | * #include "glpbfd.h" |
---|
317 | * void bfd_btran(BFD *bfd, double x[]); |
---|
318 | * |
---|
319 | * DESCRIPTION |
---|
320 | * |
---|
321 | * The routine bfd_btran performs backward transformation, i.e. solves |
---|
322 | * the system B'*x = b, where B' is a matrix transposed to the basis |
---|
323 | * matrix B, x is the vector of unknowns to be computed, b is the vector |
---|
324 | * of right-hand sides. |
---|
325 | * |
---|
326 | * On entry elements of the vector b should be stored in dense format |
---|
327 | * in locations x[1], ..., x[m], where m is the number of rows. On exit |
---|
328 | * the routine stores elements of the vector x in the same locations. */ |
---|
329 | |
---|
330 | void bfd_btran(BFD *bfd, double x[]) |
---|
331 | { xassert(bfd != NULL); |
---|
332 | xassert(bfd->valid); |
---|
333 | if (bfd->fhv != NULL) |
---|
334 | fhv_btran(bfd->fhv, x); |
---|
335 | else if (bfd->lpf != NULL) |
---|
336 | lpf_btran(bfd->lpf, x); |
---|
337 | else |
---|
338 | xassert(bfd != bfd); |
---|
339 | return; |
---|
340 | } |
---|
341 | |
---|
342 | /*********************************************************************** |
---|
343 | * NAME |
---|
344 | * |
---|
345 | * bfd_update_it - update LP basis factorization |
---|
346 | * |
---|
347 | * SYNOPSIS |
---|
348 | * |
---|
349 | * #include "glpbfd.h" |
---|
350 | * int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[], |
---|
351 | * const double val[]); |
---|
352 | * |
---|
353 | * DESCRIPTION |
---|
354 | * |
---|
355 | * The routine bfd_update_it updates the factorization of the basis |
---|
356 | * matrix B after replacing its j-th column by a new vector. |
---|
357 | * |
---|
358 | * The parameter j specifies the number of column of B, which has been |
---|
359 | * replaced, 1 <= j <= m, where m is the order of B. |
---|
360 | * |
---|
361 | * The parameter bh specifies the basis header entry for the new column |
---|
362 | * of B, which is the number of the new column in some original matrix. |
---|
363 | * This parameter is optional and can be specified as 0. |
---|
364 | * |
---|
365 | * Row indices and numerical values of non-zero elements of the new |
---|
366 | * column of B should be placed in locations ind[1], ..., ind[len] and |
---|
367 | * val[1], ..., val[len], resp., where len is the number of non-zeros |
---|
368 | * in the column. Neither zero nor duplicate elements are allowed. |
---|
369 | * |
---|
370 | * RETURNS |
---|
371 | * |
---|
372 | * 0 The factorization has been successfully updated. |
---|
373 | * |
---|
374 | * BFD_ESING |
---|
375 | * New basis matrix is singular within the working precision. |
---|
376 | * |
---|
377 | * BFD_ECHECK |
---|
378 | * The factorization is inaccurate. |
---|
379 | * |
---|
380 | * BFD_ELIMIT |
---|
381 | * Factorization update limit has been reached. |
---|
382 | * |
---|
383 | * BFD_EROOM |
---|
384 | * Overflow of the sparse vector area. |
---|
385 | * |
---|
386 | * In case of non-zero return code the factorization becomes invalid. |
---|
387 | * It should not be used until it has been recomputed with the routine |
---|
388 | * bfd_factorize. */ |
---|
389 | |
---|
390 | int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[], |
---|
391 | const double val[]) |
---|
392 | { int ret; |
---|
393 | xassert(bfd != NULL); |
---|
394 | xassert(bfd->valid); |
---|
395 | /* try to update the factorization */ |
---|
396 | if (bfd->fhv != NULL) |
---|
397 | { switch (fhv_update_it(bfd->fhv, j, len, ind, val)) |
---|
398 | { case 0: |
---|
399 | break; |
---|
400 | case FHV_ESING: |
---|
401 | bfd->valid = 0; |
---|
402 | ret = BFD_ESING; |
---|
403 | goto done; |
---|
404 | case FHV_ECHECK: |
---|
405 | bfd->valid = 0; |
---|
406 | ret = BFD_ECHECK; |
---|
407 | goto done; |
---|
408 | case FHV_ELIMIT: |
---|
409 | bfd->valid = 0; |
---|
410 | ret = BFD_ELIMIT; |
---|
411 | goto done; |
---|
412 | case FHV_EROOM: |
---|
413 | bfd->valid = 0; |
---|
414 | ret = BFD_EROOM; |
---|
415 | goto done; |
---|
416 | default: |
---|
417 | xassert(bfd != bfd); |
---|
418 | } |
---|
419 | } |
---|
420 | else if (bfd->lpf != NULL) |
---|
421 | { switch (lpf_update_it(bfd->lpf, j, bh, len, ind, val)) |
---|
422 | { case 0: |
---|
423 | break; |
---|
424 | case LPF_ESING: |
---|
425 | bfd->valid = 0; |
---|
426 | ret = BFD_ESING; |
---|
427 | goto done; |
---|
428 | case LPF_ELIMIT: |
---|
429 | bfd->valid = 0; |
---|
430 | ret = BFD_ELIMIT; |
---|
431 | goto done; |
---|
432 | default: |
---|
433 | xassert(bfd != bfd); |
---|
434 | } |
---|
435 | } |
---|
436 | else |
---|
437 | xassert(bfd != bfd); |
---|
438 | /* the factorization has been successfully updated */ |
---|
439 | /* increase the update count */ |
---|
440 | bfd->upd_cnt++; |
---|
441 | ret = 0; |
---|
442 | done: /* return to the calling program */ |
---|
443 | return ret; |
---|
444 | } |
---|
445 | |
---|
446 | /**********************************************************************/ |
---|
447 | |
---|
448 | int bfd_get_count(BFD *bfd) |
---|
449 | { /* determine factorization update count */ |
---|
450 | xassert(bfd != NULL); |
---|
451 | xassert(bfd->valid); |
---|
452 | return bfd->upd_cnt; |
---|
453 | } |
---|
454 | |
---|
455 | /*********************************************************************** |
---|
456 | * NAME |
---|
457 | * |
---|
458 | * bfd_delete_it - delete LP basis factorization |
---|
459 | * |
---|
460 | * SYNOPSIS |
---|
461 | * |
---|
462 | * #include "glpbfd.h" |
---|
463 | * void bfd_delete_it(BFD *bfd); |
---|
464 | * |
---|
465 | * DESCRIPTION |
---|
466 | * |
---|
467 | * The routine bfd_delete_it deletes LP basis factorization specified |
---|
468 | * by the parameter fhv and frees all memory allocated to this program |
---|
469 | * object. */ |
---|
470 | |
---|
471 | void bfd_delete_it(BFD *bfd) |
---|
472 | { xassert(bfd != NULL); |
---|
473 | if (bfd->fhv != NULL) |
---|
474 | fhv_delete_it(bfd->fhv); |
---|
475 | if (bfd->lpf != NULL) |
---|
476 | lpf_delete_it(bfd->lpf); |
---|
477 | xfree(bfd); |
---|
478 | return; |
---|
479 | } |
---|
480 | |
---|
481 | /* eof */ |
---|