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