1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpfhv.c Mon Dec 06 13:09:21 2010 +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 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 */