COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpbfd.c @ 2:4c8956a7bdf4

Last change on this file since 2:4c8956a7bdf4 was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

  • Generated files and doc/notes are removed
File size: 14.3 KB
Line 
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
25typedef 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
37struct 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
86BFD *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
110void 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
175int 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;
273done: /* 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
297void 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
330void 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
390int 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;
442done: /* return to the calling program */
443      return ret;
444}
445
446/**********************************************************************/
447
448int 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
471void 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 */
Note: See TracBrowser for help on using the repository browser.