lemon-project-template-glpk

annotate deps/glpk/src/glpfhv.c @ 11:4fc6ad2fb8a6

find changesets by author, revision, files, or words in the commit message
Test GLPK in src/main.cc
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 21:43:29 +0100
parents
children
rev   line source
alpar@9 1 /* glpfhv.c (LP basis factorization, FHV eta file version) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpfhv.h"
alpar@9 26 #include "glpenv.h"
alpar@9 27 #define xfault xerror
alpar@9 28
alpar@9 29 /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
alpar@9 30
alpar@9 31 #define M_MAX 100000000 /* = 100*10^6 */
alpar@9 32 /* maximal order of the basis matrix */
alpar@9 33
alpar@9 34 /***********************************************************************
alpar@9 35 * NAME
alpar@9 36 *
alpar@9 37 * fhv_create_it - create LP basis factorization
alpar@9 38 *
alpar@9 39 * SYNOPSIS
alpar@9 40 *
alpar@9 41 * #include "glpfhv.h"
alpar@9 42 * FHV *fhv_create_it(void);
alpar@9 43 *
alpar@9 44 * DESCRIPTION
alpar@9 45 *
alpar@9 46 * The routine fhv_create_it creates a program object, which represents
alpar@9 47 * a factorization of LP basis.
alpar@9 48 *
alpar@9 49 * RETURNS
alpar@9 50 *
alpar@9 51 * The routine fhv_create_it returns a pointer to the object created. */
alpar@9 52
alpar@9 53 FHV *fhv_create_it(void)
alpar@9 54 { FHV *fhv;
alpar@9 55 fhv = xmalloc(sizeof(FHV));
alpar@9 56 fhv->m_max = fhv->m = 0;
alpar@9 57 fhv->valid = 0;
alpar@9 58 fhv->luf = luf_create_it();
alpar@9 59 fhv->hh_max = 50;
alpar@9 60 fhv->hh_nfs = 0;
alpar@9 61 fhv->hh_ind = fhv->hh_ptr = fhv->hh_len = NULL;
alpar@9 62 fhv->p0_row = fhv->p0_col = NULL;
alpar@9 63 fhv->cc_ind = NULL;
alpar@9 64 fhv->cc_val = NULL;
alpar@9 65 fhv->upd_tol = 1e-6;
alpar@9 66 fhv->nnz_h = 0;
alpar@9 67 return fhv;
alpar@9 68 }
alpar@9 69
alpar@9 70 /***********************************************************************
alpar@9 71 * NAME
alpar@9 72 *
alpar@9 73 * fhv_factorize - compute LP basis factorization
alpar@9 74 *
alpar@9 75 * SYNOPSIS
alpar@9 76 *
alpar@9 77 * #include "glpfhv.h"
alpar@9 78 * int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
alpar@9 79 * int ind[], double val[]), void *info);
alpar@9 80 *
alpar@9 81 * DESCRIPTION
alpar@9 82 *
alpar@9 83 * The routine fhv_factorize computes the factorization of the basis
alpar@9 84 * matrix B specified by the routine col.
alpar@9 85 *
alpar@9 86 * The parameter fhv specified the basis factorization data structure
alpar@9 87 * created by the routine fhv_create_it.
alpar@9 88 *
alpar@9 89 * The parameter m specifies the order of B, m > 0.
alpar@9 90 *
alpar@9 91 * The formal routine col specifies the matrix B to be factorized. To
alpar@9 92 * obtain j-th column of A the routine fhv_factorize calls the routine
alpar@9 93 * col with the parameter j (1 <= j <= n). In response the routine col
alpar@9 94 * should store row indices and numerical values of non-zero elements
alpar@9 95 * of j-th column of B to locations ind[1,...,len] and val[1,...,len],
alpar@9 96 * respectively, where len is the number of non-zeros in j-th column
alpar@9 97 * returned on exit. Neither zero nor duplicate elements are allowed.
alpar@9 98 *
alpar@9 99 * The parameter info is a transit pointer passed to the routine col.
alpar@9 100 *
alpar@9 101 * RETURNS
alpar@9 102 *
alpar@9 103 * 0 The factorization has been successfully computed.
alpar@9 104 *
alpar@9 105 * FHV_ESING
alpar@9 106 * The specified matrix is singular within the working precision.
alpar@9 107 *
alpar@9 108 * FHV_ECOND
alpar@9 109 * The specified matrix is ill-conditioned.
alpar@9 110 *
alpar@9 111 * For more details see comments to the routine luf_factorize.
alpar@9 112 *
alpar@9 113 * ALGORITHM
alpar@9 114 *
alpar@9 115 * The routine fhv_factorize calls the routine luf_factorize (see the
alpar@9 116 * module GLPLUF), which actually computes LU-factorization of the basis
alpar@9 117 * matrix B in the form
alpar@9 118 *
alpar@9 119 * [B] = (F, V, P, Q),
alpar@9 120 *
alpar@9 121 * where F and V are such matrices that
alpar@9 122 *
alpar@9 123 * B = F * V,
alpar@9 124 *
alpar@9 125 * and P and Q are such permutation matrices that the matrix
alpar@9 126 *
alpar@9 127 * L = P * F * inv(P)
alpar@9 128 *
alpar@9 129 * is lower triangular with unity diagonal, and the matrix
alpar@9 130 *
alpar@9 131 * U = P * V * Q
alpar@9 132 *
alpar@9 133 * is upper triangular.
alpar@9 134 *
alpar@9 135 * In order to build the complete representation of the factorization
alpar@9 136 * (see formula (1) in the file glpfhv.h) the routine fhv_factorize just
alpar@9 137 * additionally sets H = I and P0 = P. */
alpar@9 138
alpar@9 139 int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j,
alpar@9 140 int ind[], double val[]), void *info)
alpar@9 141 { int ret;
alpar@9 142 if (m < 1)
alpar@9 143 xfault("fhv_factorize: m = %d; invalid parameter\n", m);
alpar@9 144 if (m > M_MAX)
alpar@9 145 xfault("fhv_factorize: m = %d; matrix too big\n", m);
alpar@9 146 fhv->m = m;
alpar@9 147 /* invalidate the factorization */
alpar@9 148 fhv->valid = 0;
alpar@9 149 /* allocate/reallocate arrays, if necessary */
alpar@9 150 if (fhv->hh_ind == NULL)
alpar@9 151 fhv->hh_ind = xcalloc(1+fhv->hh_max, sizeof(int));
alpar@9 152 if (fhv->hh_ptr == NULL)
alpar@9 153 fhv->hh_ptr = xcalloc(1+fhv->hh_max, sizeof(int));
alpar@9 154 if (fhv->hh_len == NULL)
alpar@9 155 fhv->hh_len = xcalloc(1+fhv->hh_max, sizeof(int));
alpar@9 156 if (fhv->m_max < m)
alpar@9 157 { if (fhv->p0_row != NULL) xfree(fhv->p0_row);
alpar@9 158 if (fhv->p0_col != NULL) xfree(fhv->p0_col);
alpar@9 159 if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);
alpar@9 160 if (fhv->cc_val != NULL) xfree(fhv->cc_val);
alpar@9 161 fhv->m_max = m + 100;
alpar@9 162 fhv->p0_row = xcalloc(1+fhv->m_max, sizeof(int));
alpar@9 163 fhv->p0_col = xcalloc(1+fhv->m_max, sizeof(int));
alpar@9 164 fhv->cc_ind = xcalloc(1+fhv->m_max, sizeof(int));
alpar@9 165 fhv->cc_val = xcalloc(1+fhv->m_max, sizeof(double));
alpar@9 166 }
alpar@9 167 /* try to factorize the basis matrix */
alpar@9 168 switch (luf_factorize(fhv->luf, m, col, info))
alpar@9 169 { case 0:
alpar@9 170 break;
alpar@9 171 case LUF_ESING:
alpar@9 172 ret = FHV_ESING;
alpar@9 173 goto done;
alpar@9 174 case LUF_ECOND:
alpar@9 175 ret = FHV_ECOND;
alpar@9 176 goto done;
alpar@9 177 default:
alpar@9 178 xassert(fhv != fhv);
alpar@9 179 }
alpar@9 180 /* the basis matrix has been successfully factorized */
alpar@9 181 fhv->valid = 1;
alpar@9 182 /* H := I */
alpar@9 183 fhv->hh_nfs = 0;
alpar@9 184 /* P0 := P */
alpar@9 185 memcpy(&fhv->p0_row[1], &fhv->luf->pp_row[1], sizeof(int) * m);
alpar@9 186 memcpy(&fhv->p0_col[1], &fhv->luf->pp_col[1], sizeof(int) * m);
alpar@9 187 /* currently H has no factors */
alpar@9 188 fhv->nnz_h = 0;
alpar@9 189 ret = 0;
alpar@9 190 done: /* return to the calling program */
alpar@9 191 return ret;
alpar@9 192 }
alpar@9 193
alpar@9 194 /***********************************************************************
alpar@9 195 * NAME
alpar@9 196 *
alpar@9 197 * fhv_h_solve - solve system H*x = b or H'*x = b
alpar@9 198 *
alpar@9 199 * SYNOPSIS
alpar@9 200 *
alpar@9 201 * #include "glpfhv.h"
alpar@9 202 * void fhv_h_solve(FHV *fhv, int tr, double x[]);
alpar@9 203 *
alpar@9 204 * DESCRIPTION
alpar@9 205 *
alpar@9 206 * The routine fhv_h_solve solves either the system H*x = b (if the
alpar@9 207 * flag tr is zero) or the system H'*x = b (if the flag tr is non-zero),
alpar@9 208 * where the matrix H is a component of the factorization specified by
alpar@9 209 * the parameter fhv, H' is a matrix transposed to H.
alpar@9 210 *
alpar@9 211 * On entry the array x should contain elements of the right-hand side
alpar@9 212 * vector b in locations x[1], ..., x[m], where m is the order of the
alpar@9 213 * matrix H. On exit this array will contain elements of the solution
alpar@9 214 * vector x in the same locations. */
alpar@9 215
alpar@9 216 void fhv_h_solve(FHV *fhv, int tr, double x[])
alpar@9 217 { int nfs = fhv->hh_nfs;
alpar@9 218 int *hh_ind = fhv->hh_ind;
alpar@9 219 int *hh_ptr = fhv->hh_ptr;
alpar@9 220 int *hh_len = fhv->hh_len;
alpar@9 221 int *sv_ind = fhv->luf->sv_ind;
alpar@9 222 double *sv_val = fhv->luf->sv_val;
alpar@9 223 int i, k, beg, end, ptr;
alpar@9 224 double temp;
alpar@9 225 if (!fhv->valid)
alpar@9 226 xfault("fhv_h_solve: the factorization is not valid\n");
alpar@9 227 if (!tr)
alpar@9 228 { /* solve the system H*x = b */
alpar@9 229 for (k = 1; k <= nfs; k++)
alpar@9 230 { i = hh_ind[k];
alpar@9 231 temp = x[i];
alpar@9 232 beg = hh_ptr[k];
alpar@9 233 end = beg + hh_len[k] - 1;
alpar@9 234 for (ptr = beg; ptr <= end; ptr++)
alpar@9 235 temp -= sv_val[ptr] * x[sv_ind[ptr]];
alpar@9 236 x[i] = temp;
alpar@9 237 }
alpar@9 238 }
alpar@9 239 else
alpar@9 240 { /* solve the system H'*x = b */
alpar@9 241 for (k = nfs; k >= 1; k--)
alpar@9 242 { i = hh_ind[k];
alpar@9 243 temp = x[i];
alpar@9 244 if (temp == 0.0) continue;
alpar@9 245 beg = hh_ptr[k];
alpar@9 246 end = beg + hh_len[k] - 1;
alpar@9 247 for (ptr = beg; ptr <= end; ptr++)
alpar@9 248 x[sv_ind[ptr]] -= sv_val[ptr] * temp;
alpar@9 249 }
alpar@9 250 }
alpar@9 251 return;
alpar@9 252 }
alpar@9 253
alpar@9 254 /***********************************************************************
alpar@9 255 * NAME
alpar@9 256 *
alpar@9 257 * fhv_ftran - perform forward transformation (solve system B*x = b)
alpar@9 258 *
alpar@9 259 * SYNOPSIS
alpar@9 260 *
alpar@9 261 * #include "glpfhv.h"
alpar@9 262 * void fhv_ftran(FHV *fhv, double x[]);
alpar@9 263 *
alpar@9 264 * DESCRIPTION
alpar@9 265 *
alpar@9 266 * The routine fhv_ftran performs forward transformation, i.e. solves
alpar@9 267 * the system B*x = b, where B is the basis matrix, x is the vector of
alpar@9 268 * unknowns to be computed, b is the vector of right-hand sides.
alpar@9 269 *
alpar@9 270 * On entry elements of the vector b should be stored in dense format
alpar@9 271 * in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@9 272 * the routine stores elements of the vector x in the same locations. */
alpar@9 273
alpar@9 274 void fhv_ftran(FHV *fhv, double x[])
alpar@9 275 { int *pp_row = fhv->luf->pp_row;
alpar@9 276 int *pp_col = fhv->luf->pp_col;
alpar@9 277 int *p0_row = fhv->p0_row;
alpar@9 278 int *p0_col = fhv->p0_col;
alpar@9 279 if (!fhv->valid)
alpar@9 280 xfault("fhv_ftran: the factorization is not valid\n");
alpar@9 281 /* B = F*H*V, therefore inv(B) = inv(V)*inv(H)*inv(F) */
alpar@9 282 fhv->luf->pp_row = p0_row;
alpar@9 283 fhv->luf->pp_col = p0_col;
alpar@9 284 luf_f_solve(fhv->luf, 0, x);
alpar@9 285 fhv->luf->pp_row = pp_row;
alpar@9 286 fhv->luf->pp_col = pp_col;
alpar@9 287 fhv_h_solve(fhv, 0, x);
alpar@9 288 luf_v_solve(fhv->luf, 0, x);
alpar@9 289 return;
alpar@9 290 }
alpar@9 291
alpar@9 292 /***********************************************************************
alpar@9 293 * NAME
alpar@9 294 *
alpar@9 295 * fhv_btran - perform backward transformation (solve system B'*x = b)
alpar@9 296 *
alpar@9 297 * SYNOPSIS
alpar@9 298 *
alpar@9 299 * #include "glpfhv.h"
alpar@9 300 * void fhv_btran(FHV *fhv, double x[]);
alpar@9 301 *
alpar@9 302 * DESCRIPTION
alpar@9 303 *
alpar@9 304 * The routine fhv_btran performs backward transformation, i.e. solves
alpar@9 305 * the system B'*x = b, where B' is a matrix transposed to the basis
alpar@9 306 * matrix B, x is the vector of unknowns to be computed, b is the vector
alpar@9 307 * of right-hand sides.
alpar@9 308 *
alpar@9 309 * On entry elements of the vector b should be stored in dense format
alpar@9 310 * in locations x[1], ..., x[m], where m is the number of rows. On exit
alpar@9 311 * the routine stores elements of the vector x in the same locations. */
alpar@9 312
alpar@9 313 void fhv_btran(FHV *fhv, double x[])
alpar@9 314 { int *pp_row = fhv->luf->pp_row;
alpar@9 315 int *pp_col = fhv->luf->pp_col;
alpar@9 316 int *p0_row = fhv->p0_row;
alpar@9 317 int *p0_col = fhv->p0_col;
alpar@9 318 if (!fhv->valid)
alpar@9 319 xfault("fhv_btran: the factorization is not valid\n");
alpar@9 320 /* B = F*H*V, therefore inv(B') = inv(F')*inv(H')*inv(V') */
alpar@9 321 luf_v_solve(fhv->luf, 1, x);
alpar@9 322 fhv_h_solve(fhv, 1, x);
alpar@9 323 fhv->luf->pp_row = p0_row;
alpar@9 324 fhv->luf->pp_col = p0_col;
alpar@9 325 luf_f_solve(fhv->luf, 1, x);
alpar@9 326 fhv->luf->pp_row = pp_row;
alpar@9 327 fhv->luf->pp_col = pp_col;
alpar@9 328 return;
alpar@9 329 }
alpar@9 330
alpar@9 331 /***********************************************************************
alpar@9 332 * NAME
alpar@9 333 *
alpar@9 334 * fhv_update_it - update LP basis factorization
alpar@9 335 *
alpar@9 336 * SYNOPSIS
alpar@9 337 *
alpar@9 338 * #include "glpfhv.h"
alpar@9 339 * int fhv_update_it(FHV *fhv, int j, int len, const int ind[],
alpar@9 340 * const double val[]);
alpar@9 341 *
alpar@9 342 * DESCRIPTION
alpar@9 343 *
alpar@9 344 * The routine fhv_update_it updates the factorization of the basis
alpar@9 345 * matrix B after replacing its j-th column by a new vector.
alpar@9 346 *
alpar@9 347 * The parameter j specifies the number of column of B, which has been
alpar@9 348 * replaced, 1 <= j <= m, where m is the order of B.
alpar@9 349 *
alpar@9 350 * Row indices and numerical values of non-zero elements of the new
alpar@9 351 * column of B should be placed in locations ind[1], ..., ind[len] and
alpar@9 352 * val[1], ..., val[len], resp., where len is the number of non-zeros
alpar@9 353 * in the column. Neither zero nor duplicate elements are allowed.
alpar@9 354 *
alpar@9 355 * RETURNS
alpar@9 356 *
alpar@9 357 * 0 The factorization has been successfully updated.
alpar@9 358 *
alpar@9 359 * FHV_ESING
alpar@9 360 * The adjacent basis matrix is structurally singular, since after
alpar@9 361 * changing j-th column of matrix V by the new column (see algorithm
alpar@9 362 * below) the case k1 > k2 occured.
alpar@9 363 *
alpar@9 364 * FHV_ECHECK
alpar@9 365 * The factorization is inaccurate, since after transforming k2-th
alpar@9 366 * row of matrix U = P*V*Q, its diagonal element u[k2,k2] is zero or
alpar@9 367 * close to zero,
alpar@9 368 *
alpar@9 369 * FHV_ELIMIT
alpar@9 370 * Maximal number of H factors has been reached.
alpar@9 371 *
alpar@9 372 * FHV_EROOM
alpar@9 373 * Overflow of the sparse vector area.
alpar@9 374 *
alpar@9 375 * In case of non-zero return code the factorization becomes invalid.
alpar@9 376 * It should not be used until it has been recomputed with the routine
alpar@9 377 * fhv_factorize.
alpar@9 378 *
alpar@9 379 * ALGORITHM
alpar@9 380 *
alpar@9 381 * The routine fhv_update_it is based on the transformation proposed by
alpar@9 382 * Forrest and Tomlin.
alpar@9 383 *
alpar@9 384 * Let j-th column of the basis matrix B have been replaced by new
alpar@9 385 * column B[j]. In order to keep the equality B = F*H*V j-th column of
alpar@9 386 * matrix V should be replaced by the column inv(F*H)*B[j].
alpar@9 387 *
alpar@9 388 * From the standpoint of matrix U = P*V*Q, replacement of j-th column
alpar@9 389 * of matrix V is equivalent to replacement of k1-th column of matrix U,
alpar@9 390 * where k1 is determined by permutation matrix Q. Thus, matrix U loses
alpar@9 391 * its upper triangular form and becomes the following:
alpar@9 392 *
alpar@9 393 * 1 k1 k2 m
alpar@9 394 * 1 x x * x x x x x x x
alpar@9 395 * . x * x x x x x x x
alpar@9 396 * k1 . . * x x x x x x x
alpar@9 397 * . . * x x x x x x x
alpar@9 398 * . . * . x x x x x x
alpar@9 399 * . . * . . x x x x x
alpar@9 400 * . . * . . . x x x x
alpar@9 401 * k2 . . * . . . . x x x
alpar@9 402 * . . . . . . . . x x
alpar@9 403 * m . . . . . . . . . x
alpar@9 404 *
alpar@9 405 * where row index k2 corresponds to the lowest non-zero element of
alpar@9 406 * k1-th column.
alpar@9 407 *
alpar@9 408 * The routine moves rows and columns k1+1, k1+2, ..., k2 of matrix U
alpar@9 409 * by one position to the left and upwards and moves k1-th row and k1-th
alpar@9 410 * column to position k2. As the result of such symmetric permutations
alpar@9 411 * matrix U becomes the following:
alpar@9 412 *
alpar@9 413 * 1 k1 k2 m
alpar@9 414 * 1 x x x x x x x * x x
alpar@9 415 * . x x x x x x * x x
alpar@9 416 * k1 . . x x x x x * x x
alpar@9 417 * . . . x x x x * x x
alpar@9 418 * . . . . x x x * x x
alpar@9 419 * . . . . . x x * x x
alpar@9 420 * . . . . . . x * x x
alpar@9 421 * k2 . . x x x x x * x x
alpar@9 422 * . . . . . . . . x x
alpar@9 423 * m . . . . . . . . . x
alpar@9 424 *
alpar@9 425 * Then the routine performs gaussian elimination to eliminate elements
alpar@9 426 * u[k2,k1], u[k2,k1+1], ..., u[k2,k2-1] using diagonal elements
alpar@9 427 * u[k1,k1], u[k1+1,k1+1], ..., u[k2-1,k2-1] as pivots in the same way
alpar@9 428 * as described in comments to the routine luf_factorize (see the module
alpar@9 429 * GLPLUF). Note that actually all operations are performed on matrix V,
alpar@9 430 * not on matrix U. During the elimination process the routine permutes
alpar@9 431 * neither rows nor columns, so only k2-th row of matrix U is changed.
alpar@9 432 *
alpar@9 433 * To keep the main equality B = F*H*V, each time when the routine
alpar@9 434 * applies elementary gaussian transformation to the transformed row of
alpar@9 435 * matrix V (which corresponds to k2-th row of matrix U), it also adds
alpar@9 436 * a new element (gaussian multiplier) to the current row-like factor
alpar@9 437 * of matrix H, which corresponds to the transformed row of matrix V. */
alpar@9 438
alpar@9 439 int fhv_update_it(FHV *fhv, int j, int len, const int ind[],
alpar@9 440 const double val[])
alpar@9 441 { int m = fhv->m;
alpar@9 442 LUF *luf = fhv->luf;
alpar@9 443 int *vr_ptr = luf->vr_ptr;
alpar@9 444 int *vr_len = luf->vr_len;
alpar@9 445 int *vr_cap = luf->vr_cap;
alpar@9 446 double *vr_piv = luf->vr_piv;
alpar@9 447 int *vc_ptr = luf->vc_ptr;
alpar@9 448 int *vc_len = luf->vc_len;
alpar@9 449 int *vc_cap = luf->vc_cap;
alpar@9 450 int *pp_row = luf->pp_row;
alpar@9 451 int *pp_col = luf->pp_col;
alpar@9 452 int *qq_row = luf->qq_row;
alpar@9 453 int *qq_col = luf->qq_col;
alpar@9 454 int *sv_ind = luf->sv_ind;
alpar@9 455 double *sv_val = luf->sv_val;
alpar@9 456 double *work = luf->work;
alpar@9 457 double eps_tol = luf->eps_tol;
alpar@9 458 int *hh_ind = fhv->hh_ind;
alpar@9 459 int *hh_ptr = fhv->hh_ptr;
alpar@9 460 int *hh_len = fhv->hh_len;
alpar@9 461 int *p0_row = fhv->p0_row;
alpar@9 462 int *p0_col = fhv->p0_col;
alpar@9 463 int *cc_ind = fhv->cc_ind;
alpar@9 464 double *cc_val = fhv->cc_val;
alpar@9 465 double upd_tol = fhv->upd_tol;
alpar@9 466 int i, i_beg, i_end, i_ptr, j_beg, j_end, j_ptr, k, k1, k2, p, q,
alpar@9 467 p_beg, p_end, p_ptr, ptr, ret;
alpar@9 468 double f, temp;
alpar@9 469 if (!fhv->valid)
alpar@9 470 xfault("fhv_update_it: the factorization is not valid\n");
alpar@9 471 if (!(1 <= j && j <= m))
alpar@9 472 xfault("fhv_update_it: j = %d; column number out of range\n",
alpar@9 473 j);
alpar@9 474 /* check if the new factor of matrix H can be created */
alpar@9 475 if (fhv->hh_nfs == fhv->hh_max)
alpar@9 476 { /* maximal number of updates has been reached */
alpar@9 477 fhv->valid = 0;
alpar@9 478 ret = FHV_ELIMIT;
alpar@9 479 goto done;
alpar@9 480 }
alpar@9 481 /* convert new j-th column of B to dense format */
alpar@9 482 for (i = 1; i <= m; i++)
alpar@9 483 cc_val[i] = 0.0;
alpar@9 484 for (k = 1; k <= len; k++)
alpar@9 485 { i = ind[k];
alpar@9 486 if (!(1 <= i && i <= m))
alpar@9 487 xfault("fhv_update_it: ind[%d] = %d; row number out of rang"
alpar@9 488 "e\n", k, i);
alpar@9 489 if (cc_val[i] != 0.0)
alpar@9 490 xfault("fhv_update_it: ind[%d] = %d; duplicate row index no"
alpar@9 491 "t allowed\n", k, i);
alpar@9 492 if (val[k] == 0.0)
alpar@9 493 xfault("fhv_update_it: val[%d] = %g; zero element not allow"
alpar@9 494 "ed\n", k, val[k]);
alpar@9 495 cc_val[i] = val[k];
alpar@9 496 }
alpar@9 497 /* new j-th column of V := inv(F * H) * (new B[j]) */
alpar@9 498 fhv->luf->pp_row = p0_row;
alpar@9 499 fhv->luf->pp_col = p0_col;
alpar@9 500 luf_f_solve(fhv->luf, 0, cc_val);
alpar@9 501 fhv->luf->pp_row = pp_row;
alpar@9 502 fhv->luf->pp_col = pp_col;
alpar@9 503 fhv_h_solve(fhv, 0, cc_val);
alpar@9 504 /* convert new j-th column of V to sparse format */
alpar@9 505 len = 0;
alpar@9 506 for (i = 1; i <= m; i++)
alpar@9 507 { temp = cc_val[i];
alpar@9 508 if (temp == 0.0 || fabs(temp) < eps_tol) continue;
alpar@9 509 len++, cc_ind[len] = i, cc_val[len] = temp;
alpar@9 510 }
alpar@9 511 /* clear old content of j-th column of matrix V */
alpar@9 512 j_beg = vc_ptr[j];
alpar@9 513 j_end = j_beg + vc_len[j] - 1;
alpar@9 514 for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
alpar@9 515 { /* get row index of v[i,j] */
alpar@9 516 i = sv_ind[j_ptr];
alpar@9 517 /* find v[i,j] in the i-th row */
alpar@9 518 i_beg = vr_ptr[i];
alpar@9 519 i_end = i_beg + vr_len[i] - 1;
alpar@9 520 for (i_ptr = i_beg; sv_ind[i_ptr] != j; i_ptr++) /* nop */;
alpar@9 521 xassert(i_ptr <= i_end);
alpar@9 522 /* remove v[i,j] from the i-th row */
alpar@9 523 sv_ind[i_ptr] = sv_ind[i_end];
alpar@9 524 sv_val[i_ptr] = sv_val[i_end];
alpar@9 525 vr_len[i]--;
alpar@9 526 }
alpar@9 527 /* now j-th column of matrix V is empty */
alpar@9 528 luf->nnz_v -= vc_len[j];
alpar@9 529 vc_len[j] = 0;
alpar@9 530 /* add new elements of j-th column of matrix V to corresponding
alpar@9 531 row lists; determine indices k1 and k2 */
alpar@9 532 k1 = qq_row[j], k2 = 0;
alpar@9 533 for (ptr = 1; ptr <= len; ptr++)
alpar@9 534 { /* get row index of v[i,j] */
alpar@9 535 i = cc_ind[ptr];
alpar@9 536 /* at least one unused location is needed in i-th row */
alpar@9 537 if (vr_len[i] + 1 > vr_cap[i])
alpar@9 538 { if (luf_enlarge_row(luf, i, vr_len[i] + 10))
alpar@9 539 { /* overflow of the sparse vector area */
alpar@9 540 fhv->valid = 0;
alpar@9 541 luf->new_sva = luf->sv_size + luf->sv_size;
alpar@9 542 xassert(luf->new_sva > luf->sv_size);
alpar@9 543 ret = FHV_EROOM;
alpar@9 544 goto done;
alpar@9 545 }
alpar@9 546 }
alpar@9 547 /* add v[i,j] to i-th row */
alpar@9 548 i_ptr = vr_ptr[i] + vr_len[i];
alpar@9 549 sv_ind[i_ptr] = j;
alpar@9 550 sv_val[i_ptr] = cc_val[ptr];
alpar@9 551 vr_len[i]++;
alpar@9 552 /* adjust index k2 */
alpar@9 553 if (k2 < pp_col[i]) k2 = pp_col[i];
alpar@9 554 }
alpar@9 555 /* capacity of j-th column (which is currently empty) should be
alpar@9 556 not less than len locations */
alpar@9 557 if (vc_cap[j] < len)
alpar@9 558 { if (luf_enlarge_col(luf, j, len))
alpar@9 559 { /* overflow of the sparse vector area */
alpar@9 560 fhv->valid = 0;
alpar@9 561 luf->new_sva = luf->sv_size + luf->sv_size;
alpar@9 562 xassert(luf->new_sva > luf->sv_size);
alpar@9 563 ret = FHV_EROOM;
alpar@9 564 goto done;
alpar@9 565 }
alpar@9 566 }
alpar@9 567 /* add new elements of matrix V to j-th column list */
alpar@9 568 j_ptr = vc_ptr[j];
alpar@9 569 memmove(&sv_ind[j_ptr], &cc_ind[1], len * sizeof(int));
alpar@9 570 memmove(&sv_val[j_ptr], &cc_val[1], len * sizeof(double));
alpar@9 571 vc_len[j] = len;
alpar@9 572 luf->nnz_v += len;
alpar@9 573 /* if k1 > k2, diagonal element u[k2,k2] of matrix U is zero and
alpar@9 574 therefore the adjacent basis matrix is structurally singular */
alpar@9 575 if (k1 > k2)
alpar@9 576 { fhv->valid = 0;
alpar@9 577 ret = FHV_ESING;
alpar@9 578 goto done;
alpar@9 579 }
alpar@9 580 /* perform implicit symmetric permutations of rows and columns of
alpar@9 581 matrix U */
alpar@9 582 i = pp_row[k1], j = qq_col[k1];
alpar@9 583 for (k = k1; k < k2; k++)
alpar@9 584 { pp_row[k] = pp_row[k+1], pp_col[pp_row[k]] = k;
alpar@9 585 qq_col[k] = qq_col[k+1], qq_row[qq_col[k]] = k;
alpar@9 586 }
alpar@9 587 pp_row[k2] = i, pp_col[i] = k2;
alpar@9 588 qq_col[k2] = j, qq_row[j] = k2;
alpar@9 589 /* now i-th row of the matrix V is k2-th row of matrix U; since
alpar@9 590 no pivoting is used, only this row will be transformed */
alpar@9 591 /* copy elements of i-th row of matrix V to the working array and
alpar@9 592 remove these elements from matrix V */
alpar@9 593 for (j = 1; j <= m; j++) work[j] = 0.0;
alpar@9 594 i_beg = vr_ptr[i];
alpar@9 595 i_end = i_beg + vr_len[i] - 1;
alpar@9 596 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
alpar@9 597 { /* get column index of v[i,j] */
alpar@9 598 j = sv_ind[i_ptr];
alpar@9 599 /* store v[i,j] to the working array */
alpar@9 600 work[j] = sv_val[i_ptr];
alpar@9 601 /* find v[i,j] in the j-th column */
alpar@9 602 j_beg = vc_ptr[j];
alpar@9 603 j_end = j_beg + vc_len[j] - 1;
alpar@9 604 for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++) /* nop */;
alpar@9 605 xassert(j_ptr <= j_end);
alpar@9 606 /* remove v[i,j] from the j-th column */
alpar@9 607 sv_ind[j_ptr] = sv_ind[j_end];
alpar@9 608 sv_val[j_ptr] = sv_val[j_end];
alpar@9 609 vc_len[j]--;
alpar@9 610 }
alpar@9 611 /* now i-th row of matrix V is empty */
alpar@9 612 luf->nnz_v -= vr_len[i];
alpar@9 613 vr_len[i] = 0;
alpar@9 614 /* create the next row-like factor of the matrix H; this factor
alpar@9 615 corresponds to i-th (transformed) row */
alpar@9 616 fhv->hh_nfs++;
alpar@9 617 hh_ind[fhv->hh_nfs] = i;
alpar@9 618 /* hh_ptr[] will be set later */
alpar@9 619 hh_len[fhv->hh_nfs] = 0;
alpar@9 620 /* up to (k2 - k1) free locations are needed to add new elements
alpar@9 621 to the non-trivial row of the row-like factor */
alpar@9 622 if (luf->sv_end - luf->sv_beg < k2 - k1)
alpar@9 623 { luf_defrag_sva(luf);
alpar@9 624 if (luf->sv_end - luf->sv_beg < k2 - k1)
alpar@9 625 { /* overflow of the sparse vector area */
alpar@9 626 fhv->valid = luf->valid = 0;
alpar@9 627 luf->new_sva = luf->sv_size + luf->sv_size;
alpar@9 628 xassert(luf->new_sva > luf->sv_size);
alpar@9 629 ret = FHV_EROOM;
alpar@9 630 goto done;
alpar@9 631 }
alpar@9 632 }
alpar@9 633 /* eliminate subdiagonal elements of matrix U */
alpar@9 634 for (k = k1; k < k2; k++)
alpar@9 635 { /* v[p,q] = u[k,k] */
alpar@9 636 p = pp_row[k], q = qq_col[k];
alpar@9 637 /* this is the crucial point, where even tiny non-zeros should
alpar@9 638 not be dropped */
alpar@9 639 if (work[q] == 0.0) continue;
alpar@9 640 /* compute gaussian multiplier f = v[i,q] / v[p,q] */
alpar@9 641 f = work[q] / vr_piv[p];
alpar@9 642 /* perform gaussian transformation:
alpar@9 643 (i-th row) := (i-th row) - f * (p-th row)
alpar@9 644 in order to eliminate v[i,q] = u[k2,k] */
alpar@9 645 p_beg = vr_ptr[p];
alpar@9 646 p_end = p_beg + vr_len[p] - 1;
alpar@9 647 for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
alpar@9 648 work[sv_ind[p_ptr]] -= f * sv_val[p_ptr];
alpar@9 649 /* store new element (gaussian multiplier that corresponds to
alpar@9 650 p-th row) in the current row-like factor */
alpar@9 651 luf->sv_end--;
alpar@9 652 sv_ind[luf->sv_end] = p;
alpar@9 653 sv_val[luf->sv_end] = f;
alpar@9 654 hh_len[fhv->hh_nfs]++;
alpar@9 655 }
alpar@9 656 /* set pointer to the current row-like factor of the matrix H
alpar@9 657 (if no elements were added to this factor, it is unity matrix
alpar@9 658 and therefore can be discarded) */
alpar@9 659 if (hh_len[fhv->hh_nfs] == 0)
alpar@9 660 fhv->hh_nfs--;
alpar@9 661 else
alpar@9 662 { hh_ptr[fhv->hh_nfs] = luf->sv_end;
alpar@9 663 fhv->nnz_h += hh_len[fhv->hh_nfs];
alpar@9 664 }
alpar@9 665 /* store new pivot which corresponds to u[k2,k2] */
alpar@9 666 vr_piv[i] = work[qq_col[k2]];
alpar@9 667 /* new elements of i-th row of matrix V (which are non-diagonal
alpar@9 668 elements u[k2,k2+1], ..., u[k2,m] of matrix U = P*V*Q) now are
alpar@9 669 contained in the working array; add them to matrix V */
alpar@9 670 len = 0;
alpar@9 671 for (k = k2+1; k <= m; k++)
alpar@9 672 { /* get column index and value of v[i,j] = u[k2,k] */
alpar@9 673 j = qq_col[k];
alpar@9 674 temp = work[j];
alpar@9 675 /* if v[i,j] is close to zero, skip it */
alpar@9 676 if (fabs(temp) < eps_tol) continue;
alpar@9 677 /* at least one unused location is needed in j-th column */
alpar@9 678 if (vc_len[j] + 1 > vc_cap[j])
alpar@9 679 { if (luf_enlarge_col(luf, j, vc_len[j] + 10))
alpar@9 680 { /* overflow of the sparse vector area */
alpar@9 681 fhv->valid = 0;
alpar@9 682 luf->new_sva = luf->sv_size + luf->sv_size;
alpar@9 683 xassert(luf->new_sva > luf->sv_size);
alpar@9 684 ret = FHV_EROOM;
alpar@9 685 goto done;
alpar@9 686 }
alpar@9 687 }
alpar@9 688 /* add v[i,j] to j-th column */
alpar@9 689 j_ptr = vc_ptr[j] + vc_len[j];
alpar@9 690 sv_ind[j_ptr] = i;
alpar@9 691 sv_val[j_ptr] = temp;
alpar@9 692 vc_len[j]++;
alpar@9 693 /* also store v[i,j] to the auxiliary array */
alpar@9 694 len++, cc_ind[len] = j, cc_val[len] = temp;
alpar@9 695 }
alpar@9 696 /* capacity of i-th row (which is currently empty) should be not
alpar@9 697 less than len locations */
alpar@9 698 if (vr_cap[i] < len)
alpar@9 699 { if (luf_enlarge_row(luf, i, len))
alpar@9 700 { /* overflow of the sparse vector area */
alpar@9 701 fhv->valid = 0;
alpar@9 702 luf->new_sva = luf->sv_size + luf->sv_size;
alpar@9 703 xassert(luf->new_sva > luf->sv_size);
alpar@9 704 ret = FHV_EROOM;
alpar@9 705 goto done;
alpar@9 706 }
alpar@9 707 }
alpar@9 708 /* add new elements to i-th row list */
alpar@9 709 i_ptr = vr_ptr[i];
alpar@9 710 memmove(&sv_ind[i_ptr], &cc_ind[1], len * sizeof(int));
alpar@9 711 memmove(&sv_val[i_ptr], &cc_val[1], len * sizeof(double));
alpar@9 712 vr_len[i] = len;
alpar@9 713 luf->nnz_v += len;
alpar@9 714 /* updating is finished; check that diagonal element u[k2,k2] is
alpar@9 715 not very small in absolute value among other elements in k2-th
alpar@9 716 row and k2-th column of matrix U = P*V*Q */
alpar@9 717 /* temp = max(|u[k2,*]|, |u[*,k2]|) */
alpar@9 718 temp = 0.0;
alpar@9 719 /* walk through k2-th row of U which is i-th row of V */
alpar@9 720 i = pp_row[k2];
alpar@9 721 i_beg = vr_ptr[i];
alpar@9 722 i_end = i_beg + vr_len[i] - 1;
alpar@9 723 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
alpar@9 724 if (temp < fabs(sv_val[i_ptr])) temp = fabs(sv_val[i_ptr]);
alpar@9 725 /* walk through k2-th column of U which is j-th column of V */
alpar@9 726 j = qq_col[k2];
alpar@9 727 j_beg = vc_ptr[j];
alpar@9 728 j_end = j_beg + vc_len[j] - 1;
alpar@9 729 for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
alpar@9 730 if (temp < fabs(sv_val[j_ptr])) temp = fabs(sv_val[j_ptr]);
alpar@9 731 /* check that u[k2,k2] is not very small */
alpar@9 732 if (fabs(vr_piv[i]) < upd_tol * temp)
alpar@9 733 { /* the factorization seems to be inaccurate and therefore must
alpar@9 734 be recomputed */
alpar@9 735 fhv->valid = 0;
alpar@9 736 ret = FHV_ECHECK;
alpar@9 737 goto done;
alpar@9 738 }
alpar@9 739 /* the factorization has been successfully updated */
alpar@9 740 ret = 0;
alpar@9 741 done: /* return to the calling program */
alpar@9 742 return ret;
alpar@9 743 }
alpar@9 744
alpar@9 745 /***********************************************************************
alpar@9 746 * NAME
alpar@9 747 *
alpar@9 748 * fhv_delete_it - delete LP basis factorization
alpar@9 749 *
alpar@9 750 * SYNOPSIS
alpar@9 751 *
alpar@9 752 * #include "glpfhv.h"
alpar@9 753 * void fhv_delete_it(FHV *fhv);
alpar@9 754 *
alpar@9 755 * DESCRIPTION
alpar@9 756 *
alpar@9 757 * The routine fhv_delete_it deletes LP basis factorization specified
alpar@9 758 * by the parameter fhv and frees all memory allocated to this program
alpar@9 759 * object. */
alpar@9 760
alpar@9 761 void fhv_delete_it(FHV *fhv)
alpar@9 762 { luf_delete_it(fhv->luf);
alpar@9 763 if (fhv->hh_ind != NULL) xfree(fhv->hh_ind);
alpar@9 764 if (fhv->hh_ptr != NULL) xfree(fhv->hh_ptr);
alpar@9 765 if (fhv->hh_len != NULL) xfree(fhv->hh_len);
alpar@9 766 if (fhv->p0_row != NULL) xfree(fhv->p0_row);
alpar@9 767 if (fhv->p0_col != NULL) xfree(fhv->p0_col);
alpar@9 768 if (fhv->cc_ind != NULL) xfree(fhv->cc_ind);
alpar@9 769 if (fhv->cc_val != NULL) xfree(fhv->cc_val);
alpar@9 770 xfree(fhv);
alpar@9 771 return;
alpar@9 772 }
alpar@9 773
alpar@9 774 /* eof */