1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpfhv.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,774 @@
     1.4 +/* glpfhv.c (LP basis factorization, FHV eta file version) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpfhv.h"
    1.29 +#include "glpenv.h"
    1.30 +#define xfault xerror
    1.31 +
    1.32 +/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
    1.33 +
    1.34 +#define M_MAX 100000000 /* = 100*10^6 */
    1.35 +/* maximal order of the basis matrix */
    1.36 +
    1.37 +/***********************************************************************
    1.38 +*  NAME
    1.39 +*
    1.40 +*  fhv_create_it - create LP basis factorization
    1.41 +*
    1.42 +*  SYNOPSIS
    1.43 +*
    1.44 +*  #include "glpfhv.h"
    1.45 +*  FHV *fhv_create_it(void);
    1.46 +*
    1.47 +*  DESCRIPTION
    1.48 +*
    1.49 +*  The routine fhv_create_it creates a program object, which represents
    1.50 +*  a factorization of LP basis.
    1.51 +*
    1.52 +*  RETURNS
    1.53 +*
    1.54 +*  The routine fhv_create_it returns a pointer to the object created. */
    1.55 +
    1.56 +FHV *fhv_create_it(void)
    1.57 +{     FHV *fhv;
    1.58 +      fhv = xmalloc(sizeof(FHV));
    1.59 +      fhv->m_max = fhv->m = 0;
    1.60 +      fhv->valid = 0;
    1.61 +      fhv->luf = luf_create_it();
    1.62 +      fhv->hh_max = 50;
    1.63 +      fhv->hh_nfs = 0;
    1.64 +      fhv->hh_ind = fhv->hh_ptr = fhv->hh_len = NULL;
    1.65 +      fhv->p0_row = fhv->p0_col = NULL;
    1.66 +      fhv->cc_ind = NULL;
    1.67 +      fhv->cc_val = NULL;
    1.68 +      fhv->upd_tol = 1e-6;
    1.69 +      fhv->nnz_h = 0;
    1.70 +      return fhv;
    1.71 +}
    1.72 +
    1.73 +/***********************************************************************
    1.74 +*  NAME
    1.75 +*
    1.76 +*  fhv_factorize - compute LP basis factorization
    1.77 +*
    1.78 +*  SYNOPSIS
    1.79 +*
    1.80 +*  #include "glpfhv.h"
    1.81 +*  int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
    1.82 +*     int ind[], double val[]), void *info);
    1.83 +*
    1.84 +*  DESCRIPTION
    1.85 +*
    1.86 +*  The routine fhv_factorize computes the factorization of the basis
    1.87 +*  matrix B specified by the routine col.
    1.88 +*
    1.89 +*  The parameter fhv specified the basis factorization data structure
    1.90 +*  created by the routine fhv_create_it.
    1.91 +*
    1.92 +*  The parameter m specifies the order of B, m > 0.
    1.93 +*
    1.94 +*  The formal routine col specifies the matrix B to be factorized. To
    1.95 +*  obtain j-th column of A the routine fhv_factorize calls the routine
    1.96 +*  col with the parameter j (1 <= j <= n). In response the routine col
    1.97 +*  should store row indices and numerical values of non-zero elements
    1.98 +*  of j-th column of B to locations ind[1,...,len] and val[1,...,len],
    1.99 +*  respectively, where len is the number of non-zeros in j-th column
   1.100 +*  returned on exit. Neither zero nor duplicate elements are allowed.
   1.101 +*
   1.102 +*  The parameter info is a transit pointer passed to the routine col.
   1.103 +*
   1.104 +*  RETURNS
   1.105 +*
   1.106 +*  0  The factorization has been successfully computed.
   1.107 +*
   1.108 +*  FHV_ESING
   1.109 +*     The specified matrix is singular within the working precision.
   1.110 +*
   1.111 +*  FHV_ECOND
   1.112 +*     The specified matrix is ill-conditioned.
   1.113 +*
   1.114 +*  For more details see comments to the routine luf_factorize.
   1.115 +*
   1.116 +*  ALGORITHM
   1.117 +*
   1.118 +*  The routine fhv_factorize calls the routine luf_factorize (see the
   1.119 +*  module GLPLUF), which actually computes LU-factorization of the basis
   1.120 +*  matrix B in the form
   1.121 +*
   1.122 +*     [B] = (F, V, P, Q),
   1.123 +*
   1.124 +*  where F and V are such matrices that
   1.125 +*
   1.126 +*     B = F * V,
   1.127 +*
   1.128 +*  and P and Q are such permutation matrices that the matrix
   1.129 +*
   1.130 +*     L = P * F * inv(P)
   1.131 +*
   1.132 +*  is lower triangular with unity diagonal, and the matrix
   1.133 +*
   1.134 +*     U = P * V * Q
   1.135 +*
   1.136 +*  is upper triangular.
   1.137 +*
   1.138 +*  In order to build the complete representation of the factorization
   1.139 +*  (see formula (1) in the file glpfhv.h) the routine fhv_factorize just
   1.140 +*  additionally sets H = I and P0 = P. */
   1.141 +
   1.142 +int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
   1.143 +      int ind[], double val[]), void *info)
   1.144 +{     int ret;
   1.145 +      if (m < 1)
   1.146 +         xfault("fhv_factorize: m = %d; invalid parameter\n", m);
   1.147 +      if (m > M_MAX)
   1.148 +         xfault("fhv_factorize: m = %d; matrix too big\n", m);
   1.149 +      fhv->m = m;
   1.150 +      /* invalidate the factorization */
   1.151 +      fhv->valid = 0;
   1.152 +      /* allocate/reallocate arrays, if necessary */
   1.153 +      if (fhv->hh_ind == NULL)
   1.154 +         fhv->hh_ind = xcalloc(1+fhv->hh_max, sizeof(int));
   1.155 +      if (fhv->hh_ptr == NULL)
   1.156 +         fhv->hh_ptr = xcalloc(1+fhv->hh_max, sizeof(int));
   1.157 +      if (fhv->hh_len == NULL)
   1.158 +         fhv->hh_len = xcalloc(1+fhv->hh_max, sizeof(int));
   1.159 +      if (fhv->m_max < m)
   1.160 +      {  if (fhv->p0_row != NULL) xfree(fhv->p0_row);
   1.161 +         if (fhv->p0_col != NULL) xfree(fhv->p0_col);
   1.162 +         if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);
   1.163 +         if (fhv->cc_val != NULL) xfree(fhv->cc_val);
   1.164 +         fhv->m_max = m + 100;
   1.165 +         fhv->p0_row = xcalloc(1+fhv->m_max, sizeof(int));
   1.166 +         fhv->p0_col = xcalloc(1+fhv->m_max, sizeof(int));
   1.167 +         fhv->cc_ind = xcalloc(1+fhv->m_max, sizeof(int));
   1.168 +         fhv->cc_val = xcalloc(1+fhv->m_max, sizeof(double));
   1.169 +      }
   1.170 +      /* try to factorize the basis matrix */
   1.171 +      switch (luf_factorize(fhv->luf, m, col, info))
   1.172 +      {  case 0:
   1.173 +            break;
   1.174 +         case LUF_ESING:
   1.175 +            ret = FHV_ESING;
   1.176 +            goto done;
   1.177 +         case LUF_ECOND:
   1.178 +            ret = FHV_ECOND;
   1.179 +            goto done;
   1.180 +         default:
   1.181 +            xassert(fhv != fhv);
   1.182 +      }
   1.183 +      /* the basis matrix has been successfully factorized */
   1.184 +      fhv->valid = 1;
   1.185 +      /* H := I */
   1.186 +      fhv->hh_nfs = 0;
   1.187 +      /* P0 := P */
   1.188 +      memcpy(&fhv->p0_row[1], &fhv->luf->pp_row[1], sizeof(int) * m);
   1.189 +      memcpy(&fhv->p0_col[1], &fhv->luf->pp_col[1], sizeof(int) * m);
   1.190 +      /* currently H has no factors */
   1.191 +      fhv->nnz_h = 0;
   1.192 +      ret = 0;
   1.193 +done: /* return to the calling program */
   1.194 +      return ret;
   1.195 +}
   1.196 +
   1.197 +/***********************************************************************
   1.198 +*  NAME
   1.199 +*
   1.200 +*  fhv_h_solve - solve system H*x = b or H'*x = b
   1.201 +*
   1.202 +*  SYNOPSIS
   1.203 +*
   1.204 +*  #include "glpfhv.h"
   1.205 +*  void fhv_h_solve(FHV *fhv, int tr, double x[]);
   1.206 +*
   1.207 +*  DESCRIPTION
   1.208 +*
   1.209 +*  The routine fhv_h_solve solves either the system H*x = b (if the
   1.210 +*  flag tr is zero) or the system H'*x = b (if the flag tr is non-zero),
   1.211 +*  where the matrix H is a component of the factorization specified by
   1.212 +*  the parameter fhv, H' is a matrix transposed to H.
   1.213 +*
   1.214 +*  On entry the array x should contain elements of the right-hand side
   1.215 +*  vector b in locations x[1], ..., x[m], where m is the order of the
   1.216 +*  matrix H. On exit this array will contain elements of the solution
   1.217 +*  vector x in the same locations. */
   1.218 +
   1.219 +void fhv_h_solve(FHV *fhv, int tr, double x[])
   1.220 +{     int nfs = fhv->hh_nfs;
   1.221 +      int *hh_ind = fhv->hh_ind;
   1.222 +      int *hh_ptr = fhv->hh_ptr;
   1.223 +      int *hh_len = fhv->hh_len;
   1.224 +      int *sv_ind = fhv->luf->sv_ind;
   1.225 +      double *sv_val = fhv->luf->sv_val;
   1.226 +      int i, k, beg, end, ptr;
   1.227 +      double temp;
   1.228 +      if (!fhv->valid)
   1.229 +         xfault("fhv_h_solve: the factorization is not valid\n");
   1.230 +      if (!tr)
   1.231 +      {  /* solve the system H*x = b */
   1.232 +         for (k = 1; k <= nfs; k++)
   1.233 +         {  i = hh_ind[k];
   1.234 +            temp = x[i];
   1.235 +            beg = hh_ptr[k];
   1.236 +            end = beg + hh_len[k] - 1;
   1.237 +            for (ptr = beg; ptr <= end; ptr++)
   1.238 +               temp -= sv_val[ptr] * x[sv_ind[ptr]];
   1.239 +            x[i] = temp;
   1.240 +         }
   1.241 +      }
   1.242 +      else
   1.243 +      {  /* solve the system H'*x = b */
   1.244 +         for (k = nfs; k >= 1; k--)
   1.245 +         {  i = hh_ind[k];
   1.246 +            temp = x[i];
   1.247 +            if (temp == 0.0) continue;
   1.248 +            beg = hh_ptr[k];
   1.249 +            end = beg + hh_len[k] - 1;
   1.250 +            for (ptr = beg; ptr <= end; ptr++)
   1.251 +               x[sv_ind[ptr]] -= sv_val[ptr] * temp;
   1.252 +         }
   1.253 +      }
   1.254 +      return;
   1.255 +}
   1.256 +
   1.257 +/***********************************************************************
   1.258 +*  NAME
   1.259 +*
   1.260 +*  fhv_ftran - perform forward transformation (solve system B*x = b)
   1.261 +*
   1.262 +*  SYNOPSIS
   1.263 +*
   1.264 +*  #include "glpfhv.h"
   1.265 +*  void fhv_ftran(FHV *fhv, double x[]);
   1.266 +*
   1.267 +*  DESCRIPTION
   1.268 +*
   1.269 +*  The routine fhv_ftran performs forward transformation, i.e. solves
   1.270 +*  the system B*x = b, where B is the basis matrix, x is the vector of
   1.271 +*  unknowns to be computed, b is the vector of right-hand sides.
   1.272 +*
   1.273 +*  On entry elements of the vector b should be stored in dense format
   1.274 +*  in locations x[1], ..., x[m], where m is the number of rows. On exit
   1.275 +*  the routine stores elements of the vector x in the same locations. */
   1.276 +
   1.277 +void fhv_ftran(FHV *fhv, double x[])
   1.278 +{     int *pp_row = fhv->luf->pp_row;
   1.279 +      int *pp_col = fhv->luf->pp_col;
   1.280 +      int *p0_row = fhv->p0_row;
   1.281 +      int *p0_col = fhv->p0_col;
   1.282 +      if (!fhv->valid)
   1.283 +         xfault("fhv_ftran: the factorization is not valid\n");
   1.284 +      /* B = F*H*V, therefore inv(B) = inv(V)*inv(H)*inv(F) */
   1.285 +      fhv->luf->pp_row = p0_row;
   1.286 +      fhv->luf->pp_col = p0_col;
   1.287 +      luf_f_solve(fhv->luf, 0, x);
   1.288 +      fhv->luf->pp_row = pp_row;
   1.289 +      fhv->luf->pp_col = pp_col;
   1.290 +      fhv_h_solve(fhv, 0, x);
   1.291 +      luf_v_solve(fhv->luf, 0, x);
   1.292 +      return;
   1.293 +}
   1.294 +
   1.295 +/***********************************************************************
   1.296 +*  NAME
   1.297 +*
   1.298 +*  fhv_btran - perform backward transformation (solve system B'*x = b)
   1.299 +*
   1.300 +*  SYNOPSIS
   1.301 +*
   1.302 +*  #include "glpfhv.h"
   1.303 +*  void fhv_btran(FHV *fhv, double x[]);
   1.304 +*
   1.305 +*  DESCRIPTION
   1.306 +*
   1.307 +*  The routine fhv_btran performs backward transformation, i.e. solves
   1.308 +*  the system B'*x = b, where B' is a matrix transposed to the basis
   1.309 +*  matrix B, x is the vector of unknowns to be computed, b is the vector
   1.310 +*  of right-hand sides.
   1.311 +*
   1.312 +*  On entry elements of the vector b should be stored in dense format
   1.313 +*  in locations x[1], ..., x[m], where m is the number of rows. On exit
   1.314 +*  the routine stores elements of the vector x in the same locations. */
   1.315 +
   1.316 +void fhv_btran(FHV *fhv, double x[])
   1.317 +{     int *pp_row = fhv->luf->pp_row;
   1.318 +      int *pp_col = fhv->luf->pp_col;
   1.319 +      int *p0_row = fhv->p0_row;
   1.320 +      int *p0_col = fhv->p0_col;
   1.321 +      if (!fhv->valid)
   1.322 +         xfault("fhv_btran: the factorization is not valid\n");
   1.323 +      /* B = F*H*V, therefore inv(B') = inv(F')*inv(H')*inv(V') */
   1.324 +      luf_v_solve(fhv->luf, 1, x);
   1.325 +      fhv_h_solve(fhv, 1, x);
   1.326 +      fhv->luf->pp_row = p0_row;
   1.327 +      fhv->luf->pp_col = p0_col;
   1.328 +      luf_f_solve(fhv->luf, 1, x);
   1.329 +      fhv->luf->pp_row = pp_row;
   1.330 +      fhv->luf->pp_col = pp_col;
   1.331 +      return;
   1.332 +}
   1.333 +
   1.334 +/***********************************************************************
   1.335 +*  NAME
   1.336 +*
   1.337 +*  fhv_update_it - update LP basis factorization
   1.338 +*
   1.339 +*  SYNOPSIS
   1.340 +*
   1.341 +*  #include "glpfhv.h"
   1.342 +*  int fhv_update_it(FHV *fhv, int j, int len, const int ind[],
   1.343 +*     const double val[]);
   1.344 +*
   1.345 +*  DESCRIPTION
   1.346 +*
   1.347 +*  The routine fhv_update_it updates the factorization of the basis
   1.348 +*  matrix B after replacing its j-th column by a new vector.
   1.349 +*
   1.350 +*  The parameter j specifies the number of column of B, which has been
   1.351 +*  replaced, 1 <= j <= m, where m is the order of B.
   1.352 +*
   1.353 +*  Row indices and numerical values of non-zero elements of the new
   1.354 +*  column of B should be placed in locations ind[1], ..., ind[len] and
   1.355 +*  val[1], ..., val[len], resp., where len is the number of non-zeros
   1.356 +*  in the column. Neither zero nor duplicate elements are allowed.
   1.357 +*
   1.358 +*  RETURNS
   1.359 +*
   1.360 +*  0  The factorization has been successfully updated.
   1.361 +*
   1.362 +*  FHV_ESING
   1.363 +*     The adjacent basis matrix is structurally singular, since after
   1.364 +*     changing j-th column of matrix V by the new column (see algorithm
   1.365 +*     below) the case k1 > k2 occured.
   1.366 +*
   1.367 +*  FHV_ECHECK
   1.368 +*     The factorization is inaccurate, since after transforming k2-th
   1.369 +*     row of matrix U = P*V*Q, its diagonal element u[k2,k2] is zero or
   1.370 +*     close to zero,
   1.371 +*
   1.372 +*  FHV_ELIMIT
   1.373 +*     Maximal number of H factors has been reached.
   1.374 +*
   1.375 +*  FHV_EROOM
   1.376 +*     Overflow of the sparse vector area.
   1.377 +*
   1.378 +*  In case of non-zero return code the factorization becomes invalid.
   1.379 +*  It should not be used until it has been recomputed with the routine
   1.380 +*  fhv_factorize.
   1.381 +*
   1.382 +*  ALGORITHM
   1.383 +*
   1.384 +*  The routine fhv_update_it is based on the transformation proposed by
   1.385 +*  Forrest and Tomlin.
   1.386 +*
   1.387 +*  Let j-th column of the basis matrix B have been replaced by new
   1.388 +*  column B[j]. In order to keep the equality B = F*H*V j-th column of
   1.389 +*  matrix V should be replaced by the column inv(F*H)*B[j].
   1.390 +*
   1.391 +*  From the standpoint of matrix U = P*V*Q, replacement of j-th column
   1.392 +*  of matrix V is equivalent to replacement of k1-th column of matrix U,
   1.393 +*  where k1 is determined by permutation matrix Q. Thus, matrix U loses
   1.394 +*  its upper triangular form and becomes the following:
   1.395 +*
   1.396 +*         1   k1       k2   m
   1.397 +*     1   x x * x x x x x x x
   1.398 +*         . x * x x x x x x x
   1.399 +*     k1  . . * x x x x x x x
   1.400 +*         . . * x x x x x x x
   1.401 +*         . . * . x x x x x x
   1.402 +*         . . * . . x x x x x
   1.403 +*         . . * . . . x x x x
   1.404 +*     k2  . . * . . . . x x x
   1.405 +*         . . . . . . . . x x
   1.406 +*     m   . . . . . . . . . x
   1.407 +*
   1.408 +*  where row index k2 corresponds to the lowest non-zero element of
   1.409 +*  k1-th column.
   1.410 +*
   1.411 +*  The routine moves rows and columns k1+1, k1+2, ..., k2 of matrix U
   1.412 +*  by one position to the left and upwards and moves k1-th row and k1-th
   1.413 +*  column to position k2. As the result of such symmetric permutations
   1.414 +*  matrix U becomes the following:
   1.415 +*
   1.416 +*         1   k1       k2   m
   1.417 +*     1   x x x x x x x * x x
   1.418 +*         . x x x x x x * x x
   1.419 +*     k1  . . x x x x x * x x
   1.420 +*         . . . x x x x * x x
   1.421 +*         . . . . x x x * x x
   1.422 +*         . . . . . x x * x x
   1.423 +*         . . . . . . x * x x
   1.424 +*     k2  . . x x x x x * x x
   1.425 +*         . . . . . . . . x x
   1.426 +*     m   . . . . . . . . . x
   1.427 +*
   1.428 +*  Then the routine performs gaussian elimination to eliminate elements
   1.429 +*  u[k2,k1], u[k2,k1+1], ..., u[k2,k2-1] using diagonal elements
   1.430 +*  u[k1,k1], u[k1+1,k1+1], ..., u[k2-1,k2-1] as pivots in the same way
   1.431 +*  as described in comments to the routine luf_factorize (see the module
   1.432 +*  GLPLUF). Note that actually all operations are performed on matrix V,
   1.433 +*  not on matrix U. During the elimination process the routine permutes
   1.434 +*  neither rows nor columns, so only k2-th row of matrix U is changed.
   1.435 +*
   1.436 +*  To keep the main equality B = F*H*V, each time when the routine
   1.437 +*  applies elementary gaussian transformation to the transformed row of
   1.438 +*  matrix V (which corresponds to k2-th row of matrix U), it also adds
   1.439 +*  a new element (gaussian multiplier) to the current row-like factor
   1.440 +*  of matrix H, which corresponds to the transformed row of matrix V. */
   1.441 +
   1.442 +int fhv_update_it(FHV *fhv, int j, int len, const int ind[],
   1.443 +      const double val[])
   1.444 +{     int m = fhv->m;
   1.445 +      LUF *luf = fhv->luf;
   1.446 +      int *vr_ptr = luf->vr_ptr;
   1.447 +      int *vr_len = luf->vr_len;
   1.448 +      int *vr_cap = luf->vr_cap;
   1.449 +      double *vr_piv = luf->vr_piv;
   1.450 +      int *vc_ptr = luf->vc_ptr;
   1.451 +      int *vc_len = luf->vc_len;
   1.452 +      int *vc_cap = luf->vc_cap;
   1.453 +      int *pp_row = luf->pp_row;
   1.454 +      int *pp_col = luf->pp_col;
   1.455 +      int *qq_row = luf->qq_row;
   1.456 +      int *qq_col = luf->qq_col;
   1.457 +      int *sv_ind = luf->sv_ind;
   1.458 +      double *sv_val = luf->sv_val;
   1.459 +      double *work = luf->work;
   1.460 +      double eps_tol = luf->eps_tol;
   1.461 +      int *hh_ind = fhv->hh_ind;
   1.462 +      int *hh_ptr = fhv->hh_ptr;
   1.463 +      int *hh_len = fhv->hh_len;
   1.464 +      int *p0_row = fhv->p0_row;
   1.465 +      int *p0_col = fhv->p0_col;
   1.466 +      int *cc_ind = fhv->cc_ind;
   1.467 +      double *cc_val = fhv->cc_val;
   1.468 +      double upd_tol = fhv->upd_tol;
   1.469 +      int i, i_beg, i_end, i_ptr, j_beg, j_end, j_ptr, k, k1, k2, p, q,
   1.470 +         p_beg, p_end, p_ptr, ptr, ret;
   1.471 +      double f, temp;
   1.472 +      if (!fhv->valid)
   1.473 +         xfault("fhv_update_it: the factorization is not valid\n");
   1.474 +      if (!(1 <= j && j <= m))
   1.475 +         xfault("fhv_update_it: j = %d; column number out of range\n",
   1.476 +            j);
   1.477 +      /* check if the new factor of matrix H can be created */
   1.478 +      if (fhv->hh_nfs == fhv->hh_max)
   1.479 +      {  /* maximal number of updates has been reached */
   1.480 +         fhv->valid = 0;
   1.481 +         ret = FHV_ELIMIT;
   1.482 +         goto done;
   1.483 +      }
   1.484 +      /* convert new j-th column of B to dense format */
   1.485 +      for (i = 1; i <= m; i++)
   1.486 +         cc_val[i] = 0.0;
   1.487 +      for (k = 1; k <= len; k++)
   1.488 +      {  i = ind[k];
   1.489 +         if (!(1 <= i && i <= m))
   1.490 +            xfault("fhv_update_it: ind[%d] = %d; row number out of rang"
   1.491 +               "e\n", k, i);
   1.492 +         if (cc_val[i] != 0.0)
   1.493 +            xfault("fhv_update_it: ind[%d] = %d; duplicate row index no"
   1.494 +               "t allowed\n", k, i);
   1.495 +         if (val[k] == 0.0)
   1.496 +            xfault("fhv_update_it: val[%d] = %g; zero element not allow"
   1.497 +               "ed\n", k, val[k]);
   1.498 +         cc_val[i] = val[k];
   1.499 +      }
   1.500 +      /* new j-th column of V := inv(F * H) * (new B[j]) */
   1.501 +      fhv->luf->pp_row = p0_row;
   1.502 +      fhv->luf->pp_col = p0_col;
   1.503 +      luf_f_solve(fhv->luf, 0, cc_val);
   1.504 +      fhv->luf->pp_row = pp_row;
   1.505 +      fhv->luf->pp_col = pp_col;
   1.506 +      fhv_h_solve(fhv, 0, cc_val);
   1.507 +      /* convert new j-th column of V to sparse format */
   1.508 +      len = 0;
   1.509 +      for (i = 1; i <= m; i++)
   1.510 +      {  temp = cc_val[i];
   1.511 +         if (temp == 0.0 || fabs(temp) < eps_tol) continue;
   1.512 +         len++, cc_ind[len] = i, cc_val[len] = temp;
   1.513 +      }
   1.514 +      /* clear old content of j-th column of matrix V */
   1.515 +      j_beg = vc_ptr[j];
   1.516 +      j_end = j_beg + vc_len[j] - 1;
   1.517 +      for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
   1.518 +      {  /* get row index of v[i,j] */
   1.519 +         i = sv_ind[j_ptr];
   1.520 +         /* find v[i,j] in the i-th row */
   1.521 +         i_beg = vr_ptr[i];
   1.522 +         i_end = i_beg + vr_len[i] - 1;
   1.523 +         for (i_ptr = i_beg; sv_ind[i_ptr] != j; i_ptr++) /* nop */;
   1.524 +         xassert(i_ptr <= i_end);
   1.525 +         /* remove v[i,j] from the i-th row */
   1.526 +         sv_ind[i_ptr] = sv_ind[i_end];
   1.527 +         sv_val[i_ptr] = sv_val[i_end];
   1.528 +         vr_len[i]--;
   1.529 +      }
   1.530 +      /* now j-th column of matrix V is empty */
   1.531 +      luf->nnz_v -= vc_len[j];
   1.532 +      vc_len[j] = 0;
   1.533 +      /* add new elements of j-th column of matrix V to corresponding
   1.534 +         row lists; determine indices k1 and k2 */
   1.535 +      k1 = qq_row[j], k2 = 0;
   1.536 +      for (ptr = 1; ptr <= len; ptr++)
   1.537 +      {  /* get row index of v[i,j] */
   1.538 +         i = cc_ind[ptr];
   1.539 +         /* at least one unused location is needed in i-th row */
   1.540 +         if (vr_len[i] + 1 > vr_cap[i])
   1.541 +         {  if (luf_enlarge_row(luf, i, vr_len[i] + 10))
   1.542 +            {  /* overflow of the sparse vector area */
   1.543 +               fhv->valid = 0;
   1.544 +               luf->new_sva = luf->sv_size + luf->sv_size;
   1.545 +               xassert(luf->new_sva > luf->sv_size);
   1.546 +               ret = FHV_EROOM;
   1.547 +               goto done;
   1.548 +            }
   1.549 +         }
   1.550 +         /* add v[i,j] to i-th row */
   1.551 +         i_ptr = vr_ptr[i] + vr_len[i];
   1.552 +         sv_ind[i_ptr] = j;
   1.553 +         sv_val[i_ptr] = cc_val[ptr];
   1.554 +         vr_len[i]++;
   1.555 +         /* adjust index k2 */
   1.556 +         if (k2 < pp_col[i]) k2 = pp_col[i];
   1.557 +      }
   1.558 +      /* capacity of j-th column (which is currently empty) should be
   1.559 +         not less than len locations */
   1.560 +      if (vc_cap[j] < len)
   1.561 +      {  if (luf_enlarge_col(luf, j, len))
   1.562 +         {  /* overflow of the sparse vector area */
   1.563 +            fhv->valid = 0;
   1.564 +            luf->new_sva = luf->sv_size + luf->sv_size;
   1.565 +            xassert(luf->new_sva > luf->sv_size);
   1.566 +            ret = FHV_EROOM;
   1.567 +            goto done;
   1.568 +         }
   1.569 +      }
   1.570 +      /* add new elements of matrix V to j-th column list */
   1.571 +      j_ptr = vc_ptr[j];
   1.572 +      memmove(&sv_ind[j_ptr], &cc_ind[1], len * sizeof(int));
   1.573 +      memmove(&sv_val[j_ptr], &cc_val[1], len * sizeof(double));
   1.574 +      vc_len[j] = len;
   1.575 +      luf->nnz_v += len;
   1.576 +      /* if k1 > k2, diagonal element u[k2,k2] of matrix U is zero and
   1.577 +         therefore the adjacent basis matrix is structurally singular */
   1.578 +      if (k1 > k2)
   1.579 +      {  fhv->valid = 0;
   1.580 +         ret = FHV_ESING;
   1.581 +         goto done;
   1.582 +      }
   1.583 +      /* perform implicit symmetric permutations of rows and columns of
   1.584 +         matrix U */
   1.585 +      i = pp_row[k1], j = qq_col[k1];
   1.586 +      for (k = k1; k < k2; k++)
   1.587 +      {  pp_row[k] = pp_row[k+1], pp_col[pp_row[k]] = k;
   1.588 +         qq_col[k] = qq_col[k+1], qq_row[qq_col[k]] = k;
   1.589 +      }
   1.590 +      pp_row[k2] = i, pp_col[i] = k2;
   1.591 +      qq_col[k2] = j, qq_row[j] = k2;
   1.592 +      /* now i-th row of the matrix V is k2-th row of matrix U; since
   1.593 +         no pivoting is used, only this row will be transformed */
   1.594 +      /* copy elements of i-th row of matrix V to the working array and
   1.595 +         remove these elements from matrix V */
   1.596 +      for (j = 1; j <= m; j++) work[j] = 0.0;
   1.597 +      i_beg = vr_ptr[i];
   1.598 +      i_end = i_beg + vr_len[i] - 1;
   1.599 +      for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
   1.600 +      {  /* get column index of v[i,j] */
   1.601 +         j = sv_ind[i_ptr];
   1.602 +         /* store v[i,j] to the working array */
   1.603 +         work[j] = sv_val[i_ptr];
   1.604 +         /* find v[i,j] in the j-th column */
   1.605 +         j_beg = vc_ptr[j];
   1.606 +         j_end = j_beg + vc_len[j] - 1;
   1.607 +         for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++) /* nop */;
   1.608 +         xassert(j_ptr <= j_end);
   1.609 +         /* remove v[i,j] from the j-th column */
   1.610 +         sv_ind[j_ptr] = sv_ind[j_end];
   1.611 +         sv_val[j_ptr] = sv_val[j_end];
   1.612 +         vc_len[j]--;
   1.613 +      }
   1.614 +      /* now i-th row of matrix V is empty */
   1.615 +      luf->nnz_v -= vr_len[i];
   1.616 +      vr_len[i] = 0;
   1.617 +      /* create the next row-like factor of the matrix H; this factor
   1.618 +         corresponds to i-th (transformed) row */
   1.619 +      fhv->hh_nfs++;
   1.620 +      hh_ind[fhv->hh_nfs] = i;
   1.621 +      /* hh_ptr[] will be set later */
   1.622 +      hh_len[fhv->hh_nfs] = 0;
   1.623 +      /* up to (k2 - k1) free locations are needed to add new elements
   1.624 +         to the non-trivial row of the row-like factor */
   1.625 +      if (luf->sv_end - luf->sv_beg < k2 - k1)
   1.626 +      {  luf_defrag_sva(luf);
   1.627 +         if (luf->sv_end - luf->sv_beg < k2 - k1)
   1.628 +         {  /* overflow of the sparse vector area */
   1.629 +            fhv->valid = luf->valid = 0;
   1.630 +            luf->new_sva = luf->sv_size + luf->sv_size;
   1.631 +            xassert(luf->new_sva > luf->sv_size);
   1.632 +            ret = FHV_EROOM;
   1.633 +            goto done;
   1.634 +         }
   1.635 +      }
   1.636 +      /* eliminate subdiagonal elements of matrix U */
   1.637 +      for (k = k1; k < k2; k++)
   1.638 +      {  /* v[p,q] = u[k,k] */
   1.639 +         p = pp_row[k], q = qq_col[k];
   1.640 +         /* this is the crucial point, where even tiny non-zeros should
   1.641 +            not be dropped */
   1.642 +         if (work[q] == 0.0) continue;
   1.643 +         /* compute gaussian multiplier f = v[i,q] / v[p,q] */
   1.644 +         f = work[q] / vr_piv[p];
   1.645 +         /* perform gaussian transformation:
   1.646 +            (i-th row) := (i-th row) - f * (p-th row)
   1.647 +            in order to eliminate v[i,q] = u[k2,k] */
   1.648 +         p_beg = vr_ptr[p];
   1.649 +         p_end = p_beg + vr_len[p] - 1;
   1.650 +         for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
   1.651 +            work[sv_ind[p_ptr]] -= f * sv_val[p_ptr];
   1.652 +         /* store new element (gaussian multiplier that corresponds to
   1.653 +            p-th row) in the current row-like factor */
   1.654 +         luf->sv_end--;
   1.655 +         sv_ind[luf->sv_end] = p;
   1.656 +         sv_val[luf->sv_end] = f;
   1.657 +         hh_len[fhv->hh_nfs]++;
   1.658 +      }
   1.659 +      /* set pointer to the current row-like factor of the matrix H
   1.660 +         (if no elements were added to this factor, it is unity matrix
   1.661 +         and therefore can be discarded) */
   1.662 +      if (hh_len[fhv->hh_nfs] == 0)
   1.663 +         fhv->hh_nfs--;
   1.664 +      else
   1.665 +      {  hh_ptr[fhv->hh_nfs] = luf->sv_end;
   1.666 +         fhv->nnz_h += hh_len[fhv->hh_nfs];
   1.667 +      }
   1.668 +      /* store new pivot which corresponds to u[k2,k2] */
   1.669 +      vr_piv[i] = work[qq_col[k2]];
   1.670 +      /* new elements of i-th row of matrix V (which are non-diagonal
   1.671 +         elements u[k2,k2+1], ..., u[k2,m] of matrix U = P*V*Q) now are
   1.672 +         contained in the working array; add them to matrix V */
   1.673 +      len = 0;
   1.674 +      for (k = k2+1; k <= m; k++)
   1.675 +      {  /* get column index and value of v[i,j] = u[k2,k] */
   1.676 +         j = qq_col[k];
   1.677 +         temp = work[j];
   1.678 +         /* if v[i,j] is close to zero, skip it */
   1.679 +         if (fabs(temp) < eps_tol) continue;
   1.680 +         /* at least one unused location is needed in j-th column */
   1.681 +         if (vc_len[j] + 1 > vc_cap[j])
   1.682 +         {  if (luf_enlarge_col(luf, j, vc_len[j] + 10))
   1.683 +            {  /* overflow of the sparse vector area */
   1.684 +               fhv->valid = 0;
   1.685 +               luf->new_sva = luf->sv_size + luf->sv_size;
   1.686 +               xassert(luf->new_sva > luf->sv_size);
   1.687 +               ret = FHV_EROOM;
   1.688 +               goto done;
   1.689 +            }
   1.690 +         }
   1.691 +         /* add v[i,j] to j-th column */
   1.692 +         j_ptr = vc_ptr[j] + vc_len[j];
   1.693 +         sv_ind[j_ptr] = i;
   1.694 +         sv_val[j_ptr] = temp;
   1.695 +         vc_len[j]++;
   1.696 +         /* also store v[i,j] to the auxiliary array */
   1.697 +         len++, cc_ind[len] = j, cc_val[len] = temp;
   1.698 +      }
   1.699 +      /* capacity of i-th row (which is currently empty) should be not
   1.700 +         less than len locations */
   1.701 +      if (vr_cap[i] < len)
   1.702 +      {  if (luf_enlarge_row(luf, i, len))
   1.703 +         {  /* overflow of the sparse vector area */
   1.704 +            fhv->valid = 0;
   1.705 +            luf->new_sva = luf->sv_size + luf->sv_size;
   1.706 +            xassert(luf->new_sva > luf->sv_size);
   1.707 +            ret = FHV_EROOM;
   1.708 +            goto done;
   1.709 +         }
   1.710 +      }
   1.711 +      /* add new elements to i-th row list */
   1.712 +      i_ptr = vr_ptr[i];
   1.713 +      memmove(&sv_ind[i_ptr], &cc_ind[1], len * sizeof(int));
   1.714 +      memmove(&sv_val[i_ptr], &cc_val[1], len * sizeof(double));
   1.715 +      vr_len[i] = len;
   1.716 +      luf->nnz_v += len;
   1.717 +      /* updating is finished; check that diagonal element u[k2,k2] is
   1.718 +         not very small in absolute value among other elements in k2-th
   1.719 +         row and k2-th column of matrix U = P*V*Q */
   1.720 +      /* temp = max(|u[k2,*]|, |u[*,k2]|) */
   1.721 +      temp = 0.0;
   1.722 +      /* walk through k2-th row of U which is i-th row of V */
   1.723 +      i = pp_row[k2];
   1.724 +      i_beg = vr_ptr[i];
   1.725 +      i_end = i_beg + vr_len[i] - 1;
   1.726 +      for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
   1.727 +         if (temp < fabs(sv_val[i_ptr])) temp = fabs(sv_val[i_ptr]);
   1.728 +      /* walk through k2-th column of U which is j-th column of V */
   1.729 +      j = qq_col[k2];
   1.730 +      j_beg = vc_ptr[j];
   1.731 +      j_end = j_beg + vc_len[j] - 1;
   1.732 +      for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
   1.733 +         if (temp < fabs(sv_val[j_ptr])) temp = fabs(sv_val[j_ptr]);
   1.734 +      /* check that u[k2,k2] is not very small */
   1.735 +      if (fabs(vr_piv[i]) < upd_tol * temp)
   1.736 +      {  /* the factorization seems to be inaccurate and therefore must
   1.737 +            be recomputed */
   1.738 +         fhv->valid = 0;
   1.739 +         ret = FHV_ECHECK;
   1.740 +         goto done;
   1.741 +      }
   1.742 +      /* the factorization has been successfully updated */
   1.743 +      ret = 0;
   1.744 +done: /* return to the calling program */
   1.745 +      return ret;
   1.746 +}
   1.747 +
   1.748 +/***********************************************************************
   1.749 +*  NAME
   1.750 +*
   1.751 +*  fhv_delete_it - delete LP basis factorization
   1.752 +*
   1.753 +*  SYNOPSIS
   1.754 +*
   1.755 +*  #include "glpfhv.h"
   1.756 +*  void fhv_delete_it(FHV *fhv);
   1.757 +*
   1.758 +*  DESCRIPTION
   1.759 +*
   1.760 +*  The routine fhv_delete_it deletes LP basis factorization specified
   1.761 +*  by the parameter fhv and frees all memory allocated to this program
   1.762 +*  object. */
   1.763 +
   1.764 +void fhv_delete_it(FHV *fhv)
   1.765 +{     luf_delete_it(fhv->luf);
   1.766 +      if (fhv->hh_ind != NULL) xfree(fhv->hh_ind);
   1.767 +      if (fhv->hh_ptr != NULL) xfree(fhv->hh_ptr);
   1.768 +      if (fhv->hh_len != NULL) xfree(fhv->hh_len);
   1.769 +      if (fhv->p0_row != NULL) xfree(fhv->p0_row);
   1.770 +      if (fhv->p0_col != NULL) xfree(fhv->p0_col);
   1.771 +      if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);
   1.772 +      if (fhv->cc_val != NULL) xfree(fhv->cc_val);
   1.773 +      xfree(fhv);
   1.774 +      return;
   1.775 +}
   1.776 +
   1.777 +/* eof */