1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpbfd.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,481 @@
     1.4 +/* glpbfd.c (LP basis factorization driver) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +typedef struct BFD BFD;
    1.29 +
    1.30 +#define GLPBFD_PRIVATE
    1.31 +#include "glpapi.h"
    1.32 +#include "glpfhv.h"
    1.33 +#include "glplpf.h"
    1.34 +
    1.35 +/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
    1.36 +
    1.37 +#define M_MAX 100000000 /* = 100*10^6 */
    1.38 +/* maximal order of the basis matrix */
    1.39 +
    1.40 +struct BFD
    1.41 +{     /* LP basis factorization */
    1.42 +      int valid;
    1.43 +      /* factorization is valid only if this flag is set */
    1.44 +      int type;
    1.45 +      /* factorization type:
    1.46 +         GLP_BF_FT - LUF + Forrest-Tomlin
    1.47 +         GLP_BF_BG - LUF + Schur compl. + Bartels-Golub
    1.48 +         GLP_BF_GR - LUF + Schur compl. + Givens rotation */
    1.49 +      FHV *fhv;
    1.50 +      /* LP basis factorization (GLP_BF_FT) */
    1.51 +      LPF *lpf;
    1.52 +      /* LP basis factorization (GLP_BF_BG, GLP_BF_GR) */
    1.53 +      int lu_size;      /* luf.sv_size */
    1.54 +      double piv_tol;   /* luf.piv_tol */
    1.55 +      int piv_lim;      /* luf.piv_lim */
    1.56 +      int suhl;         /* luf.suhl */
    1.57 +      double eps_tol;   /* luf.eps_tol */
    1.58 +      double max_gro;   /* luf.max_gro */
    1.59 +      int nfs_max;      /* fhv.hh_max */
    1.60 +      double upd_tol;   /* fhv.upd_tol */
    1.61 +      int nrs_max;      /* lpf.n_max */
    1.62 +      int rs_size;      /* lpf.v_size */
    1.63 +      /* internal control parameters */
    1.64 +      int upd_lim;
    1.65 +      /* the factorization update limit */
    1.66 +      int upd_cnt;
    1.67 +      /* the factorization update count */
    1.68 +};
    1.69 +
    1.70 +/***********************************************************************
    1.71 +*  NAME
    1.72 +*
    1.73 +*  bfd_create_it - create LP basis factorization
    1.74 +*
    1.75 +*  SYNOPSIS
    1.76 +*
    1.77 +*  #include "glpbfd.h"
    1.78 +*  BFD *bfd_create_it(void);
    1.79 +*
    1.80 +*  DESCRIPTION
    1.81 +*
    1.82 +*  The routine bfd_create_it creates a program object, which represents
    1.83 +*  a factorization of LP basis.
    1.84 +*
    1.85 +*  RETURNS
    1.86 +*
    1.87 +*  The routine bfd_create_it returns a pointer to the object created. */
    1.88 +
    1.89 +BFD *bfd_create_it(void)
    1.90 +{     BFD *bfd;
    1.91 +      bfd = xmalloc(sizeof(BFD));
    1.92 +      bfd->valid = 0;
    1.93 +      bfd->type = GLP_BF_FT;
    1.94 +      bfd->fhv = NULL;
    1.95 +      bfd->lpf = NULL;
    1.96 +      bfd->lu_size = 0;
    1.97 +      bfd->piv_tol = 0.10;
    1.98 +      bfd->piv_lim = 4;
    1.99 +      bfd->suhl = 1;
   1.100 +      bfd->eps_tol = 1e-15;
   1.101 +      bfd->max_gro = 1e+10;
   1.102 +      bfd->nfs_max = 100;
   1.103 +      bfd->upd_tol = 1e-6;
   1.104 +      bfd->nrs_max = 100;
   1.105 +      bfd->rs_size = 1000;
   1.106 +      bfd->upd_lim = -1;
   1.107 +      bfd->upd_cnt = 0;
   1.108 +      return bfd;
   1.109 +}
   1.110 +
   1.111 +/**********************************************************************/
   1.112 +
   1.113 +void bfd_set_parm(BFD *bfd, const void *_parm)
   1.114 +{     /* change LP basis factorization control parameters */
   1.115 +      const glp_bfcp *parm = _parm;
   1.116 +      xassert(bfd != NULL);
   1.117 +      bfd->type = parm->type;
   1.118 +      bfd->lu_size = parm->lu_size;
   1.119 +      bfd->piv_tol = parm->piv_tol;
   1.120 +      bfd->piv_lim = parm->piv_lim;
   1.121 +      bfd->suhl = parm->suhl;
   1.122 +      bfd->eps_tol = parm->eps_tol;
   1.123 +      bfd->max_gro = parm->max_gro;
   1.124 +      bfd->nfs_max = parm->nfs_max;
   1.125 +      bfd->upd_tol = parm->upd_tol;
   1.126 +      bfd->nrs_max = parm->nrs_max;
   1.127 +      bfd->rs_size = parm->rs_size;
   1.128 +      return;
   1.129 +}
   1.130 +
   1.131 +/***********************************************************************
   1.132 +*  NAME
   1.133 +*
   1.134 +*  bfd_factorize - compute LP basis factorization
   1.135 +*
   1.136 +*  SYNOPSIS
   1.137 +*
   1.138 +*  #include "glpbfd.h"
   1.139 +*  int bfd_factorize(BFD *bfd, int m, int bh[], int (*col)(void *info,
   1.140 +*     int j, int ind[], double val[]), void *info);
   1.141 +*
   1.142 +*  DESCRIPTION
   1.143 +*
   1.144 +*  The routine bfd_factorize computes the factorization of the basis
   1.145 +*  matrix B specified by the routine col.
   1.146 +*
   1.147 +*  The parameter bfd specified the basis factorization data structure
   1.148 +*  created with the routine bfd_create_it.
   1.149 +*
   1.150 +*  The parameter m specifies the order of B, m > 0.
   1.151 +*
   1.152 +*  The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
   1.153 +*  number of j-th column of B in some original matrix. The array bh is
   1.154 +*  optional and can be specified as NULL.
   1.155 +*
   1.156 +*  The formal routine col specifies the matrix B to be factorized. To
   1.157 +*  obtain j-th column of A the routine bfd_factorize calls the routine
   1.158 +*  col with the parameter j (1 <= j <= n). In response the routine col
   1.159 +*  should store row indices and numerical values of non-zero elements
   1.160 +*  of j-th column of B to locations ind[1,...,len] and val[1,...,len],
   1.161 +*  respectively, where len is the number of non-zeros in j-th column
   1.162 +*  returned on exit. Neither zero nor duplicate elements are allowed.
   1.163 +*
   1.164 +*  The parameter info is a transit pointer passed to the routine col.
   1.165 +*
   1.166 +*  RETURNS
   1.167 +*
   1.168 +*  0  The factorization has been successfully computed.
   1.169 +*
   1.170 +*  BFD_ESING
   1.171 +*     The specified matrix is singular within the working precision.
   1.172 +*
   1.173 +*  BFD_ECOND
   1.174 +*     The specified matrix is ill-conditioned.
   1.175 +*
   1.176 +*  For more details see comments to the routine luf_factorize. */
   1.177 +
   1.178 +int bfd_factorize(BFD *bfd, int m, const int bh[], int (*col)
   1.179 +      (void *info, int j, int ind[], double val[]), void *info)
   1.180 +{     LUF *luf;
   1.181 +      int nov, ret;
   1.182 +      xassert(bfd != NULL);
   1.183 +      xassert(1 <= m && m <= M_MAX);
   1.184 +      /* invalidate the factorization */
   1.185 +      bfd->valid = 0;
   1.186 +      /* create the factorization, if necessary */
   1.187 +      nov = 0;
   1.188 +      switch (bfd->type)
   1.189 +      {  case GLP_BF_FT:
   1.190 +            if (bfd->lpf != NULL)
   1.191 +               lpf_delete_it(bfd->lpf), bfd->lpf = NULL;
   1.192 +            if (bfd->fhv == NULL)
   1.193 +               bfd->fhv = fhv_create_it(), nov = 1;
   1.194 +            break;
   1.195 +         case GLP_BF_BG:
   1.196 +         case GLP_BF_GR:
   1.197 +            if (bfd->fhv != NULL)
   1.198 +               fhv_delete_it(bfd->fhv), bfd->fhv = NULL;
   1.199 +            if (bfd->lpf == NULL)
   1.200 +               bfd->lpf = lpf_create_it(), nov = 1;
   1.201 +            break;
   1.202 +         default:
   1.203 +            xassert(bfd != bfd);
   1.204 +      }
   1.205 +      /* set control parameters specific to LUF */
   1.206 +      if (bfd->fhv != NULL)
   1.207 +         luf = bfd->fhv->luf;
   1.208 +      else if (bfd->lpf != NULL)
   1.209 +         luf = bfd->lpf->luf;
   1.210 +      else
   1.211 +         xassert(bfd != bfd);
   1.212 +      if (nov) luf->new_sva = bfd->lu_size;
   1.213 +      luf->piv_tol = bfd->piv_tol;
   1.214 +      luf->piv_lim = bfd->piv_lim;
   1.215 +      luf->suhl = bfd->suhl;
   1.216 +      luf->eps_tol = bfd->eps_tol;
   1.217 +      luf->max_gro = bfd->max_gro;
   1.218 +      /* set control parameters specific to FHV */
   1.219 +      if (bfd->fhv != NULL)
   1.220 +      {  if (nov) bfd->fhv->hh_max = bfd->nfs_max;
   1.221 +         bfd->fhv->upd_tol = bfd->upd_tol;
   1.222 +      }
   1.223 +      /* set control parameters specific to LPF */
   1.224 +      if (bfd->lpf != NULL)
   1.225 +      {  if (nov) bfd->lpf->n_max = bfd->nrs_max;
   1.226 +         if (nov) bfd->lpf->v_size = bfd->rs_size;
   1.227 +      }
   1.228 +      /* try to factorize the basis matrix */
   1.229 +      if (bfd->fhv != NULL)
   1.230 +      {  switch (fhv_factorize(bfd->fhv, m, col, info))
   1.231 +         {  case 0:
   1.232 +               break;
   1.233 +            case FHV_ESING:
   1.234 +               ret = BFD_ESING;
   1.235 +               goto done;
   1.236 +            case FHV_ECOND:
   1.237 +               ret = BFD_ECOND;
   1.238 +               goto done;
   1.239 +            default:
   1.240 +               xassert(bfd != bfd);
   1.241 +         }
   1.242 +      }
   1.243 +      else if (bfd->lpf != NULL)
   1.244 +      {  switch (lpf_factorize(bfd->lpf, m, bh, col, info))
   1.245 +         {  case 0:
   1.246 +               /* set the Schur complement update type */
   1.247 +               switch (bfd->type)
   1.248 +               {  case GLP_BF_BG:
   1.249 +                     /* Bartels-Golub update */
   1.250 +                     bfd->lpf->scf->t_opt = SCF_TBG;
   1.251 +                     break;
   1.252 +                  case GLP_BF_GR:
   1.253 +                     /* Givens rotation update */
   1.254 +                     bfd->lpf->scf->t_opt = SCF_TGR;
   1.255 +                     break;
   1.256 +                  default:
   1.257 +                     xassert(bfd != bfd);
   1.258 +               }
   1.259 +               break;
   1.260 +            case LPF_ESING:
   1.261 +               ret = BFD_ESING;
   1.262 +               goto done;
   1.263 +            case LPF_ECOND:
   1.264 +               ret = BFD_ECOND;
   1.265 +               goto done;
   1.266 +            default:
   1.267 +               xassert(bfd != bfd);
   1.268 +         }
   1.269 +      }
   1.270 +      else
   1.271 +         xassert(bfd != bfd);
   1.272 +      /* the basis matrix has been successfully factorized */
   1.273 +      bfd->valid = 1;
   1.274 +      bfd->upd_cnt = 0;
   1.275 +      ret = 0;
   1.276 +done: /* return to the calling program */
   1.277 +      return ret;
   1.278 +}
   1.279 +
   1.280 +/***********************************************************************
   1.281 +*  NAME
   1.282 +*
   1.283 +*  bfd_ftran - perform forward transformation (solve system B*x = b)
   1.284 +*
   1.285 +*  SYNOPSIS
   1.286 +*
   1.287 +*  #include "glpbfd.h"
   1.288 +*  void bfd_ftran(BFD *bfd, double x[]);
   1.289 +*
   1.290 +*  DESCRIPTION
   1.291 +*
   1.292 +*  The routine bfd_ftran performs forward transformation, i.e. solves
   1.293 +*  the system B*x = b, where B is the basis matrix, x is the vector of
   1.294 +*  unknowns to be computed, b is the vector of right-hand sides.
   1.295 +*
   1.296 +*  On entry elements of the vector b should be stored in dense format
   1.297 +*  in locations x[1], ..., x[m], where m is the number of rows. On exit
   1.298 +*  the routine stores elements of the vector x in the same locations. */
   1.299 +
   1.300 +void bfd_ftran(BFD *bfd, double x[])
   1.301 +{     xassert(bfd != NULL);
   1.302 +      xassert(bfd->valid);
   1.303 +      if (bfd->fhv != NULL)
   1.304 +         fhv_ftran(bfd->fhv, x);
   1.305 +      else if (bfd->lpf != NULL)
   1.306 +         lpf_ftran(bfd->lpf, x);
   1.307 +      else
   1.308 +         xassert(bfd != bfd);
   1.309 +      return;
   1.310 +}
   1.311 +
   1.312 +/***********************************************************************
   1.313 +*  NAME
   1.314 +*
   1.315 +*  bfd_btran - perform backward transformation (solve system B'*x = b)
   1.316 +*
   1.317 +*  SYNOPSIS
   1.318 +*
   1.319 +*  #include "glpbfd.h"
   1.320 +*  void bfd_btran(BFD *bfd, double x[]);
   1.321 +*
   1.322 +*  DESCRIPTION
   1.323 +*
   1.324 +*  The routine bfd_btran performs backward transformation, i.e. solves
   1.325 +*  the system B'*x = b, where B' is a matrix transposed to the basis
   1.326 +*  matrix B, x is the vector of unknowns to be computed, b is the vector
   1.327 +*  of right-hand sides.
   1.328 +*
   1.329 +*  On entry elements of the vector b should be stored in dense format
   1.330 +*  in locations x[1], ..., x[m], where m is the number of rows. On exit
   1.331 +*  the routine stores elements of the vector x in the same locations. */
   1.332 +
   1.333 +void bfd_btran(BFD *bfd, double x[])
   1.334 +{     xassert(bfd != NULL);
   1.335 +      xassert(bfd->valid);
   1.336 +      if (bfd->fhv != NULL)
   1.337 +         fhv_btran(bfd->fhv, x);
   1.338 +      else if (bfd->lpf != NULL)
   1.339 +         lpf_btran(bfd->lpf, x);
   1.340 +      else
   1.341 +         xassert(bfd != bfd);
   1.342 +      return;
   1.343 +}
   1.344 +
   1.345 +/***********************************************************************
   1.346 +*  NAME
   1.347 +*
   1.348 +*  bfd_update_it - update LP basis factorization
   1.349 +*
   1.350 +*  SYNOPSIS
   1.351 +*
   1.352 +*  #include "glpbfd.h"
   1.353 +*  int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[],
   1.354 +*     const double val[]);
   1.355 +*
   1.356 +*  DESCRIPTION
   1.357 +*
   1.358 +*  The routine bfd_update_it updates the factorization of the basis
   1.359 +*  matrix B after replacing its j-th column by a new vector.
   1.360 +*
   1.361 +*  The parameter j specifies the number of column of B, which has been
   1.362 +*  replaced, 1 <= j <= m, where m is the order of B.
   1.363 +*
   1.364 +*  The parameter bh specifies the basis header entry for the new column
   1.365 +*  of B, which is the number of the new column in some original matrix.
   1.366 +*  This parameter is optional and can be specified as 0.
   1.367 +*
   1.368 +*  Row indices and numerical values of non-zero elements of the new
   1.369 +*  column of B should be placed in locations ind[1], ..., ind[len] and
   1.370 +*  val[1], ..., val[len], resp., where len is the number of non-zeros
   1.371 +*  in the column. Neither zero nor duplicate elements are allowed.
   1.372 +*
   1.373 +*  RETURNS
   1.374 +*
   1.375 +*  0  The factorization has been successfully updated.
   1.376 +*
   1.377 +*  BFD_ESING
   1.378 +*     New basis matrix is singular within the working precision.
   1.379 +*
   1.380 +*  BFD_ECHECK
   1.381 +*     The factorization is inaccurate.
   1.382 +*
   1.383 +*  BFD_ELIMIT
   1.384 +*     Factorization update limit has been reached.
   1.385 +*
   1.386 +*  BFD_EROOM
   1.387 +*     Overflow of the sparse vector area.
   1.388 +*
   1.389 +*  In case of non-zero return code the factorization becomes invalid.
   1.390 +*  It should not be used until it has been recomputed with the routine
   1.391 +*  bfd_factorize. */
   1.392 +
   1.393 +int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[],
   1.394 +      const double val[])
   1.395 +{     int ret;
   1.396 +      xassert(bfd != NULL);
   1.397 +      xassert(bfd->valid);
   1.398 +      /* try to update the factorization */
   1.399 +      if (bfd->fhv != NULL)
   1.400 +      {  switch (fhv_update_it(bfd->fhv, j, len, ind, val))
   1.401 +         {  case 0:
   1.402 +               break;
   1.403 +            case FHV_ESING:
   1.404 +               bfd->valid = 0;
   1.405 +               ret = BFD_ESING;
   1.406 +               goto done;
   1.407 +            case FHV_ECHECK:
   1.408 +               bfd->valid = 0;
   1.409 +               ret = BFD_ECHECK;
   1.410 +               goto done;
   1.411 +            case FHV_ELIMIT:
   1.412 +               bfd->valid = 0;
   1.413 +               ret = BFD_ELIMIT;
   1.414 +               goto done;
   1.415 +            case FHV_EROOM:
   1.416 +               bfd->valid = 0;
   1.417 +               ret = BFD_EROOM;
   1.418 +               goto done;
   1.419 +            default:
   1.420 +               xassert(bfd != bfd);
   1.421 +         }
   1.422 +      }
   1.423 +      else if (bfd->lpf != NULL)
   1.424 +      {  switch (lpf_update_it(bfd->lpf, j, bh, len, ind, val))
   1.425 +         {  case 0:
   1.426 +               break;
   1.427 +            case LPF_ESING:
   1.428 +               bfd->valid = 0;
   1.429 +               ret = BFD_ESING;
   1.430 +               goto done;
   1.431 +            case LPF_ELIMIT:
   1.432 +               bfd->valid = 0;
   1.433 +               ret = BFD_ELIMIT;
   1.434 +               goto done;
   1.435 +            default:
   1.436 +               xassert(bfd != bfd);
   1.437 +         }
   1.438 +      }
   1.439 +      else
   1.440 +         xassert(bfd != bfd);
   1.441 +      /* the factorization has been successfully updated */
   1.442 +      /* increase the update count */
   1.443 +      bfd->upd_cnt++;
   1.444 +      ret = 0;
   1.445 +done: /* return to the calling program */
   1.446 +      return ret;
   1.447 +}
   1.448 +
   1.449 +/**********************************************************************/
   1.450 +
   1.451 +int bfd_get_count(BFD *bfd)
   1.452 +{     /* determine factorization update count */
   1.453 +      xassert(bfd != NULL);
   1.454 +      xassert(bfd->valid);
   1.455 +      return bfd->upd_cnt;
   1.456 +}
   1.457 +
   1.458 +/***********************************************************************
   1.459 +*  NAME
   1.460 +*
   1.461 +*  bfd_delete_it - delete LP basis factorization
   1.462 +*
   1.463 +*  SYNOPSIS
   1.464 +*
   1.465 +*  #include "glpbfd.h"
   1.466 +*  void bfd_delete_it(BFD *bfd);
   1.467 +*
   1.468 +*  DESCRIPTION
   1.469 +*
   1.470 +*  The routine bfd_delete_it deletes LP basis factorization specified
   1.471 +*  by the parameter fhv and frees all memory allocated to this program
   1.472 +*  object. */
   1.473 +
   1.474 +void bfd_delete_it(BFD *bfd)
   1.475 +{     xassert(bfd != NULL);
   1.476 +      if (bfd->fhv != NULL)
   1.477 +         fhv_delete_it(bfd->fhv);
   1.478 +      if (bfd->lpf != NULL)
   1.479 +         lpf_delete_it(bfd->lpf);
   1.480 +      xfree(bfd);
   1.481 +      return;
   1.482 +}
   1.483 +
   1.484 +/* eof */