glpfhv.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

/* glpfhv.c (LP basis factorization, FHV eta file version) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[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 <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpfhv.h"
#include "glplib.h"
#define xfault xerror

/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */

#define M_MAX 100000000 /* = 100*10^6 */
/* maximal order of the basis matrix */

/***********************************************************************
*  NAME
*
*  fhv_create_it - create LP basis factorization
*
*  SYNOPSIS
*
*  #include "glpfhv.h"
*  FHV *fhv_create_it(void);
*
*  DESCRIPTION
*
*  The routine fhv_create_it creates a program object, which represents
*  a factorization of LP basis.
*
*  RETURNS
*
*  The routine fhv_create_it returns a pointer to the object created. */

FHV *fhv_create_it(void)
{     FHV *fhv;
      fhv = xmalloc(sizeof(FHV));
      fhv->m_max = fhv->m = 0;
      fhv->valid = 0;
      fhv->luf = luf_create_it();
      fhv->hh_max = 50;
      fhv->hh_nfs = 0;
      fhv->hh_ind = fhv->hh_ptr = fhv->hh_len = NULL;
      fhv->p0_row = fhv->p0_col = NULL;
      fhv->cc_ind = NULL;
      fhv->cc_val = NULL;
      fhv->upd_tol = 1e-6;
      fhv->nnz_h = 0;
      return fhv;
}

/***********************************************************************
*  NAME
*
*  fhv_factorize - compute LP basis factorization
*
*  SYNOPSIS
*
*  #include "glpfhv.h"
*  int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
*     int ind[], double val[]), void *info);
*
*  DESCRIPTION
*
*  The routine fhv_factorize computes the factorization of the basis
*  matrix B specified by the routine col.
*
*  The parameter fhv specified the basis factorization data structure
*  created by the routine fhv_create_it.
*
*  The parameter m specifies the order of B, m > 0.
*
*  The formal routine col specifies the matrix B to be factorized. To
*  obtain j-th column of A the routine fhv_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.
*
*  FHV_ESING
*     The specified matrix is singular within the working precision.
*
*  FHV_ECOND
*     The specified matrix is ill-conditioned.
*
*  For more details see comments to the routine luf_factorize.
*
*  ALGORITHM
*
*  The routine fhv_factorize calls the routine luf_factorize (see the
*  module GLPLUF), which actually computes LU-factorization of the basis
*  matrix B in the form
*
*     [B] = (F, V, P, Q),
*
*  where F and V are such matrices that
*
*     B = F * V,
*
*  and P and Q are such permutation matrices that the matrix
*
*     L = P * F * inv(P)
*
*  is lower triangular with unity diagonal, and the matrix
*
*     U = P * V * Q
*
*  is upper triangular.
*
*  In order to build the complete representation of the factorization
*  (see formula (1) in the file glpfhv.h) the routine fhv_factorize just
*  additionally sets H = I and P0 = P. */

int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
      int ind[], double val[]), void *info)
{     int ret;
      if (m < 1)
         xfault("fhv_factorize: m = %d; invalid parameter\n", m);
      if (m > M_MAX)
         xfault("fhv_factorize: m = %d; matrix too big\n", m);
      fhv->m = m;
      /* invalidate the factorization */
      fhv->valid = 0;
      /* allocate/reallocate arrays, if necessary */
      if (fhv->hh_ind == NULL)
         fhv->hh_ind = xcalloc(1+fhv->hh_max, sizeof(int));
      if (fhv->hh_ptr == NULL)
         fhv->hh_ptr = xcalloc(1+fhv->hh_max, sizeof(int));
      if (fhv->hh_len == NULL)
         fhv->hh_len = xcalloc(1+fhv->hh_max, sizeof(int));
      if (fhv->m_max < m)
      {  if (fhv->p0_row != NULL) xfree(fhv->p0_row);
         if (fhv->p0_col != NULL) xfree(fhv->p0_col);
         if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);
         if (fhv->cc_val != NULL) xfree(fhv->cc_val);
         fhv->m_max = m + 100;
         fhv->p0_row = xcalloc(1+fhv->m_max, sizeof(int));
         fhv->p0_col = xcalloc(1+fhv->m_max, sizeof(int));
         fhv->cc_ind = xcalloc(1+fhv->m_max, sizeof(int));
         fhv->cc_val = xcalloc(1+fhv->m_max, sizeof(double));
      }
      /* try to factorize the basis matrix */
      switch (luf_factorize(fhv->luf, m, col, info))
      {  case 0:
            break;
         case LUF_ESING:
            ret = FHV_ESING;
            goto done;
         case LUF_ECOND:
            ret = FHV_ECOND;
            goto done;
         default:
            xassert(fhv != fhv);
      }
      /* the basis matrix has been successfully factorized */
      fhv->valid = 1;
      /* H := I */
      fhv->hh_nfs = 0;
      /* P0 := P */
      memcpy(&fhv->p0_row[1], &fhv->luf->pp_row[1], sizeof(int) * m);
      memcpy(&fhv->p0_col[1], &fhv->luf->pp_col[1], sizeof(int) * m);
      /* currently H has no factors */
      fhv->nnz_h = 0;
      ret = 0;
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  NAME
*
*  fhv_h_solve - solve system H*x = b or H'*x = b
*
*  SYNOPSIS
*
*  #include "glpfhv.h"
*  void fhv_h_solve(FHV *fhv, int tr, double x[]);
*
*  DESCRIPTION
*
*  The routine fhv_h_solve solves either the system H*x = b (if the
*  flag tr is zero) or the system H'*x = b (if the flag tr is non-zero),
*  where the matrix H is a component of the factorization specified by
*  the parameter fhv, H' is a matrix transposed to H.
*
*  On entry the array x should contain elements of the right-hand side
*  vector b in locations x[1], ..., x[m], where m is the order of the
*  matrix H. On exit this array will contain elements of the solution
*  vector x in the same locations. */

