diff -r d59bea55db9b -r c445c931472f src/glpbfd.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/glpbfd.c Mon Dec 06 13:09:21 2010 +0100 @@ -0,0 +1,481 @@ +/* glpbfd.c (LP basis factorization driver) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics, +* Moscow Aviation Institute, Moscow, Russia. All rights reserved. +* E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +typedef struct BFD BFD; + +#define GLPBFD_PRIVATE +#include "glpapi.h" +#include "glpfhv.h" +#include "glplpf.h" + +/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ + +#define M_MAX 100000000 /* = 100*10^6 */ +/* maximal order of the basis matrix */ + +struct BFD +{ /* LP basis factorization */ + int valid; + /* factorization is valid only if this flag is set */ + int type; + /* factorization type: + GLP_BF_FT - LUF + Forrest-Tomlin + GLP_BF_BG - LUF + Schur compl. + Bartels-Golub + GLP_BF_GR - LUF + Schur compl. + Givens rotation */ + FHV *fhv; + /* LP basis factorization (GLP_BF_FT) */ + LPF *lpf; + /* LP basis factorization (GLP_BF_BG, GLP_BF_GR) */ + int lu_size; /* luf.sv_size */ + double piv_tol; /* luf.piv_tol */ + int piv_lim; /* luf.piv_lim */ + int suhl; /* luf.suhl */ + double eps_tol; /* luf.eps_tol */ + double max_gro; /* luf.max_gro */ + int nfs_max; /* fhv.hh_max */ + double upd_tol; /* fhv.upd_tol */ + int nrs_max; /* lpf.n_max */ + int rs_size; /* lpf.v_size */ + /* internal control parameters */ + int upd_lim; + /* the factorization update limit */ + int upd_cnt; + /* the factorization update count */ +}; + +/*********************************************************************** +* NAME +* +* bfd_create_it - create LP basis factorization +* +* SYNOPSIS +* +* #include "glpbfd.h" +* BFD *bfd_create_it(void); +* +* DESCRIPTION +* +* The routine bfd_create_it creates a program object, which represents +* a factorization of LP basis. +* +* RETURNS +* +* The routine bfd_create_it returns a pointer to the object created. */ + +BFD *bfd_create_it(void) +{ BFD *bfd; + bfd = xmalloc(sizeof(BFD)); + bfd->valid = 0; + bfd->type = GLP_BF_FT; + bfd->fhv = NULL; + bfd->lpf = NULL; + bfd->lu_size = 0; + bfd->piv_tol = 0.10; + bfd->piv_lim = 4; + bfd->suhl = 1; + bfd->eps_tol = 1e-15; + bfd->max_gro = 1e+10; + bfd->nfs_max = 100; + bfd->upd_tol = 1e-6; + bfd->nrs_max = 100; + bfd->rs_size = 1000; + bfd->upd_lim = -1; + bfd->upd_cnt = 0; + return bfd; +} + +/**********************************************************************/ + +void bfd_set_parm(BFD *bfd, const void *_parm) +{ /* change LP basis factorization control parameters */ + const glp_bfcp *parm = _parm; + xassert(bfd != NULL); + bfd->type = parm->type; + bfd->lu_size = parm->lu_size; + bfd->piv_tol = parm->piv_tol; + bfd->piv_lim = parm->piv_lim; + bfd->suhl = parm->suhl; + bfd->eps_tol = parm->eps_tol; + bfd->max_gro = parm->max_gro; + bfd->nfs_max = parm->nfs_max; + bfd->upd_tol = parm->upd_tol; + bfd->nrs_max = parm->nrs_max; + bfd->rs_size = parm->rs_size; + return; +} + +/*********************************************************************** +* NAME +* +* bfd_factorize - compute LP basis factorization +* +* SYNOPSIS +* +* #include "glpbfd.h" +* int bfd_factorize(BFD *bfd, int m, int bh[], int (*col)(void *info, +* int j, int ind[], double val[]), void *info); +* +* DESCRIPTION +* +* The routine bfd_factorize computes the factorization of the basis +* matrix B specified by the routine col. +* +* The parameter bfd specified the basis factorization data structure +* created with the routine bfd_create_it. +* +* The parameter m specifies the order of B, m > 0. +* +* The array bh specifies the basis header: bh[j], 1 <= j <= m, is the +* number of j-th column of B in some original matrix. The array bh is +* optional and can be specified as NULL. +* +* The formal routine col specifies the matrix B to be factorized. To +* obtain j-th column of A the routine bfd_factorize calls the routine +* col with the parameter j (1 <= j <= n). In response the routine col +* should store row indices and numerical values of non-zero elements +* of j-th column of B to locations ind[1,...,len] and val[1,...,len], +* respectively, where len is the number of non-zeros in j-th column +* returned on exit. Neither zero nor duplicate elements are allowed. +* +* The parameter info is a transit pointer passed to the routine col. +* +* RETURNS +* +* 0 The factorization has been successfully computed. +* +* BFD_ESING +* The specified matrix is singular within the working precision. +* +* BFD_ECOND +* The specified matrix is ill-conditioned. +* +* For more details see comments to the routine luf_factorize. */ + +int bfd_factorize(BFD *bfd, int m, const int bh[], int (*col) + (void *info, int j, int ind[], double val[]), void *info) +{ LUF *luf; + int nov, ret; + xassert(bfd != NULL); + xassert(1 <= m && m <= M_MAX); + /* invalidate the factorization */ + bfd->valid = 0; + /* create the factorization, if necessary */ + nov = 0; + switch (bfd->type) + { case GLP_BF_FT: + if (bfd->lpf != NULL) + lpf_delete_it(bfd->lpf), bfd->lpf = NULL; + if (bfd->fhv == NULL) + bfd->fhv = fhv_create_it(), nov = 1; + break; + case GLP_BF_BG: + case GLP_BF_GR: + if (bfd->fhv != NULL) + fhv_delete_it(bfd->fhv), bfd->fhv = NULL; + if (bfd->lpf == NULL) + bfd->lpf = lpf_create_it(), nov = 1; + break; + default: + xassert(bfd != bfd); + } + /* set control parameters specific to LUF */ + if (bfd->fhv != NULL) + luf = bfd->fhv->luf; + else if (bfd->lpf != NULL) + luf = bfd->lpf->luf; + else + xassert(bfd != bfd); + if (nov) luf->new_sva = bfd->lu_size; + luf->piv_tol = bfd->piv_tol; + luf->piv_lim = bfd->piv_lim; + luf->suhl = bfd->suhl; + luf->eps_tol = bfd->eps_tol; + luf->max_gro = bfd->max_gro; + /* set control parameters specific to FHV */ + if (bfd->fhv != NULL) + { if (nov) bfd->fhv->hh_max = bfd->nfs_max; + bfd->fhv->upd_tol = bfd->upd_tol; + } + /* set control parameters specific to LPF */ + if (bfd->lpf != NULL) + { if (nov) bfd->lpf->n_max = bfd->nrs_max; + if (nov) bfd->lpf->v_size = bfd->rs_size; + } + /* try to factorize the basis matrix */ + if (bfd->fhv != NULL) + { switch (fhv_factorize(bfd->fhv, m, col, info)) + { case 0: + break; + case FHV_ESING: + ret = BFD_ESING; + goto done; + case FHV_ECOND: + ret = BFD_ECOND; + goto done; + default: + xassert(bfd != bfd); + } + } + else if (bfd->lpf != NULL) + { switch (lpf_factorize(bfd->lpf, m, bh, col, info)) + { case 0: + /* set the Schur complement update type */ + switch (bfd->type) + { case GLP_BF_BG: + /* Bartels-Golub update */ + bfd->lpf->scf->t_opt = SCF_TBG; + break; + case GLP_BF_GR: + /* Givens rotation update */ + bfd->lpf->scf->t_opt = SCF_TGR; + break; + default: + xassert(bfd != bfd); + } + break; + case LPF_ESING: + ret = BFD_ESING; + goto done; + case LPF_ECOND: + ret = BFD_ECOND; + goto done; + default: + xassert(bfd != bfd); + } + } + else + xassert(bfd != bfd); + /* the basis matrix has been successfully factorized */ + bfd->valid = 1; + bfd->upd_cnt = 0; + ret = 0; +done: /* return to the calling program */ + return ret; +} + +/*********************************************************************** +* NAME +* +* bfd_ftran - perform forward transformation (solve system B*x = b) +* +* SYNOPSIS +* +* #include "glpbfd.h" +* void bfd_ftran(BFD *bfd, double x[]); +* +* DESCRIPTION +* +* The routine bfd_ftran performs forward transformation, i.e. solves +* the system B*x = b, where B is the basis matrix, x is the vector of +* unknowns to be computed, b is the vector of right-hand sides. +* +* On entry elements of the vector b should be stored in dense format +* in locations x[1], ..., x[m], where m is the number of rows. On exit +* the routine stores elements of the vector x in the same locations. */ + +void bfd_ftran(BFD *bfd, double x[]) +{ xassert(bfd != NULL); + xassert(bfd->valid); + if (bfd->fhv != NULL) + fhv_ftran(bfd->fhv, x); + else if (bfd->lpf != NULL) + lpf_ftran(bfd->lpf, x); + else + xassert(bfd != bfd); + return; +} + +/*********************************************************************** +* NAME +* +* bfd_btran - perform backward transformation (solve system B'*x = b) +* +* SYNOPSIS +* +* #include "glpbfd.h" +* void bfd_btran(BFD *bfd, double x[]); +* +* DESCRIPTION +* +* The routine bfd_btran performs backward transformation, i.e. solves +* the system B'*x = b, where B' is a matrix transposed to the basis +* matrix B, x is the vector of unknowns to be computed, b is the vector +* of right-hand sides. +* +* On entry elements of the vector b should be stored in dense format +* in locations x[1], ..., x[m], where m is the number of rows. On exit +* the routine stores elements of the vector x in the same locations. */ + +void bfd_btran(BFD *bfd, double x[]) +{ xassert(bfd != NULL); + xassert(bfd->valid); + if (bfd->fhv != NULL) + fhv_btran(bfd->fhv, x); + else if (bfd->lpf != NULL) + lpf_btran(bfd->lpf, x); + else + xassert(bfd != bfd); + return; +} + +/*********************************************************************** +* NAME +* +* bfd_update_it - update LP basis factorization +* +* SYNOPSIS +* +* #include "glpbfd.h" +* int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[], +* const double val[]); +* +* DESCRIPTION +* +* The routine bfd_update_it updates the factorization of the basis +* matrix B after replacing its j-th column by a new vector. +* +* The parameter j specifies the number of column of B, which has been +* replaced, 1 <= j <= m, where m is the order of B. +* +* The parameter bh specifies the basis header entry for the new column +* of B, which is the number of the new column in some original matrix. +* This parameter is optional and can be specified as 0. +* +* Row indices and numerical values of non-zero elements of the new +* column of B should be placed in locations ind[1], ..., ind[len] and +* val[1], ..., val[len], resp., where len is the number of non-zeros +* in the column. Neither zero nor duplicate elements are allowed. +* +* RETURNS +* +* 0 The factorization has been successfully updated. +* +* BFD_ESING +* New basis matrix is singular within the working precision. +* +* BFD_ECHECK +* The factorization is inaccurate. +* +* BFD_ELIMIT +* Factorization update limit has been reached. +* +* BFD_EROOM +* Overflow of the sparse vector area. +* +* In case of non-zero return code the factorization becomes invalid. +* It should not be used until it has been recomputed with the routine +* bfd_factorize. */ + +int bfd_update_it(BFD *bfd, int j, int bh, int len, const int ind[], + const double val[]) +{ int ret; + xassert(bfd != NULL); + xassert(bfd->valid); + /* try to update the factorization */ + if (bfd->fhv != NULL) + { switch (fhv_update_it(bfd->fhv, j, len, ind, val)) + { case 0: + break; + case FHV_ESING: + bfd->valid = 0; + ret = BFD_ESING; + goto done; + case FHV_ECHECK: + bfd->valid = 0; + ret = BFD_ECHECK; + goto done; + case FHV_ELIMIT: + bfd->valid = 0; + ret = BFD_ELIMIT; + goto done; + case FHV_EROOM: + bfd->valid = 0; + ret = BFD_EROOM; + goto done; + default: + xassert(bfd != bfd); + } + } + else if (bfd->lpf != NULL) + { switch (lpf_update_it(bfd->lpf, j, bh, len, ind, val)) + { case 0: + break; + case LPF_ESING: + bfd->valid = 0; + ret = BFD_ESING; + goto done; + case LPF_ELIMIT: + bfd->valid = 0; + ret = BFD_ELIMIT; + goto done; + default: + xassert(bfd != bfd); + } + } + else + xassert(bfd != bfd); + /* the factorization has been successfully updated */ + /* increase the update count */ + bfd->upd_cnt++; + ret = 0; +done: /* return to the calling program */ + return ret; +} + +/**********************************************************************/ + +int bfd_get_count(BFD *bfd) +{ /* determine factorization update count */ + xassert(bfd != NULL); + xassert(bfd->valid); + return bfd->upd_cnt; +} + +/*********************************************************************** +* NAME +* +* bfd_delete_it - delete LP basis factorization +* +* SYNOPSIS +* +* #include "glpbfd.h" +* void bfd_delete_it(BFD *bfd); +* +* DESCRIPTION +* +* The routine bfd_delete_it deletes LP basis factorization specified +* by the parameter fhv and frees all memory allocated to this program +* object. */ + +void bfd_delete_it(BFD *bfd) +{ xassert(bfd != NULL); + if (bfd->fhv != NULL) + fhv_delete_it(bfd->fhv); + if (bfd->lpf != NULL) + lpf_delete_it(bfd->lpf); + xfree(bfd); + return; +} + +/* eof */