/* 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: <mao@gnu.org>.

 *

 *  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 <http://www.gnu.org/licenses/>.

 ***********************************************************************/

 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 */