|
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 */ |