void fhv_h_solve(FHV *fhv, int tr, double x[])
{     int nfs = fhv->hh_nfs;
      int *hh_ind = fhv->hh_ind;
      int *hh_ptr = fhv->hh_ptr;
      int *hh_len = fhv->hh_len;
      int *sv_ind = fhv->luf->sv_ind;
      double *sv_val = fhv->luf->sv_val;
      int i, k, beg, end, ptr;
      double temp;
      if (!fhv->valid)
         xfault("fhv_h_solve: the factorization is not valid\n");
      if (!tr)
      {  /* solve the system H*x = b */
         for (k = 1; k <= nfs; k++)
         {  i = hh_ind[k];
            temp = x[i];
            beg = hh_ptr[k];
            end = beg + hh_len[k] - 1;
            for (ptr = beg; ptr <= end; ptr++)
               temp -= sv_val[ptr] * x[sv_ind[ptr]];
            x[i] = temp;
         }
      }
      else
      {  /* solve the system H'*x = b */
         for (k = nfs; k >= 1; k--)
         {  i = hh_ind[k];
            temp = x[i];
            if (temp == 0.0) continue;
            beg = hh_ptr[k];
            end = beg + hh_len[k] - 1;
            for (ptr = beg; ptr <= end; ptr++)
               x[sv_ind[ptr]] -= sv_val[ptr] * temp;
         }
      }
      return;
}

/***********************************************************************
*  NAME
*
*  fhv_ftran - perform forward transformation (solve system B*x = b)
*
*  SYNOPSIS
*
*  #include "glpfhv.h"
*  void fhv_ftran(FHV *fhv, double x[]);
*
*  DESCRIPTION
*
*  The routine fhv_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 fhv_ftran(FHV *fhv, double x[])
{     int *pp_row = fhv->luf->pp_row;
      int *pp_col = fhv->luf->pp_col;
      int *p0_row = fhv->p0_row;
      int *p0_col = fhv->p0_col;
      if (!fhv->valid)
         xfault("fhv_ftran: the factorization is not valid\n");
      /* B = F*H*V, therefore inv(B) = inv(V)*inv(H)*inv(F) */
      fhv->luf->pp_row = p0_row;
      fhv->luf->pp_col = p0_col;
      luf_f_solve(fhv->luf, 0, x);
      fhv->luf->pp_row = pp_row;
      fhv->luf->pp_col = pp_col;
      fhv_h_solve(fhv, 0, x);
      luf_v_solve(fhv->luf, 0, x);
      return;
}

/***********************************************************************
*  NAME
*
*  fhv_btran - perform backward transformation (solve system B'*x = b)
*
*  SYNOPSIS
*
*  #include "glpfhv.h"
*  void fhv_btran(FHV *fhv, double x[]);
*
*  DESCRIPTION
*
*  The routine fhv_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 fhv_btran(FHV *fhv, double x[])
{     int *pp_row = fhv->luf->pp_row;
      int *pp_col = fhv->luf->pp_col;
      int *p0_row = fhv->p0_row;
      int *p0_col = fhv->p0_col;
      if (!fhv->valid)
         xfault("fhv_btran: the factorization is not valid\n");
      /* B = F*H*V, therefore inv(B') = inv(F')*inv(H')*inv(V') */
      luf_v_solve(fhv->luf, 1, x);
      fhv_h_solve(fhv, 1, x);
      fhv->luf->pp_row = p0_row;
      fhv->luf->pp_col = p0_col;
      luf_f_solve(fhv->luf, 1, x);
      fhv->luf->pp_row = pp_row;
      fhv->luf->pp_col = pp_col;
      return;
}

/***********************************************************************
*  NAME
*
*  fhv_update_it - update LP basis factorization
*
*  SYNOPSIS
*
*  #include "glpfhv.h"
*  int fhv_update_it(FHV *fhv, int j, int len, const int ind[],
*     const double val[]);
*
*  DESCRIPTION
*
*  The routine fhv_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.
*
*  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.
*
*  FHV_ESING
*     The adjacent basis matrix is structurally singular, since after
*     changing j-th column of matrix V by the new column (see algorithm
*     below) the case k1 > k2 occured.
*
*  FHV_ECHECK
*     The factorization is inaccurate, since after transforming k2-th
*     row of matrix U = P*V*Q, its diagonal element u[k2,k2] is zero or
*     close to zero,
*
*  FHV_ELIMIT
*     Maximal number of H factors has been reached.
*
*  FHV_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
*  fhv_factorize.
*
*  ALGORITHM
*
*  The routine fhv_update_it is based on the transformation proposed by
*  Forrest and Tomlin.
*
*  Let j-th column of the basis matrix B have been replaced by new
*  column B[j]. In order to keep the equality B = F*H*V j-th column of
*  matrix V should be replaced by the column inv(F*H)*B[j].
*
*  From the standpoint of matrix U = P*V*Q, replacement of j-th column