[1] | 1 | /* glplpf.c (LP basis factorization, Schur complement 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 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 "glplpf.h" |
---|
| 26 | #include "glpenv.h" |
---|
| 27 | #define xfault xerror |
---|
| 28 | |
---|
| 29 | #define _GLPLPF_DEBUG 0 |
---|
| 30 | |
---|
| 31 | /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ |
---|
| 32 | |
---|
| 33 | #define M_MAX 100000000 /* = 100*10^6 */ |
---|
| 34 | /* maximal order of the basis matrix */ |
---|
| 35 | |
---|
| 36 | /*********************************************************************** |
---|
| 37 | * NAME |
---|
| 38 | * |
---|
| 39 | * lpf_create_it - create LP basis factorization |
---|
| 40 | * |
---|
| 41 | * SYNOPSIS |
---|
| 42 | * |
---|
| 43 | * #include "glplpf.h" |
---|
| 44 | * LPF *lpf_create_it(void); |
---|
| 45 | * |
---|
| 46 | * DESCRIPTION |
---|
| 47 | * |
---|
| 48 | * The routine lpf_create_it creates a program object, which represents |
---|
| 49 | * a factorization of LP basis. |
---|
| 50 | * |
---|
| 51 | * RETURNS |
---|
| 52 | * |
---|
| 53 | * The routine lpf_create_it returns a pointer to the object created. */ |
---|
| 54 | |
---|
| 55 | LPF *lpf_create_it(void) |
---|
| 56 | { LPF *lpf; |
---|
| 57 | #if _GLPLPF_DEBUG |
---|
| 58 | xprintf("lpf_create_it: warning: debug mode enabled\n"); |
---|
| 59 | #endif |
---|
| 60 | lpf = xmalloc(sizeof(LPF)); |
---|
| 61 | lpf->valid = 0; |
---|
| 62 | lpf->m0_max = lpf->m0 = 0; |
---|
| 63 | lpf->luf = luf_create_it(); |
---|
| 64 | lpf->m = 0; |
---|
| 65 | lpf->B = NULL; |
---|
| 66 | lpf->n_max = 50; |
---|
| 67 | lpf->n = 0; |
---|
| 68 | lpf->R_ptr = lpf->R_len = NULL; |
---|
| 69 | lpf->S_ptr = lpf->S_len = NULL; |
---|
| 70 | lpf->scf = NULL; |
---|
| 71 | lpf->P_row = lpf->P_col = NULL; |
---|
| 72 | lpf->Q_row = lpf->Q_col = NULL; |
---|
| 73 | lpf->v_size = 1000; |
---|
| 74 | lpf->v_ptr = 0; |
---|
| 75 | lpf->v_ind = NULL; |
---|
| 76 | lpf->v_val = NULL; |
---|
| 77 | lpf->work1 = lpf->work2 = NULL; |
---|
| 78 | return lpf; |
---|
| 79 | } |
---|
| 80 | |
---|
| 81 | /*********************************************************************** |
---|
| 82 | * NAME |
---|
| 83 | * |
---|
| 84 | * lpf_factorize - compute LP basis factorization |
---|
| 85 | * |
---|
| 86 | * SYNOPSIS |
---|
| 87 | * |
---|
| 88 | * #include "glplpf.h" |
---|
| 89 | * int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col) |
---|
| 90 | * (void *info, int j, int ind[], double val[]), void *info); |
---|
| 91 | * |
---|
| 92 | * DESCRIPTION |
---|
| 93 | * |
---|
| 94 | * The routine lpf_factorize computes the factorization of the basis |
---|
| 95 | * matrix B specified by the routine col. |
---|
| 96 | * |
---|
| 97 | * The parameter lpf specified the basis factorization data structure |
---|
| 98 | * created with the routine lpf_create_it. |
---|
| 99 | * |
---|
| 100 | * The parameter m specifies the order of B, m > 0. |
---|
| 101 | * |
---|
| 102 | * The array bh specifies the basis header: bh[j], 1 <= j <= m, is the |
---|
| 103 | * number of j-th column of B in some original matrix. The array bh is |
---|
| 104 | * optional and can be specified as NULL. |
---|
| 105 | * |
---|
| 106 | * The formal routine col specifies the matrix B to be factorized. To |
---|
| 107 | * obtain j-th column of A the routine lpf_factorize calls the routine |
---|
| 108 | * col with the parameter j (1 <= j <= n). In response the routine col |
---|
| 109 | * should store row indices and numerical values of non-zero elements |
---|
| 110 | * of j-th column of B to locations ind[1,...,len] and val[1,...,len], |
---|
| 111 | * respectively, where len is the number of non-zeros in j-th column |
---|
| 112 | * returned on exit. Neither zero nor duplicate elements are allowed. |
---|
| 113 | * |
---|
| 114 | * The parameter info is a transit pointer passed to the routine col. |
---|
| 115 | * |
---|
| 116 | * RETURNS |
---|
| 117 | * |
---|
| 118 | * 0 The factorization has been successfully computed. |
---|
| 119 | * |
---|
| 120 | * LPF_ESING |
---|
| 121 | * The specified matrix is singular within the working precision. |
---|
| 122 | * |
---|
| 123 | * LPF_ECOND |
---|
| 124 | * The specified matrix is ill-conditioned. |
---|
| 125 | * |
---|
| 126 | * For more details see comments to the routine luf_factorize. */ |
---|
| 127 | |
---|
| 128 | int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col) |
---|
| 129 | (void *info, int j, int ind[], double val[]), void *info) |
---|
| 130 | { int k, ret; |
---|
| 131 | #if _GLPLPF_DEBUG |
---|
| 132 | int i, j, len, *ind; |
---|
| 133 | double *B, *val; |
---|
| 134 | #endif |
---|
| 135 | xassert(bh == bh); |
---|
| 136 | if (m < 1) |
---|
| 137 | xfault("lpf_factorize: m = %d; invalid parameter\n", m); |
---|
| 138 | if (m > M_MAX) |
---|
| 139 | xfault("lpf_factorize: m = %d; matrix too big\n", m); |
---|
| 140 | lpf->m0 = lpf->m = m; |
---|
| 141 | /* invalidate the factorization */ |
---|
| 142 | lpf->valid = 0; |
---|
| 143 | /* allocate/reallocate arrays, if necessary */ |
---|
| 144 | if (lpf->R_ptr == NULL) |
---|
| 145 | lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int)); |
---|
| 146 | if (lpf->R_len == NULL) |
---|
| 147 | lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int)); |
---|
| 148 | if (lpf->S_ptr == NULL) |
---|
| 149 | lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int)); |
---|
| 150 | if (lpf->S_len == NULL) |
---|
| 151 | lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int)); |
---|
| 152 | if (lpf->scf == NULL) |
---|
| 153 | lpf->scf = scf_create_it(lpf->n_max); |
---|
| 154 | if (lpf->v_ind == NULL) |
---|
| 155 | lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int)); |
---|
| 156 | if (lpf->v_val == NULL) |
---|
| 157 | lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double)); |
---|
| 158 | if (lpf->m0_max < m) |
---|
| 159 | { if (lpf->P_row != NULL) xfree(lpf->P_row); |
---|
| 160 | if (lpf->P_col != NULL) xfree(lpf->P_col); |
---|
| 161 | if (lpf->Q_row != NULL) xfree(lpf->Q_row); |
---|
| 162 | if (lpf->Q_col != NULL) xfree(lpf->Q_col); |
---|
| 163 | if (lpf->work1 != NULL) xfree(lpf->work1); |
---|
| 164 | if (lpf->work2 != NULL) xfree(lpf->work2); |
---|
| 165 | lpf->m0_max = m + 100; |
---|
| 166 | lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
---|
| 167 | lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
---|
| 168 | lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
---|
| 169 | lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); |
---|
| 170 | lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double)); |
---|
| 171 | lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double)); |
---|
| 172 | } |
---|
| 173 | /* try to factorize the basis matrix */ |
---|
| 174 | switch (luf_factorize(lpf->luf, m, col, info)) |
---|
| 175 | { case 0: |
---|
| 176 | break; |
---|
| 177 | case LUF_ESING: |
---|
| 178 | ret = LPF_ESING; |
---|
| 179 | goto done; |
---|
| 180 | case LUF_ECOND: |
---|
| 181 | ret = LPF_ECOND; |
---|
| 182 | goto done; |
---|
| 183 | default: |
---|
| 184 | xassert(lpf != lpf); |
---|
| 185 | } |
---|
| 186 | /* the basis matrix has been successfully factorized */ |
---|
| 187 | lpf->valid = 1; |
---|
| 188 | #if _GLPLPF_DEBUG |
---|
| 189 | /* store the basis matrix for debugging */ |
---|
| 190 | if (lpf->B != NULL) xfree(lpf->B); |
---|
| 191 | xassert(m <= 32767); |
---|
| 192 | lpf->B = B = xcalloc(1+m*m, sizeof(double)); |
---|
| 193 | ind = xcalloc(1+m, sizeof(int)); |
---|
| 194 | val = xcalloc(1+m, sizeof(double)); |
---|
| 195 | for (k = 1; k <= m * m; k++) |
---|
| 196 | B[k] = 0.0; |
---|
| 197 | for (j = 1; j <= m; j++) |
---|
| 198 | { len = col(info, j, ind, val); |
---|
| 199 | xassert(0 <= len && len <= m); |
---|
| 200 | for (k = 1; k <= len; k++) |
---|
| 201 | { i = ind[k]; |
---|
| 202 | xassert(1 <= i && i <= m); |
---|
| 203 | xassert(B[(i - 1) * m + j] == 0.0); |
---|
| 204 | xassert(val[k] != 0.0); |
---|
| 205 | B[(i - 1) * m + j] = val[k]; |
---|
| 206 | } |
---|
| 207 | } |
---|
| 208 | xfree(ind); |
---|
| 209 | xfree(val); |
---|
| 210 | #endif |
---|
| 211 | /* B = B0, so there are no additional rows/columns */ |
---|
| 212 | lpf->n = 0; |
---|
| 213 | /* reset the Schur complement factorization */ |
---|
| 214 | scf_reset_it(lpf->scf); |
---|
| 215 | /* P := Q := I */ |
---|
| 216 | for (k = 1; k <= m; k++) |
---|
| 217 | { lpf->P_row[k] = lpf->P_col[k] = k; |
---|
| 218 | lpf->Q_row[k] = lpf->Q_col[k] = k; |
---|
| 219 | } |
---|
| 220 | /* make all SVA locations free */ |
---|
| 221 | lpf->v_ptr = 1; |
---|
| 222 | ret = 0; |
---|
| 223 | done: /* return to the calling program */ |
---|
| 224 | return ret; |
---|
| 225 | } |
---|
| 226 | |
---|
| 227 | /*********************************************************************** |
---|
| 228 | * The routine r_prod computes the product y := y + alpha * R * x, |
---|
| 229 | * where x is a n-vector, alpha is a scalar, y is a m0-vector. |
---|
| 230 | * |
---|
| 231 | * Since matrix R is available by columns, the product is computed as |
---|
| 232 | * a linear combination: |
---|
| 233 | * |
---|
| 234 | * y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]), |
---|
| 235 | * |
---|
| 236 | * where R[j] is j-th column of R. */ |
---|
| 237 | |
---|
| 238 | static void r_prod(LPF *lpf, double y[], double a, const double x[]) |
---|
| 239 | { int n = lpf->n; |
---|
| 240 | int *R_ptr = lpf->R_ptr; |
---|
| 241 | int *R_len = lpf->R_len; |
---|
| 242 | int *v_ind = lpf->v_ind; |
---|
| 243 | double *v_val = lpf->v_val; |
---|
| 244 | int j, beg, end, ptr; |
---|
| 245 | double t; |
---|
| 246 | for (j = 1; j <= n; j++) |
---|
| 247 | { if (x[j] == 0.0) continue; |
---|
| 248 | /* y := y + alpha * R[j] * x[j] */ |
---|
| 249 | t = a * x[j]; |
---|
| 250 | beg = R_ptr[j]; |
---|
| 251 | end = beg + R_len[j]; |
---|
| 252 | for (ptr = beg; ptr < end; ptr++) |
---|
| 253 | y[v_ind[ptr]] += t * v_val[ptr]; |
---|
| 254 | } |
---|
| 255 | return; |
---|
| 256 | } |
---|
| 257 | |
---|
| 258 | /*********************************************************************** |
---|
| 259 | * The routine rt_prod computes the product y := y + alpha * R' * x, |
---|
| 260 | * where R' is a matrix transposed to R, x is a m0-vector, alpha is a |
---|
| 261 | * scalar, y is a n-vector. |
---|
| 262 | * |
---|
| 263 | * Since matrix R is available by columns, the product components are |
---|
| 264 | * computed as inner products: |
---|
| 265 | * |
---|
| 266 | * y[j] := y[j] + alpha * (j-th column of R) * x |
---|
| 267 | * |
---|
| 268 | * for j = 1, 2, ..., n. */ |
---|
| 269 | |
---|
| 270 | static void rt_prod(LPF *lpf, double y[], double a, const double x[]) |
---|
| 271 | { int n = lpf->n; |
---|
| 272 | int *R_ptr = lpf->R_ptr; |
---|
| 273 | int *R_len = lpf->R_len; |
---|
| 274 | int *v_ind = lpf->v_ind; |
---|
| 275 | double *v_val = lpf->v_val; |
---|
| 276 | int j, beg, end, ptr; |
---|
| 277 | double t; |
---|
| 278 | for (j = 1; j <= n; j++) |
---|
| 279 | { /* t := (j-th column of R) * x */ |
---|
| 280 | t = 0.0; |
---|
| 281 | beg = R_ptr[j]; |
---|
| 282 | end = beg + R_len[j]; |
---|
| 283 | for (ptr = beg; ptr < end; ptr++) |
---|
| 284 | t += v_val[ptr] * x[v_ind[ptr]]; |
---|
| 285 | /* y[j] := y[j] + alpha * t */ |
---|
| 286 | y[j] += a * t; |
---|
| 287 | } |
---|
| 288 | return; |
---|
| 289 | } |
---|
| 290 | |
---|
| 291 | /*********************************************************************** |
---|
| 292 | * The routine s_prod computes the product y := y + alpha * S * x, |
---|
| 293 | * where x is a m0-vector, alpha is a scalar, y is a n-vector. |
---|
| 294 | * |
---|
| 295 | * Since matrix S is available by rows, the product components are |
---|
| 296 | * computed as inner products: |
---|
| 297 | * |
---|
| 298 | * y[i] = y[i] + alpha * (i-th row of S) * x |
---|
| 299 | * |
---|
| 300 | * for i = 1, 2, ..., n. */ |
---|
| 301 | |
---|
| 302 | static void s_prod(LPF *lpf, double y[], double a, const double x[]) |
---|
| 303 | { int n = lpf->n; |
---|
| 304 | int *S_ptr = lpf->S_ptr; |
---|
| 305 | int *S_len = lpf->S_len; |
---|
| 306 | int *v_ind = lpf->v_ind; |
---|
| 307 | double *v_val = lpf->v_val; |
---|
| 308 | int i, beg, end, ptr; |
---|
| 309 | double t; |
---|
| 310 | for (i = 1; i <= n; i++) |
---|
| 311 | { /* t := (i-th row of S) * x */ |
---|
| 312 | t = 0.0; |
---|
| 313 | beg = S_ptr[i]; |
---|
| 314 | end = beg + S_len[i]; |
---|
| 315 | for (ptr = beg; ptr < end; ptr++) |
---|
| 316 | t += v_val[ptr] * x[v_ind[ptr]]; |
---|
| 317 | /* y[i] := y[i] + alpha * t */ |
---|
| 318 | y[i] += a * t; |
---|
| 319 | } |
---|
| 320 | return; |
---|
| 321 | } |
---|
| 322 | |
---|
| 323 | /*********************************************************************** |
---|
| 324 | * The routine st_prod computes the product y := y + alpha * S' * x, |
---|
| 325 | * where S' is a matrix transposed to S, x is a n-vector, alpha is a |
---|
| 326 | * scalar, y is m0-vector. |
---|
| 327 | * |
---|
| 328 | * Since matrix R is available by rows, the product is computed as a |
---|
| 329 | * linear combination: |
---|
| 330 | * |
---|
| 331 | * y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]), |
---|
| 332 | * |
---|
| 333 | * where S'[i] is i-th row of S. */ |
---|
| 334 | |
---|
| 335 | static void st_prod(LPF *lpf, double y[], double a, const double x[]) |
---|
| 336 | { int n = lpf->n; |
---|
| 337 | int *S_ptr = lpf->S_ptr; |
---|
| 338 | int *S_len = lpf->S_len; |
---|
| 339 | int *v_ind = lpf->v_ind; |
---|
| 340 | double *v_val = lpf->v_val; |
---|
| 341 | int i, beg, end, ptr; |
---|
| 342 | double t; |
---|
| 343 | for (i = 1; i <= n; i++) |
---|
| 344 | { if (x[i] == 0.0) continue; |
---|
| 345 | /* y := y + alpha * S'[i] * x[i] */ |
---|
| 346 | t = a * x[i]; |
---|
| 347 | beg = S_ptr[i]; |
---|
| 348 | end = beg + S_len[i]; |
---|
| 349 | for (ptr = beg; ptr < end; ptr++) |
---|
| 350 | y[v_ind[ptr]] += t * v_val[ptr]; |
---|
| 351 | } |
---|
| 352 | return; |
---|
| 353 | } |
---|
| 354 | |
---|
| 355 | #if _GLPLPF_DEBUG |
---|
| 356 | /*********************************************************************** |
---|
| 357 | * The routine check_error computes the maximal relative error between |
---|
| 358 | * left- and right-hand sides for the system B * x = b (if tr is zero) |
---|
| 359 | * or B' * x = b (if tr is non-zero), where B' is a matrix transposed |
---|
| 360 | * to B. (This routine is intended for debugging only.) */ |
---|
| 361 | |
---|
| 362 | static void check_error(LPF *lpf, int tr, const double x[], |
---|
| 363 | const double b[]) |
---|
| 364 | { int m = lpf->m; |
---|
| 365 | double *B = lpf->B; |
---|
| 366 | int i, j; |
---|
| 367 | double d, dmax = 0.0, s, t, tmax; |
---|
| 368 | for (i = 1; i <= m; i++) |
---|
| 369 | { s = 0.0; |
---|
| 370 | tmax = 1.0; |
---|
| 371 | for (j = 1; j <= m; j++) |
---|
| 372 | { if (!tr) |
---|
| 373 | t = B[m * (i - 1) + j] * x[j]; |
---|
| 374 | else |
---|
| 375 | t = B[m * (j - 1) + i] * x[j]; |
---|
| 376 | if (tmax < fabs(t)) tmax = fabs(t); |
---|
| 377 | s += t; |
---|
| 378 | } |
---|
| 379 | d = fabs(s - b[i]) / tmax; |
---|
| 380 | if (dmax < d) dmax = d; |
---|
| 381 | } |
---|
| 382 | if (dmax > 1e-8) |
---|
| 383 | xprintf("%s: dmax = %g; relative error too large\n", |
---|
| 384 | !tr ? "lpf_ftran" : "lpf_btran", dmax); |
---|
| 385 | return; |
---|
| 386 | } |
---|
| 387 | #endif |
---|
| 388 | |
---|
| 389 | /*********************************************************************** |
---|
| 390 | * NAME |
---|
| 391 | * |
---|
| 392 | * lpf_ftran - perform forward transformation (solve system B*x = b) |
---|
| 393 | * |
---|
| 394 | * SYNOPSIS |
---|
| 395 | * |
---|
| 396 | * #include "glplpf.h" |
---|
| 397 | * void lpf_ftran(LPF *lpf, double x[]); |
---|
| 398 | * |
---|
| 399 | * DESCRIPTION |
---|
| 400 | * |
---|
| 401 | * The routine lpf_ftran performs forward transformation, i.e. solves |
---|
| 402 | * the system B*x = b, where B is the basis matrix, x is the vector of |
---|
| 403 | * unknowns to be computed, b is the vector of right-hand sides. |
---|
| 404 | * |
---|
| 405 | * On entry elements of the vector b should be stored in dense format |
---|
| 406 | * in locations x[1], ..., x[m], where m is the number of rows. On exit |
---|
| 407 | * the routine stores elements of the vector x in the same locations. |
---|
| 408 | * |
---|
| 409 | * BACKGROUND |
---|
| 410 | * |
---|
| 411 | * Solution of the system B * x = b can be obtained by solving the |
---|
| 412 | * following augmented system: |
---|
| 413 | * |
---|
| 414 | * ( B F^) ( x ) ( b ) |
---|
| 415 | * ( ) ( ) = ( ) |
---|
| 416 | * ( G^ H^) ( y ) ( 0 ) |
---|
| 417 | * |
---|
| 418 | * which, using the main equality, can be written as follows: |
---|
| 419 | * |
---|
| 420 | * ( L0 0 ) ( U0 R ) ( x ) ( b ) |
---|
| 421 | * P ( ) ( ) Q ( ) = ( ) |
---|
| 422 | * ( S I ) ( 0 C ) ( y ) ( 0 ) |
---|
| 423 | * |
---|
| 424 | * therefore, |
---|
| 425 | * |
---|
| 426 | * ( x ) ( U0 R )-1 ( L0 0 )-1 ( b ) |
---|
| 427 | * ( ) = Q' ( ) ( ) P' ( ) |
---|
| 428 | * ( y ) ( 0 C ) ( S I ) ( 0 ) |
---|
| 429 | * |
---|
| 430 | * Thus, computing the solution includes the following steps: |
---|
| 431 | * |
---|
| 432 | * 1. Compute |
---|
| 433 | * |
---|
| 434 | * ( f ) ( b ) |
---|
| 435 | * ( ) = P' ( ) |
---|
| 436 | * ( g ) ( 0 ) |
---|
| 437 | * |
---|
| 438 | * 2. Solve the system |
---|
| 439 | * |
---|
| 440 | * ( f1 ) ( L0 0 )-1 ( f ) ( L0 0 ) ( f1 ) ( f ) |
---|
| 441 | * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
---|
| 442 | * ( g1 ) ( S I ) ( g ) ( S I ) ( g1 ) ( g ) |
---|
| 443 | * |
---|
| 444 | * from which it follows that: |
---|
| 445 | * |
---|
| 446 | * { L0 * f1 = f f1 = inv(L0) * f |
---|
| 447 | * { => |
---|
| 448 | * { S * f1 + g1 = g g1 = g - S * f1 |
---|
| 449 | * |
---|
| 450 | * 3. Solve the system |
---|
| 451 | * |
---|
| 452 | * ( f2 ) ( U0 R )-1 ( f1 ) ( U0 R ) ( f2 ) ( f1 ) |
---|
| 453 | * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
---|
| 454 | * ( g2 ) ( 0 C ) ( g1 ) ( 0 C ) ( g2 ) ( g1 ) |
---|
| 455 | * |
---|
| 456 | * from which it follows that: |
---|
| 457 | * |
---|
| 458 | * { U0 * f2 + R * g2 = f1 f2 = inv(U0) * (f1 - R * g2) |
---|
| 459 | * { => |
---|
| 460 | * { C * g2 = g1 g2 = inv(C) * g1 |
---|
| 461 | * |
---|
| 462 | * 4. Compute |
---|
| 463 | * |
---|
| 464 | * ( x ) ( f2 ) |
---|
| 465 | * ( ) = Q' ( ) |
---|
| 466 | * ( y ) ( g2 ) */ |
---|
| 467 | |
---|
| 468 | void lpf_ftran(LPF *lpf, double x[]) |
---|
| 469 | { int m0 = lpf->m0; |
---|
| 470 | int m = lpf->m; |
---|
| 471 | int n = lpf->n; |
---|
| 472 | int *P_col = lpf->P_col; |
---|
| 473 | int *Q_col = lpf->Q_col; |
---|
| 474 | double *fg = lpf->work1; |
---|
| 475 | double *f = fg; |
---|
| 476 | double *g = fg + m0; |
---|
| 477 | int i, ii; |
---|
| 478 | #if _GLPLPF_DEBUG |
---|
| 479 | double *b; |
---|
| 480 | #endif |
---|
| 481 | if (!lpf->valid) |
---|
| 482 | xfault("lpf_ftran: the factorization is not valid\n"); |
---|
| 483 | xassert(0 <= m && m <= m0 + n); |
---|
| 484 | #if _GLPLPF_DEBUG |
---|
| 485 | /* save the right-hand side vector */ |
---|
| 486 | b = xcalloc(1+m, sizeof(double)); |
---|
| 487 | for (i = 1; i <= m; i++) b[i] = x[i]; |
---|
| 488 | #endif |
---|
| 489 | /* (f g) := inv(P) * (b 0) */ |
---|
| 490 | for (i = 1; i <= m0 + n; i++) |
---|
| 491 | fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0); |
---|
| 492 | /* f1 := inv(L0) * f */ |
---|
| 493 | luf_f_solve(lpf->luf, 0, f); |
---|
| 494 | /* g1 := g - S * f1 */ |
---|
| 495 | s_prod(lpf, g, -1.0, f); |
---|
| 496 | /* g2 := inv(C) * g1 */ |
---|
| 497 | scf_solve_it(lpf->scf, 0, g); |
---|
| 498 | /* f2 := inv(U0) * (f1 - R * g2) */ |
---|
| 499 | r_prod(lpf, f, -1.0, g); |
---|
| 500 | luf_v_solve(lpf->luf, 0, f); |
---|
| 501 | /* (x y) := inv(Q) * (f2 g2) */ |
---|
| 502 | for (i = 1; i <= m; i++) |
---|
| 503 | x[i] = fg[Q_col[i]]; |
---|
| 504 | #if _GLPLPF_DEBUG |
---|
| 505 | /* check relative error in solution */ |
---|
| 506 | check_error(lpf, 0, x, b); |
---|
| 507 | xfree(b); |
---|
| 508 | #endif |
---|
| 509 | return; |
---|
| 510 | } |
---|
| 511 | |
---|
| 512 | /*********************************************************************** |
---|
| 513 | * NAME |
---|
| 514 | * |
---|
| 515 | * lpf_btran - perform backward transformation (solve system B'*x = b) |
---|
| 516 | * |
---|
| 517 | * SYNOPSIS |
---|
| 518 | * |
---|
| 519 | * #include "glplpf.h" |
---|
| 520 | * void lpf_btran(LPF *lpf, double x[]); |
---|
| 521 | * |
---|
| 522 | * DESCRIPTION |
---|
| 523 | * |
---|
| 524 | * The routine lpf_btran performs backward transformation, i.e. solves |
---|
| 525 | * the system B'*x = b, where B' is a matrix transposed to the basis |
---|
| 526 | * matrix B, x is the vector of unknowns to be computed, b is the vector |
---|
| 527 | * of right-hand sides. |
---|
| 528 | * |
---|
| 529 | * On entry elements of the vector b should be stored in dense format |
---|
| 530 | * in locations x[1], ..., x[m], where m is the number of rows. On exit |
---|
| 531 | * the routine stores elements of the vector x in the same locations. |
---|
| 532 | * |
---|
| 533 | * BACKGROUND |
---|
| 534 | * |
---|
| 535 | * Solution of the system B' * x = b, where B' is a matrix transposed |
---|
| 536 | * to B, can be obtained by solving the following augmented system: |
---|
| 537 | * |
---|
| 538 | * ( B F^)T ( x ) ( b ) |
---|
| 539 | * ( ) ( ) = ( ) |
---|
| 540 | * ( G^ H^) ( y ) ( 0 ) |
---|
| 541 | * |
---|
| 542 | * which, using the main equality, can be written as follows: |
---|
| 543 | * |
---|
| 544 | * T ( U0 R )T ( L0 0 )T T ( x ) ( b ) |
---|
| 545 | * Q ( ) ( ) P ( ) = ( ) |
---|
| 546 | * ( 0 C ) ( S I ) ( y ) ( 0 ) |
---|
| 547 | * |
---|
| 548 | * or, equivalently, as follows: |
---|
| 549 | * |
---|
| 550 | * ( U'0 0 ) ( L'0 S') ( x ) ( b ) |
---|
| 551 | * Q' ( ) ( ) P' ( ) = ( ) |
---|
| 552 | * ( R' C') ( 0 I ) ( y ) ( 0 ) |
---|
| 553 | * |
---|
| 554 | * therefore, |
---|
| 555 | * |
---|
| 556 | * ( x ) ( L'0 S')-1 ( U'0 0 )-1 ( b ) |
---|
| 557 | * ( ) = P ( ) ( ) Q ( ) |
---|
| 558 | * ( y ) ( 0 I ) ( R' C') ( 0 ) |
---|
| 559 | * |
---|
| 560 | * Thus, computing the solution includes the following steps: |
---|
| 561 | * |
---|
| 562 | * 1. Compute |
---|
| 563 | * |
---|
| 564 | * ( f ) ( b ) |
---|
| 565 | * ( ) = Q ( ) |
---|
| 566 | * ( g ) ( 0 ) |
---|
| 567 | * |
---|
| 568 | * 2. Solve the system |
---|
| 569 | * |
---|
| 570 | * ( f1 ) ( U'0 0 )-1 ( f ) ( U'0 0 ) ( f1 ) ( f ) |
---|
| 571 | * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
---|
| 572 | * ( g1 ) ( R' C') ( g ) ( R' C') ( g1 ) ( g ) |
---|
| 573 | * |
---|
| 574 | * from which it follows that: |
---|
| 575 | * |
---|
| 576 | * { U'0 * f1 = f f1 = inv(U'0) * f |
---|
| 577 | * { => |
---|
| 578 | * { R' * f1 + C' * g1 = g g1 = inv(C') * (g - R' * f1) |
---|
| 579 | * |
---|
| 580 | * 3. Solve the system |
---|
| 581 | * |
---|
| 582 | * ( f2 ) ( L'0 S')-1 ( f1 ) ( L'0 S') ( f2 ) ( f1 ) |
---|
| 583 | * ( ) = ( ) ( ) => ( ) ( ) = ( ) |
---|
| 584 | * ( g2 ) ( 0 I ) ( g1 ) ( 0 I ) ( g2 ) ( g1 ) |
---|
| 585 | * |
---|
| 586 | * from which it follows that: |
---|
| 587 | * |
---|
| 588 | * { L'0 * f2 + S' * g2 = f1 |
---|
| 589 | * { => f2 = inv(L'0) * ( f1 - S' * g2) |
---|
| 590 | * { g2 = g1 |
---|
| 591 | * |
---|
| 592 | * 4. Compute |
---|
| 593 | * |
---|
| 594 | * ( x ) ( f2 ) |
---|
| 595 | * ( ) = P ( ) |
---|
| 596 | * ( y ) ( g2 ) */ |
---|
| 597 | |
---|
| 598 | void lpf_btran(LPF *lpf, double x[]) |
---|
| 599 | { int m0 = lpf->m0; |
---|
| 600 | int m = lpf->m; |
---|
| 601 | int n = lpf->n; |
---|
| 602 | int *P_row = lpf->P_row; |
---|
| 603 | int *Q_row = lpf->Q_row; |
---|
| 604 | double *fg = lpf->work1; |
---|
| 605 | double *f = fg; |
---|
| 606 | double *g = fg + m0; |
---|
| 607 | int i, ii; |
---|
| 608 | #if _GLPLPF_DEBUG |
---|
| 609 | double *b; |
---|
| 610 | #endif |
---|
| 611 | if (!lpf->valid) |
---|
| 612 | xfault("lpf_btran: the factorization is not valid\n"); |
---|
| 613 | xassert(0 <= m && m <= m0 + n); |
---|
| 614 | #if _GLPLPF_DEBUG |
---|
| 615 | /* save the right-hand side vector */ |
---|
| 616 | b = xcalloc(1+m, sizeof(double)); |
---|
| 617 | for (i = 1; i <= m; i++) b[i] = x[i]; |
---|
| 618 | #endif |
---|
| 619 | /* (f g) := Q * (b 0) */ |
---|
| 620 | for (i = 1; i <= m0 + n; i++) |
---|
| 621 | fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0); |
---|
| 622 | /* f1 := inv(U'0) * f */ |
---|
| 623 | luf_v_solve(lpf->luf, 1, f); |
---|
| 624 | /* g1 := inv(C') * (g - R' * f1) */ |
---|
| 625 | rt_prod(lpf, g, -1.0, f); |
---|
| 626 | scf_solve_it(lpf->scf, 1, g); |
---|
| 627 | /* g2 := g1 */ |
---|
| 628 | g = g; |
---|
| 629 | /* f2 := inv(L'0) * (f1 - S' * g2) */ |
---|
| 630 | st_prod(lpf, f, -1.0, g); |
---|
| 631 | luf_f_solve(lpf->luf, 1, f); |
---|
| 632 | /* (x y) := P * (f2 g2) */ |
---|
| 633 | for (i = 1; i <= m; i++) |
---|
| 634 | x[i] = fg[P_row[i]]; |
---|
| 635 | #if _GLPLPF_DEBUG |
---|
| 636 | /* check relative error in solution */ |
---|
| 637 | check_error(lpf, 1, x, b); |
---|
| 638 | xfree(b); |
---|
| 639 | #endif |
---|
| 640 | return; |
---|
| 641 | } |
---|
| 642 | |
---|
| 643 | /*********************************************************************** |
---|
| 644 | * The routine enlarge_sva enlarges the Sparse Vector Area to new_size |
---|
| 645 | * locations by reallocating the arrays v_ind and v_val. */ |
---|
| 646 | |
---|
| 647 | static void enlarge_sva(LPF *lpf, int new_size) |
---|
| 648 | { int v_size = lpf->v_size; |
---|
| 649 | int used = lpf->v_ptr - 1; |
---|
| 650 | int *v_ind = lpf->v_ind; |
---|
| 651 | double *v_val = lpf->v_val; |
---|
| 652 | xassert(v_size < new_size); |
---|
| 653 | while (v_size < new_size) v_size += v_size; |
---|
| 654 | lpf->v_size = v_size; |
---|
| 655 | lpf->v_ind = xcalloc(1+v_size, sizeof(int)); |
---|
| 656 | lpf->v_val = xcalloc(1+v_size, sizeof(double)); |
---|
| 657 | xassert(used >= 0); |
---|
| 658 | memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int)); |
---|
| 659 | memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double)); |
---|
| 660 | xfree(v_ind); |
---|
| 661 | xfree(v_val); |
---|
| 662 | return; |
---|
| 663 | } |
---|
| 664 | |
---|
| 665 | /*********************************************************************** |
---|
| 666 | * NAME |
---|
| 667 | * |
---|
| 668 | * lpf_update_it - update LP basis factorization |
---|
| 669 | * |
---|
| 670 | * SYNOPSIS |
---|
| 671 | * |
---|
| 672 | * #include "glplpf.h" |
---|
| 673 | * int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[], |
---|
| 674 | * const double val[]); |
---|
| 675 | * |
---|
| 676 | * DESCRIPTION |
---|
| 677 | * |
---|
| 678 | * The routine lpf_update_it updates the factorization of the basis |
---|
| 679 | * matrix B after replacing its j-th column by a new vector. |
---|
| 680 | * |
---|
| 681 | * The parameter j specifies the number of column of B, which has been |
---|
| 682 | * replaced, 1 <= j <= m, where m is the order of B. |
---|
| 683 | * |
---|
| 684 | * The parameter bh specifies the basis header entry for the new column |
---|
| 685 | * of B, which is the number of the new column in some original matrix. |
---|
| 686 | * This parameter is optional and can be specified as 0. |
---|
| 687 | * |
---|
| 688 | * Row indices and numerical values of non-zero elements of the new |
---|
| 689 | * column of B should be placed in locations ind[1], ..., ind[len] and |
---|
| 690 | * val[1], ..., val[len], resp., where len is the number of non-zeros |
---|
| 691 | * in the column. Neither zero nor duplicate elements are allowed. |
---|
| 692 | * |
---|
| 693 | * RETURNS |
---|
| 694 | * |
---|
| 695 | * 0 The factorization has been successfully updated. |
---|
| 696 | * |
---|
| 697 | * LPF_ESING |
---|
| 698 | * New basis B is singular within the working precision. |
---|
| 699 | * |
---|
| 700 | * LPF_ELIMIT |
---|
| 701 | * Maximal number of additional rows and columns has been reached. |
---|
| 702 | * |
---|
| 703 | * BACKGROUND |
---|
| 704 | * |
---|
| 705 | * Let j-th column of the current basis matrix B have to be replaced by |
---|
| 706 | * a new column a. This replacement is equivalent to removing the old |
---|
| 707 | * j-th column by fixing it at zero and introducing the new column as |
---|
| 708 | * follows: |
---|
| 709 | * |
---|
| 710 | * ( B F^| a ) |
---|
| 711 | * ( B F^) ( | ) |
---|
| 712 | * ( ) ---> ( G^ H^| 0 ) |
---|
| 713 | * ( G^ H^) (-------+---) |
---|
| 714 | * ( e'j 0 | 0 ) |
---|
| 715 | * |
---|
| 716 | * where ej is a unit vector with 1 in j-th position which used to fix |
---|
| 717 | * the old j-th column of B (at zero). Then using the main equality we |
---|
| 718 | * have: |
---|
| 719 | * |
---|
| 720 | * ( B F^| a ) ( B0 F | f ) |
---|
| 721 | * ( | ) ( P 0 ) ( | ) ( Q 0 ) |
---|
| 722 | * ( G^ H^| 0 ) = ( ) ( G H | g ) ( ) = |
---|
| 723 | * (-------+---) ( 0 1 ) (-------+---) ( 0 1 ) |
---|
| 724 | * ( e'j 0 | 0 ) ( v' w'| 0 ) |
---|
| 725 | * |
---|
| 726 | * [ ( B0 F )| ( f ) ] [ ( B0 F ) | ( f ) ] |
---|
| 727 | * [ P ( )| P ( ) ] ( Q 0 ) [ P ( ) Q| P ( ) ] |
---|
| 728 | * = [ ( G H )| ( g ) ] ( ) = [ ( G H ) | ( g ) ] |
---|
| 729 | * [------------+-------- ] ( 0 1 ) [-------------+---------] |
---|
| 730 | * [ ( v' w')| 0 ] [ ( v' w') Q| 0 ] |
---|
| 731 | * |
---|
| 732 | * where: |
---|
| 733 | * |
---|
| 734 | * ( a ) ( f ) ( f ) ( a ) |
---|
| 735 | * ( ) = P ( ) => ( ) = P' * ( ) |
---|
| 736 | * ( 0 ) ( g ) ( g ) ( 0 ) |
---|
| 737 | * |
---|
| 738 | * ( ej ) ( v ) ( v ) ( ej ) |
---|
| 739 | * ( e'j 0 ) = ( v' w' ) Q => ( ) = Q' ( ) => ( ) = Q ( ) |
---|
| 740 | * ( 0 ) ( w ) ( w ) ( 0 ) |
---|
| 741 | * |
---|
| 742 | * On the other hand: |
---|
| 743 | * |
---|
| 744 | * ( B0| F f ) |
---|
| 745 | * ( P 0 ) (---+------) ( Q 0 ) ( B0 new F ) |
---|
| 746 | * ( ) ( G | H g ) ( ) = new P ( ) new Q |
---|
| 747 | * ( 0 1 ) ( | ) ( 0 1 ) ( new G new H ) |
---|
| 748 | * ( v'| w' 0 ) |
---|
| 749 | * |
---|
| 750 | * where: |
---|
| 751 | * ( G ) ( H g ) |
---|
| 752 | * new F = ( F f ), new G = ( ), new H = ( ), |
---|
| 753 | * ( v') ( w' 0 ) |
---|
| 754 | * |
---|
| 755 | * ( P 0 ) ( Q 0 ) |
---|
| 756 | * new P = ( ) , new Q = ( ) . |
---|
| 757 | * ( 0 1 ) ( 0 1 ) |
---|
| 758 | * |
---|
| 759 | * The factorization structure for the new augmented matrix remains the |
---|
| 760 | * same, therefore: |
---|
| 761 | * |
---|
| 762 | * ( B0 new F ) ( L0 0 ) ( U0 new R ) |
---|
| 763 | * new P ( ) new Q = ( ) ( ) |
---|
| 764 | * ( new G new H ) ( new S I ) ( 0 new C ) |
---|
| 765 | * |
---|
| 766 | * where: |
---|
| 767 | * |
---|
| 768 | * new F = L0 * new R => |
---|
| 769 | * |
---|
| 770 | * new R = inv(L0) * new F = inv(L0) * (F f) = ( R inv(L0)*f ) |
---|
| 771 | * |
---|
| 772 | * new G = new S * U0 => |
---|
| 773 | * |
---|
| 774 | * ( G ) ( S ) |
---|
| 775 | * new S = new G * inv(U0) = ( ) * inv(U0) = ( ) |
---|
| 776 | * ( v') ( v'*inv(U0) ) |
---|
| 777 | * |
---|
| 778 | * new H = new S * new R + new C => |
---|
| 779 | * |
---|
| 780 | * new C = new H - new S * new R = |
---|
| 781 | * |
---|
| 782 | * ( H g ) ( S ) |
---|
| 783 | * = ( ) - ( ) * ( R inv(L0)*f ) = |
---|
| 784 | * ( w' 0 ) ( v'*inv(U0) ) |
---|
| 785 | * |
---|
| 786 | * ( H - S*R g - S*inv(L0)*f ) ( C x ) |
---|
| 787 | * = ( ) = ( ) |
---|
| 788 | * ( w'- v'*inv(U0)*R -v'*inv(U0)*inv(L0)*f) ( y' z ) |
---|
| 789 | * |
---|
| 790 | * Note that new C is resulted by expanding old C with new column x, |
---|
| 791 | * row y', and diagonal element z, where: |
---|
| 792 | * |
---|
| 793 | * x = g - S * inv(L0) * f = g - S * (new column of R) |
---|
| 794 | * |
---|
| 795 | * y = w - R'* inv(U'0)* v = w - R'* (new row of S) |
---|
| 796 | * |
---|
| 797 | * z = - (new row of S) * (new column of R) |
---|
| 798 | * |
---|
| 799 | * Finally, to replace old B by new B we have to permute j-th and last |
---|
| 800 | * (just added) columns of the matrix |
---|
| 801 | * |
---|
| 802 | * ( B F^| a ) |
---|
| 803 | * ( | ) |
---|
| 804 | * ( G^ H^| 0 ) |
---|
| 805 | * (-------+---) |
---|
| 806 | * ( e'j 0 | 0 ) |
---|
| 807 | * |
---|
| 808 | * and to keep the main equality do the same for matrix Q. */ |
---|
| 809 | |
---|
| 810 | int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[], |
---|
| 811 | const double val[]) |
---|
| 812 | { int m0 = lpf->m0; |
---|
| 813 | int m = lpf->m; |
---|
| 814 | #if _GLPLPF_DEBUG |
---|
| 815 | double *B = lpf->B; |
---|
| 816 | #endif |
---|
| 817 | int n = lpf->n; |
---|
| 818 | int *R_ptr = lpf->R_ptr; |
---|
| 819 | int *R_len = lpf->R_len; |
---|
| 820 | int *S_ptr = lpf->S_ptr; |
---|
| 821 | int *S_len = lpf->S_len; |
---|
| 822 | int *P_row = lpf->P_row; |
---|
| 823 | int *P_col = lpf->P_col; |
---|
| 824 | int *Q_row = lpf->Q_row; |
---|
| 825 | int *Q_col = lpf->Q_col; |
---|
| 826 | int v_ptr = lpf->v_ptr; |
---|
| 827 | int *v_ind = lpf->v_ind; |
---|
| 828 | double *v_val = lpf->v_val; |
---|
| 829 | double *a = lpf->work2; /* new column */ |
---|
| 830 | double *fg = lpf->work1, *f = fg, *g = fg + m0; |
---|
| 831 | double *vw = lpf->work2, *v = vw, *w = vw + m0; |
---|
| 832 | double *x = g, *y = w, z; |
---|
| 833 | int i, ii, k, ret; |
---|
| 834 | xassert(bh == bh); |
---|
| 835 | if (!lpf->valid) |
---|
| 836 | xfault("lpf_update_it: the factorization is not valid\n"); |
---|
| 837 | if (!(1 <= j && j <= m)) |
---|
| 838 | xfault("lpf_update_it: j = %d; column number out of range\n", |
---|
| 839 | j); |
---|
| 840 | xassert(0 <= m && m <= m0 + n); |
---|
| 841 | /* check if the basis factorization can be expanded */ |
---|
| 842 | if (n == lpf->n_max) |
---|
| 843 | { lpf->valid = 0; |
---|
| 844 | ret = LPF_ELIMIT; |
---|
| 845 | goto done; |
---|
| 846 | } |
---|
| 847 | /* convert new j-th column of B to dense format */ |
---|
| 848 | for (i = 1; i <= m; i++) |
---|
| 849 | a[i] = 0.0; |
---|
| 850 | for (k = 1; k <= len; k++) |
---|
| 851 | { i = ind[k]; |
---|
| 852 | if (!(1 <= i && i <= m)) |
---|
| 853 | xfault("lpf_update_it: ind[%d] = %d; row number out of rang" |
---|
| 854 | "e\n", k, i); |
---|
| 855 | if (a[i] != 0.0) |
---|
| 856 | xfault("lpf_update_it: ind[%d] = %d; duplicate row index no" |
---|
| 857 | "t allowed\n", k, i); |
---|
| 858 | if (val[k] == 0.0) |
---|
| 859 | xfault("lpf_update_it: val[%d] = %g; zero element not allow" |
---|
| 860 | "ed\n", k, val[k]); |
---|
| 861 | a[i] = val[k]; |
---|
| 862 | } |
---|
| 863 | #if _GLPLPF_DEBUG |
---|
| 864 | /* change column in the basis matrix for debugging */ |
---|
| 865 | for (i = 1; i <= m; i++) |
---|
| 866 | B[(i - 1) * m + j] = a[i]; |
---|
| 867 | #endif |
---|
| 868 | /* (f g) := inv(P) * (a 0) */ |
---|
| 869 | for (i = 1; i <= m0+n; i++) |
---|
| 870 | fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0); |
---|
| 871 | /* (v w) := Q * (ej 0) */ |
---|
| 872 | for (i = 1; i <= m0+n; i++) vw[i] = 0.0; |
---|
| 873 | vw[Q_col[j]] = 1.0; |
---|
| 874 | /* f1 := inv(L0) * f (new column of R) */ |
---|
| 875 | luf_f_solve(lpf->luf, 0, f); |
---|
| 876 | /* v1 := inv(U'0) * v (new row of S) */ |
---|
| 877 | luf_v_solve(lpf->luf, 1, v); |
---|
| 878 | /* we need at most 2 * m0 available locations in the SVA to store |
---|
| 879 | new column of matrix R and new row of matrix S */ |
---|
| 880 | if (lpf->v_size < v_ptr + m0 + m0) |
---|
| 881 | { enlarge_sva(lpf, v_ptr + m0 + m0); |
---|
| 882 | v_ind = lpf->v_ind; |
---|
| 883 | v_val = lpf->v_val; |
---|
| 884 | } |
---|
| 885 | /* store new column of R */ |
---|
| 886 | R_ptr[n+1] = v_ptr; |
---|
| 887 | for (i = 1; i <= m0; i++) |
---|
| 888 | { if (f[i] != 0.0) |
---|
| 889 | v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++; |
---|
| 890 | } |
---|
| 891 | R_len[n+1] = v_ptr - lpf->v_ptr; |
---|
| 892 | lpf->v_ptr = v_ptr; |
---|
| 893 | /* store new row of S */ |
---|
| 894 | S_ptr[n+1] = v_ptr; |
---|
| 895 | for (i = 1; i <= m0; i++) |
---|
| 896 | { if (v[i] != 0.0) |
---|
| 897 | v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++; |
---|
| 898 | } |
---|
| 899 | S_len[n+1] = v_ptr - lpf->v_ptr; |
---|
| 900 | lpf->v_ptr = v_ptr; |
---|
| 901 | /* x := g - S * f1 (new column of C) */ |
---|
| 902 | s_prod(lpf, x, -1.0, f); |
---|
| 903 | /* y := w - R' * v1 (new row of C) */ |
---|
| 904 | rt_prod(lpf, y, -1.0, v); |
---|
| 905 | /* z := - v1 * f1 (new diagonal element of C) */ |
---|
| 906 | z = 0.0; |
---|
| 907 | for (i = 1; i <= m0; i++) z -= v[i] * f[i]; |
---|
| 908 | /* update factorization of new matrix C */ |
---|
| 909 | switch (scf_update_exp(lpf->scf, x, y, z)) |
---|
| 910 | { case 0: |
---|
| 911 | break; |
---|
| 912 | case SCF_ESING: |
---|
| 913 | lpf->valid = 0; |
---|
| 914 | ret = LPF_ESING; |
---|
| 915 | goto done; |
---|
| 916 | case SCF_ELIMIT: |
---|
| 917 | xassert(lpf != lpf); |
---|
| 918 | default: |
---|
| 919 | xassert(lpf != lpf); |
---|
| 920 | } |
---|
| 921 | /* expand matrix P */ |
---|
| 922 | P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1; |
---|
| 923 | /* expand matrix Q */ |
---|
| 924 | Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1; |
---|
| 925 | /* permute j-th and last (just added) column of matrix Q */ |
---|
| 926 | i = Q_col[j], ii = Q_col[m0+n+1]; |
---|
| 927 | Q_row[i] = m0+n+1, Q_col[m0+n+1] = i; |
---|
| 928 | Q_row[ii] = j, Q_col[j] = ii; |
---|
| 929 | /* increase the number of additional rows and columns */ |
---|
| 930 | lpf->n++; |
---|
| 931 | xassert(lpf->n <= lpf->n_max); |
---|
| 932 | /* the factorization has been successfully updated */ |
---|
| 933 | ret = 0; |
---|
| 934 | done: /* return to the calling program */ |
---|
| 935 | return ret; |
---|
| 936 | } |
---|
| 937 | |
---|
| 938 | /*********************************************************************** |
---|
| 939 | * NAME |
---|
| 940 | * |
---|
| 941 | * lpf_delete_it - delete LP basis factorization |
---|
| 942 | * |
---|
| 943 | * SYNOPSIS |
---|
| 944 | * |
---|
| 945 | * #include "glplpf.h" |
---|
| 946 | * void lpf_delete_it(LPF *lpf) |
---|
| 947 | * |
---|
| 948 | * DESCRIPTION |
---|
| 949 | * |
---|
| 950 | * The routine lpf_delete_it deletes LP basis factorization specified |
---|
| 951 | * by the parameter lpf and frees all memory allocated to this program |
---|
| 952 | * object. */ |
---|
| 953 | |
---|
| 954 | void lpf_delete_it(LPF *lpf) |
---|
| 955 | { luf_delete_it(lpf->luf); |
---|
| 956 | #if _GLPLPF_DEBUG |
---|
| 957 | if (lpf->B != NULL) xfree(lpf->B); |
---|
| 958 | #else |
---|
| 959 | xassert(lpf->B == NULL); |
---|
| 960 | #endif |
---|
| 961 | if (lpf->R_ptr != NULL) xfree(lpf->R_ptr); |
---|
| 962 | if (lpf->R_len != NULL) xfree(lpf->R_len); |
---|
| 963 | if (lpf->S_ptr != NULL) xfree(lpf->S_ptr); |
---|
| 964 | if (lpf->S_len != NULL) xfree(lpf->S_len); |
---|
| 965 | if (lpf->scf != NULL) scf_delete_it(lpf->scf); |
---|
| 966 | if (lpf->P_row != NULL) xfree(lpf->P_row); |
---|
| 967 | if (lpf->P_col != NULL) xfree(lpf->P_col); |
---|
| 968 | if (lpf->Q_row != NULL) xfree(lpf->Q_row); |
---|
| 969 | if (lpf->Q_col != NULL) xfree(lpf->Q_col); |
---|
| 970 | if (lpf->v_ind != NULL) xfree(lpf->v_ind); |
---|
| 971 | if (lpf->v_val != NULL) xfree(lpf->v_val); |
---|
| 972 | if (lpf->work1 != NULL) xfree(lpf->work1); |
---|
| 973 | if (lpf->work2 != NULL) xfree(lpf->work2); |
---|
| 974 | xfree(lpf); |
---|
| 975 | return; |
---|
| 976 | } |
---|
| 977 | |
---|
| 978 | /* eof */ |
---|