lemon-project-template-glpk
view deps/glpk/src/glpbfd.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line source
1 /* glpbfd.c (LP basis factorization driver) */
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, 2011 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 ***********************************************************************/
25 typedef struct BFD BFD;
27 #define GLPBFD_PRIVATE
28 #include "glpapi.h"
29 #include "glpfhv.h"
30 #include "glplpf.h"
32 /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
34 #define M_MAX 100000000 /* = 100*10^6 */
35 /* maximal order of the basis matrix */
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 };
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. */
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 }
108 /**********************************************************************/
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
446 /**********************************************************************/
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 }
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. */
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 }
481 /* eof */