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