praat/external/glpk/glpbfd.c

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