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