/farmR/src/glpbfd.c

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