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