1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glplpf.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,978 @@
1.4 +/* glplpf.c (LP basis factorization, Schur complement version) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glplpf.h"
1.29 +#include "glpenv.h"
1.30 +#define xfault xerror
1.31 +
1.32 +#define _GLPLPF_DEBUG 0
1.33 +
1.34 +/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
1.35 +
1.36 +#define M_MAX 100000000 /* = 100*10^6 */
1.37 +/* maximal order of the basis matrix */
1.38 +
1.39 +/***********************************************************************
1.40 +* NAME
1.41 +*
1.42 +* lpf_create_it - create LP basis factorization
1.43 +*
1.44 +* SYNOPSIS
1.45 +*
1.46 +* #include "glplpf.h"
1.47 +* LPF *lpf_create_it(void);
1.48 +*
1.49 +* DESCRIPTION
1.50 +*
1.51 +* The routine lpf_create_it creates a program object, which represents
1.52 +* a factorization of LP basis.
1.53 +*
1.54 +* RETURNS
1.55 +*
1.56 +* The routine lpf_create_it returns a pointer to the object created. */
1.57 +
1.58 +LPF *lpf_create_it(void)
1.59 +{ LPF *lpf;
1.60 +#if _GLPLPF_DEBUG
1.61 + xprintf("lpf_create_it: warning: debug mode enabled\n");
1.62 +#endif
1.63 + lpf = xmalloc(sizeof(LPF));
1.64 + lpf->valid = 0;
1.65 + lpf->m0_max = lpf->m0 = 0;
1.66 + lpf->luf = luf_create_it();
1.67 + lpf->m = 0;
1.68 + lpf->B = NULL;
1.69 + lpf->n_max = 50;
1.70 + lpf->n = 0;
1.71 + lpf->R_ptr = lpf->R_len = NULL;
1.72 + lpf->S_ptr = lpf->S_len = NULL;
1.73 + lpf->scf = NULL;
1.74 + lpf->P_row = lpf->P_col = NULL;
1.75 + lpf->Q_row = lpf->Q_col = NULL;
1.76 + lpf->v_size = 1000;
1.77 + lpf->v_ptr = 0;
1.78 + lpf->v_ind = NULL;
1.79 + lpf->v_val = NULL;
1.80 + lpf->work1 = lpf->work2 = NULL;
1.81 + return lpf;
1.82 +}
1.83 +
1.84 +/***********************************************************************
1.85 +* NAME
1.86 +*
1.87 +* lpf_factorize - compute LP basis factorization
1.88 +*
1.89 +* SYNOPSIS
1.90 +*
1.91 +* #include "glplpf.h"
1.92 +* int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
1.93 +* (void *info, int j, int ind[], double val[]), void *info);
1.94 +*
1.95 +* DESCRIPTION
1.96 +*
1.97 +* The routine lpf_factorize computes the factorization of the basis
1.98 +* matrix B specified by the routine col.
1.99 +*
1.100 +* The parameter lpf specified the basis factorization data structure
1.101 +* created with the routine lpf_create_it.
1.102 +*
1.103 +* The parameter m specifies the order of B, m > 0.
1.104 +*
1.105 +* The array bh specifies the basis header: bh[j], 1 <= j <= m, is the
1.106 +* number of j-th column of B in some original matrix. The array bh is
1.107 +* optional and can be specified as NULL.
1.108 +*
1.109 +* The formal routine col specifies the matrix B to be factorized. To
1.110 +* obtain j-th column of A the routine lpf_factorize calls the routine
1.111 +* col with the parameter j (1 <= j <= n). In response the routine col
1.112 +* should store row indices and numerical values of non-zero elements
1.113 +* of j-th column of B to locations ind[1,...,len] and val[1,...,len],
1.114 +* respectively, where len is the number of non-zeros in j-th column
1.115 +* returned on exit. Neither zero nor duplicate elements are allowed.
1.116 +*
1.117 +* The parameter info is a transit pointer passed to the routine col.
1.118 +*
1.119 +* RETURNS
1.120 +*
1.121 +* 0 The factorization has been successfully computed.
1.122 +*
1.123 +* LPF_ESING
1.124 +* The specified matrix is singular within the working precision.
1.125 +*
1.126 +* LPF_ECOND
1.127 +* The specified matrix is ill-conditioned.
1.128 +*
1.129 +* For more details see comments to the routine luf_factorize. */
1.130 +
1.131 +int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col)
1.132 + (void *info, int j, int ind[], double val[]), void *info)
1.133 +{ int k, ret;
1.134 +#if _GLPLPF_DEBUG
1.135 + int i, j, len, *ind;
1.136 + double *B, *val;
1.137 +#endif
1.138 + xassert(bh == bh);
1.139 + if (m < 1)
1.140 + xfault("lpf_factorize: m = %d; invalid parameter\n", m);
1.141 + if (m > M_MAX)
1.142 + xfault("lpf_factorize: m = %d; matrix too big\n", m);
1.143 + lpf->m0 = lpf->m = m;
1.144 + /* invalidate the factorization */
1.145 + lpf->valid = 0;
1.146 + /* allocate/reallocate arrays, if necessary */
1.147 + if (lpf->R_ptr == NULL)
1.148 + lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int));
1.149 + if (lpf->R_len == NULL)
1.150 + lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int));
1.151 + if (lpf->S_ptr == NULL)
1.152 + lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int));
1.153 + if (lpf->S_len == NULL)
1.154 + lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int));
1.155 + if (lpf->scf == NULL)
1.156 + lpf->scf = scf_create_it(lpf->n_max);
1.157 + if (lpf->v_ind == NULL)
1.158 + lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int));
1.159 + if (lpf->v_val == NULL)
1.160 + lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double));
1.161 + if (lpf->m0_max < m)
1.162 + { if (lpf->P_row != NULL) xfree(lpf->P_row);
1.163 + if (lpf->P_col != NULL) xfree(lpf->P_col);
1.164 + if (lpf->Q_row != NULL) xfree(lpf->Q_row);
1.165 + if (lpf->Q_col != NULL) xfree(lpf->Q_col);
1.166 + if (lpf->work1 != NULL) xfree(lpf->work1);
1.167 + if (lpf->work2 != NULL) xfree(lpf->work2);
1.168 + lpf->m0_max = m + 100;
1.169 + lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
1.170 + lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
1.171 + lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
1.172 + lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int));
1.173 + lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
1.174 + lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double));
1.175 + }
1.176 + /* try to factorize the basis matrix */
1.177 + switch (luf_factorize(lpf->luf, m, col, info))
1.178 + { case 0:
1.179 + break;
1.180 + case LUF_ESING:
1.181 + ret = LPF_ESING;
1.182 + goto done;
1.183 + case LUF_ECOND:
1.184 + ret = LPF_ECOND;
1.185 + goto done;
1.186 + default:
1.187 + xassert(lpf != lpf);
1.188 + }
1.189 + /* the basis matrix has been successfully factorized */
1.190 + lpf->valid = 1;
1.191 +#if _GLPLPF_DEBUG
1.192 + /* store the basis matrix for debugging */
1.193 + if (lpf->B != NULL) xfree(lpf->B);
1.194 + xassert(m <= 32767);
1.195 + lpf->B = B = xcalloc(1+m*m, sizeof(double));
1.196 + ind = xcalloc(1+m, sizeof(int));
1.197 + val = xcalloc(1+m, sizeof(double));
1.198 + for (k = 1; k <= m * m; k++)
1.199 + B[k] = 0.0;
1.200 + for (j = 1; j <= m; j++)
1.201 + { len = col(info, j, ind, val);
1.202 + xassert(0 <= len && len <= m);
1.203 + for (k = 1; k <= len; k++)
1.204 + { i = ind[k];
1.205 + xassert(1 <= i && i <= m);
1.206 + xassert(B[(i - 1) * m + j] == 0.0);
1.207 + xassert(val[k] != 0.0);
1.208 + B[(i - 1) * m + j] = val[k];
1.209 + }
1.210 + }
1.211 + xfree(ind);
1.212 + xfree(val);
1.213 +#endif
1.214 + /* B = B0, so there are no additional rows/columns */
1.215 + lpf->n = 0;
1.216 + /* reset the Schur complement factorization */
1.217 + scf_reset_it(lpf->scf);
1.218 + /* P := Q := I */
1.219 + for (k = 1; k <= m; k++)
1.220 + { lpf->P_row[k] = lpf->P_col[k] = k;
1.221 + lpf->Q_row[k] = lpf->Q_col[k] = k;
1.222 + }
1.223 + /* make all SVA locations free */
1.224 + lpf->v_ptr = 1;
1.225 + ret = 0;
1.226 +done: /* return to the calling program */
1.227 + return ret;
1.228 +}
1.229 +
1.230 +/***********************************************************************
1.231 +* The routine r_prod computes the product y := y + alpha * R * x,
1.232 +* where x is a n-vector, alpha is a scalar, y is a m0-vector.
1.233 +*
1.234 +* Since matrix R is available by columns, the product is computed as
1.235 +* a linear combination:
1.236 +*
1.237 +* y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]),
1.238 +*
1.239 +* where R[j] is j-th column of R. */
1.240 +
1.241 +static void r_prod(LPF *lpf, double y[], double a, const double x[])
1.242 +{ int n = lpf->n;
1.243 + int *R_ptr = lpf->R_ptr;
1.244 + int *R_len = lpf->R_len;
1.245 + int *v_ind = lpf->v_ind;
1.246 + double *v_val = lpf->v_val;
1.247 + int j, beg, end, ptr;
1.248 + double t;
1.249 + for (j = 1; j <= n; j++)
1.250 + { if (x[j] == 0.0) continue;
1.251 + /* y := y + alpha * R[j] * x[j] */
1.252 + t = a * x[j];
1.253 + beg = R_ptr[j];
1.254 + end = beg + R_len[j];
1.255 + for (ptr = beg; ptr < end; ptr++)
1.256 + y[v_ind[ptr]] += t * v_val[ptr];
1.257 + }
1.258 + return;
1.259 +}
1.260 +
1.261 +/***********************************************************************
1.262 +* The routine rt_prod computes the product y := y + alpha * R' * x,
1.263 +* where R' is a matrix transposed to R, x is a m0-vector, alpha is a
1.264 +* scalar, y is a n-vector.
1.265 +*
1.266 +* Since matrix R is available by columns, the product components are
1.267 +* computed as inner products:
1.268 +*
1.269 +* y[j] := y[j] + alpha * (j-th column of R) * x
1.270 +*
1.271 +* for j = 1, 2, ..., n. */
1.272 +
1.273 +static void rt_prod(LPF *lpf, double y[], double a, const double x[])
1.274 +{ int n = lpf->n;
1.275 + int *R_ptr = lpf->R_ptr;
1.276 + int *R_len = lpf->R_len;
1.277 + int *v_ind = lpf->v_ind;
1.278 + double *v_val = lpf->v_val;
1.279 + int j, beg, end, ptr;
1.280 + double t;
1.281 + for (j = 1; j <= n; j++)
1.282 + { /* t := (j-th column of R) * x */
1.283 + t = 0.0;
1.284 + beg = R_ptr[j];
1.285 + end = beg + R_len[j];
1.286 + for (ptr = beg; ptr < end; ptr++)
1.287 + t += v_val[ptr] * x[v_ind[ptr]];
1.288 + /* y[j] := y[j] + alpha * t */
1.289 + y[j] += a * t;
1.290 + }
1.291 + return;
1.292 +}
1.293 +
1.294 +/***********************************************************************
1.295 +* The routine s_prod computes the product y := y + alpha * S * x,
1.296 +* where x is a m0-vector, alpha is a scalar, y is a n-vector.
1.297 +*
1.298 +* Since matrix S is available by rows, the product components are
1.299 +* computed as inner products:
1.300 +*
1.301 +* y[i] = y[i] + alpha * (i-th row of S) * x
1.302 +*
1.303 +* for i = 1, 2, ..., n. */
1.304 +
1.305 +static void s_prod(LPF *lpf, double y[], double a, const double x[])
1.306 +{ int n = lpf->n;
1.307 + int *S_ptr = lpf->S_ptr;
1.308 + int *S_len = lpf->S_len;
1.309 + int *v_ind = lpf->v_ind;
1.310 + double *v_val = lpf->v_val;
1.311 + int i, beg, end, ptr;
1.312 + double t;
1.313 + for (i = 1; i <= n; i++)
1.314 + { /* t := (i-th row of S) * x */
1.315 + t = 0.0;
1.316 + beg = S_ptr[i];
1.317 + end = beg + S_len[i];
1.318 + for (ptr = beg; ptr < end; ptr++)
1.319 + t += v_val[ptr] * x[v_ind[ptr]];
1.320 + /* y[i] := y[i] + alpha * t */
1.321 + y[i] += a * t;
1.322 + }
1.323 + return;
1.324 +}
1.325 +
1.326 +/***********************************************************************
1.327 +* The routine st_prod computes the product y := y + alpha * S' * x,
1.328 +* where S' is a matrix transposed to S, x is a n-vector, alpha is a
1.329 +* scalar, y is m0-vector.
1.330 +*
1.331 +* Since matrix R is available by rows, the product is computed as a
1.332 +* linear combination:
1.333 +*
1.334 +* y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]),
1.335 +*
1.336 +* where S'[i] is i-th row of S. */
1.337 +
1.338 +static void st_prod(LPF *lpf, double y[], double a, const double x[])
1.339 +{ int n = lpf->n;
1.340 + int *S_ptr = lpf->S_ptr;
1.341 + int *S_len = lpf->S_len;
1.342 + int *v_ind = lpf->v_ind;
1.343 + double *v_val = lpf->v_val;
1.344 + int i, beg, end, ptr;
1.345 + double t;
1.346 + for (i = 1; i <= n; i++)
1.347 + { if (x[i] == 0.0) continue;
1.348 + /* y := y + alpha * S'[i] * x[i] */
1.349 + t = a * x[i];
1.350 + beg = S_ptr[i];
1.351 + end = beg + S_len[i];
1.352 + for (ptr = beg; ptr < end; ptr++)
1.353 + y[v_ind[ptr]] += t * v_val[ptr];
1.354 + }
1.355 + return;
1.356 +}
1.357 +
1.358 +#if _GLPLPF_DEBUG
1.359 +/***********************************************************************
1.360 +* The routine check_error computes the maximal relative error between
1.361 +* left- and right-hand sides for the system B * x = b (if tr is zero)
1.362 +* or B' * x = b (if tr is non-zero), where B' is a matrix transposed
1.363 +* to B. (This routine is intended for debugging only.) */
1.364 +
1.365 +static void check_error(LPF *lpf, int tr, const double x[],
1.366 + const double b[])
1.367 +{ int m = lpf->m;
1.368 + double *B = lpf->B;
1.369 + int i, j;
1.370 + double d, dmax = 0.0, s, t, tmax;
1.371 + for (i = 1; i <= m; i++)
1.372 + { s = 0.0;
1.373 + tmax = 1.0;
1.374 + for (j = 1; j <= m; j++)
1.375 + { if (!tr)
1.376 + t = B[m * (i - 1) + j] * x[j];
1.377 + else
1.378 + t = B[m * (j - 1) + i] * x[j];
1.379 + if (tmax < fabs(t)) tmax = fabs(t);
1.380 + s += t;
1.381 + }
1.382 + d = fabs(s - b[i]) / tmax;
1.383 + if (dmax < d) dmax = d;
1.384 + }
1.385 + if (dmax > 1e-8)
1.386 + xprintf("%s: dmax = %g; relative error too large\n",
1.387 + !tr ? "lpf_ftran" : "lpf_btran", dmax);
1.388 + return;
1.389 +}
1.390 +#endif
1.391 +
1.392 +/***********************************************************************
1.393 +* NAME
1.394 +*
1.395 +* lpf_ftran - perform forward transformation (solve system B*x = b)
1.396 +*
1.397 +* SYNOPSIS
1.398 +*
1.399 +* #include "glplpf.h"
1.400 +* void lpf_ftran(LPF *lpf, double x[]);
1.401 +*
1.402 +* DESCRIPTION
1.403 +*
1.404 +* The routine lpf_ftran performs forward transformation, i.e. solves
1.405 +* the system B*x = b, where B is the basis matrix, x is the vector of
1.406 +* unknowns to be computed, b is the vector of right-hand sides.
1.407 +*
1.408 +* On entry elements of the vector b should be stored in dense format
1.409 +* in locations x[1], ..., x[m], where m is the number of rows. On exit
1.410 +* the routine stores elements of the vector x in the same locations.
1.411 +*
1.412 +* BACKGROUND
1.413 +*
1.414 +* Solution of the system B * x = b can be obtained by solving the
1.415 +* following augmented system:
1.416 +*
1.417 +* ( B F^) ( x ) ( b )
1.418 +* ( ) ( ) = ( )
1.419 +* ( G^ H^) ( y ) ( 0 )
1.420 +*
1.421 +* which, using the main equality, can be written as follows:
1.422 +*
1.423 +* ( L0 0 ) ( U0 R ) ( x ) ( b )
1.424 +* P ( ) ( ) Q ( ) = ( )
1.425 +* ( S I ) ( 0 C ) ( y ) ( 0 )
1.426 +*
1.427 +* therefore,
1.428 +*
1.429 +* ( x ) ( U0 R )-1 ( L0 0 )-1 ( b )
1.430 +* ( ) = Q' ( ) ( ) P' ( )
1.431 +* ( y ) ( 0 C ) ( S I ) ( 0 )
1.432 +*
1.433 +* Thus, computing the solution includes the following steps:
1.434 +*
1.435 +* 1. Compute
1.436 +*
1.437 +* ( f ) ( b )
1.438 +* ( ) = P' ( )
1.439 +* ( g ) ( 0 )
1.440 +*
1.441 +* 2. Solve the system
1.442 +*
1.443 +* ( f1 ) ( L0 0 )-1 ( f ) ( L0 0 ) ( f1 ) ( f )
1.444 +* ( ) = ( ) ( ) => ( ) ( ) = ( )
1.445 +* ( g1 ) ( S I ) ( g ) ( S I ) ( g1 ) ( g )
1.446 +*
1.447 +* from which it follows that:
1.448 +*
1.449 +* { L0 * f1 = f f1 = inv(L0) * f
1.450 +* { =>
1.451 +* { S * f1 + g1 = g g1 = g - S * f1
1.452 +*
1.453 +* 3. Solve the system
1.454 +*
1.455 +* ( f2 ) ( U0 R )-1 ( f1 ) ( U0 R ) ( f2 ) ( f1 )
1.456 +* ( ) = ( ) ( ) => ( ) ( ) = ( )
1.457 +* ( g2 ) ( 0 C ) ( g1 ) ( 0 C ) ( g2 ) ( g1 )
1.458 +*
1.459 +* from which it follows that:
1.460 +*
1.461 +* { U0 * f2 + R * g2 = f1 f2 = inv(U0) * (f1 - R * g2)
1.462 +* { =>
1.463 +* { C * g2 = g1 g2 = inv(C) * g1
1.464 +*
1.465 +* 4. Compute
1.466 +*
1.467 +* ( x ) ( f2 )
1.468 +* ( ) = Q' ( )
1.469 +* ( y ) ( g2 ) */
1.470 +
1.471 +void lpf_ftran(LPF *lpf, double x[])
1.472 +{ int m0 = lpf->m0;
1.473 + int m = lpf->m;
1.474 + int n = lpf->n;
1.475 + int *P_col = lpf->P_col;
1.476 + int *Q_col = lpf->Q_col;
1.477 + double *fg = lpf->work1;
1.478 + double *f = fg;
1.479 + double *g = fg + m0;
1.480 + int i, ii;
1.481 +#if _GLPLPF_DEBUG
1.482 + double *b;
1.483 +#endif
1.484 + if (!lpf->valid)
1.485 + xfault("lpf_ftran: the factorization is not valid\n");
1.486 + xassert(0 <= m && m <= m0 + n);
1.487 +#if _GLPLPF_DEBUG
1.488 + /* save the right-hand side vector */
1.489 + b = xcalloc(1+m, sizeof(double));
1.490 + for (i = 1; i <= m; i++) b[i] = x[i];
1.491 +#endif
1.492 + /* (f g) := inv(P) * (b 0) */
1.493 + for (i = 1; i <= m0 + n; i++)
1.494 + fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0);
1.495 + /* f1 := inv(L0) * f */
1.496 + luf_f_solve(lpf->luf, 0, f);
1.497 + /* g1 := g - S * f1 */
1.498 + s_prod(lpf, g, -1.0, f);
1.499 + /* g2 := inv(C) * g1 */
1.500 + scf_solve_it(lpf->scf, 0, g);
1.501 + /* f2 := inv(U0) * (f1 - R * g2) */
1.502 + r_prod(lpf, f, -1.0, g);
1.503 + luf_v_solve(lpf->luf, 0, f);
1.504 + /* (x y) := inv(Q) * (f2 g2) */
1.505 + for (i = 1; i <= m; i++)
1.506 + x[i] = fg[Q_col[i]];
1.507 +#if _GLPLPF_DEBUG
1.508 + /* check relative error in solution */
1.509 + check_error(lpf, 0, x, b);
1.510 + xfree(b);
1.511 +#endif
1.512 + return;
1.513 +}
1.514 +
1.515 +/***********************************************************************
1.516 +* NAME
1.517 +*
1.518 +* lpf_btran - perform backward transformation (solve system B'*x = b)
1.519 +*
1.520 +* SYNOPSIS
1.521 +*
1.522 +* #include "glplpf.h"
1.523 +* void lpf_btran(LPF *lpf, double x[]);
1.524 +*
1.525 +* DESCRIPTION
1.526 +*
1.527 +* The routine lpf_btran performs backward transformation, i.e. solves
1.528 +* the system B'*x = b, where B' is a matrix transposed to the basis
1.529 +* matrix B, x is the vector of unknowns to be computed, b is the vector
1.530 +* of right-hand sides.
1.531 +*
1.532 +* On entry elements of the vector b should be stored in dense format
1.533 +* in locations x[1], ..., x[m], where m is the number of rows. On exit
1.534 +* the routine stores elements of the vector x in the same locations.
1.535 +*
1.536 +* BACKGROUND
1.537 +*
1.538 +* Solution of the system B' * x = b, where B' is a matrix transposed
1.539 +* to B, can be obtained by solving the following augmented system:
1.540 +*
1.541 +* ( B F^)T ( x ) ( b )
1.542 +* ( ) ( ) = ( )
1.543 +* ( G^ H^) ( y ) ( 0 )
1.544 +*
1.545 +* which, using the main equality, can be written as follows:
1.546 +*
1.547 +* T ( U0 R )T ( L0 0 )T T ( x ) ( b )
1.548 +* Q ( ) ( ) P ( ) = ( )
1.549 +* ( 0 C ) ( S I ) ( y ) ( 0 )
1.550 +*
1.551 +* or, equivalently, as follows:
1.552 +*
1.553 +* ( U'0 0 ) ( L'0 S') ( x ) ( b )
1.554 +* Q' ( ) ( ) P' ( ) = ( )
1.555 +* ( R' C') ( 0 I ) ( y ) ( 0 )
1.556 +*
1.557 +* therefore,
1.558 +*
1.559 +* ( x ) ( L'0 S')-1 ( U'0 0 )-1 ( b )
1.560 +* ( ) = P ( ) ( ) Q ( )
1.561 +* ( y ) ( 0 I ) ( R' C') ( 0 )
1.562 +*
1.563 +* Thus, computing the solution includes the following steps:
1.564 +*
1.565 +* 1. Compute
1.566 +*
1.567 +* ( f ) ( b )
1.568 +* ( ) = Q ( )
1.569 +* ( g ) ( 0 )
1.570 +*
1.571 +* 2. Solve the system
1.572 +*
1.573 +* ( f1 ) ( U'0 0 )-1 ( f ) ( U'0 0 ) ( f1 ) ( f )
1.574 +* ( ) = ( ) ( ) => ( ) ( ) = ( )
1.575 +* ( g1 ) ( R' C') ( g ) ( R' C') ( g1 ) ( g )
1.576 +*
1.577 +* from which it follows that:
1.578 +*
1.579 +* { U'0 * f1 = f f1 = inv(U'0) * f
1.580 +* { =>
1.581 +* { R' * f1 + C' * g1 = g g1 = inv(C') * (g - R' * f1)
1.582 +*
1.583 +* 3. Solve the system
1.584 +*
1.585 +* ( f2 ) ( L'0 S')-1 ( f1 ) ( L'0 S') ( f2 ) ( f1 )
1.586 +* ( ) = ( ) ( ) => ( ) ( ) = ( )
1.587 +* ( g2 ) ( 0 I ) ( g1 ) ( 0 I ) ( g2 ) ( g1 )
1.588 +*
1.589 +* from which it follows that:
1.590 +*
1.591 +* { L'0 * f2 + S' * g2 = f1
1.592 +* { => f2 = inv(L'0) * ( f1 - S' * g2)
1.593 +* { g2 = g1
1.594 +*
1.595 +* 4. Compute
1.596 +*
1.597 +* ( x ) ( f2 )
1.598 +* ( ) = P ( )
1.599 +* ( y ) ( g2 ) */
1.600 +
1.601 +void lpf_btran(LPF *lpf, double x[])
1.602 +{ int m0 = lpf->m0;
1.603 + int m = lpf->m;
1.604 + int n = lpf->n;
1.605 + int *P_row = lpf->P_row;
1.606 + int *Q_row = lpf->Q_row;
1.607 + double *fg = lpf->work1;
1.608 + double *f = fg;
1.609 + double *g = fg + m0;
1.610 + int i, ii;
1.611 +#if _GLPLPF_DEBUG
1.612 + double *b;
1.613 +#endif
1.614 + if (!lpf->valid)
1.615 + xfault("lpf_btran: the factorization is not valid\n");
1.616 + xassert(0 <= m && m <= m0 + n);
1.617 +#if _GLPLPF_DEBUG
1.618 + /* save the right-hand side vector */
1.619 + b = xcalloc(1+m, sizeof(double));
1.620 + for (i = 1; i <= m; i++) b[i] = x[i];
1.621 +#endif
1.622 + /* (f g) := Q * (b 0) */
1.623 + for (i = 1; i <= m0 + n; i++)
1.624 + fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0);
1.625 + /* f1 := inv(U'0) * f */
1.626 + luf_v_solve(lpf->luf, 1, f);
1.627 + /* g1 := inv(C') * (g - R' * f1) */
1.628 + rt_prod(lpf, g, -1.0, f);
1.629 + scf_solve_it(lpf->scf, 1, g);
1.630 + /* g2 := g1 */
1.631 + g = g;
1.632 + /* f2 := inv(L'0) * (f1 - S' * g2) */
1.633 + st_prod(lpf, f, -1.0, g);
1.634 + luf_f_solve(lpf->luf, 1, f);
1.635 + /* (x y) := P * (f2 g2) */
1.636 + for (i = 1; i <= m; i++)
1.637 + x[i] = fg[P_row[i]];
1.638 +#if _GLPLPF_DEBUG
1.639 + /* check relative error in solution */
1.640 + check_error(lpf, 1, x, b);
1.641 + xfree(b);
1.642 +#endif
1.643 + return;
1.644 +}
1.645 +
1.646 +/***********************************************************************
1.647 +* The routine enlarge_sva enlarges the Sparse Vector Area to new_size
1.648 +* locations by reallocating the arrays v_ind and v_val. */
1.649 +
1.650 +static void enlarge_sva(LPF *lpf, int new_size)
1.651 +{ int v_size = lpf->v_size;
1.652 + int used = lpf->v_ptr - 1;
1.653 + int *v_ind = lpf->v_ind;
1.654 + double *v_val = lpf->v_val;
1.655 + xassert(v_size < new_size);
1.656 + while (v_size < new_size) v_size += v_size;
1.657 + lpf->v_size = v_size;
1.658 + lpf->v_ind = xcalloc(1+v_size, sizeof(int));
1.659 + lpf->v_val = xcalloc(1+v_size, sizeof(double));
1.660 + xassert(used >= 0);
1.661 + memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int));
1.662 + memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double));
1.663 + xfree(v_ind);
1.664 + xfree(v_val);
1.665 + return;
1.666 +}
1.667 +
1.668 +/***********************************************************************
1.669 +* NAME
1.670 +*
1.671 +* lpf_update_it - update LP basis factorization
1.672 +*
1.673 +* SYNOPSIS
1.674 +*
1.675 +* #include "glplpf.h"
1.676 +* int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
1.677 +* const double val[]);
1.678 +*
1.679 +* DESCRIPTION
1.680 +*
1.681 +* The routine lpf_update_it updates the factorization of the basis
1.682 +* matrix B after replacing its j-th column by a new vector.
1.683 +*
1.684 +* The parameter j specifies the number of column of B, which has been
1.685 +* replaced, 1 <= j <= m, where m is the order of B.
1.686 +*
1.687 +* The parameter bh specifies the basis header entry for the new column
1.688 +* of B, which is the number of the new column in some original matrix.
1.689 +* This parameter is optional and can be specified as 0.
1.690 +*
1.691 +* Row indices and numerical values of non-zero elements of the new
1.692 +* column of B should be placed in locations ind[1], ..., ind[len] and
1.693 +* val[1], ..., val[len], resp., where len is the number of non-zeros
1.694 +* in the column. Neither zero nor duplicate elements are allowed.
1.695 +*
1.696 +* RETURNS
1.697 +*
1.698 +* 0 The factorization has been successfully updated.
1.699 +*
1.700 +* LPF_ESING
1.701 +* New basis B is singular within the working precision.
1.702 +*
1.703 +* LPF_ELIMIT
1.704 +* Maximal number of additional rows and columns has been reached.
1.705 +*
1.706 +* BACKGROUND
1.707 +*
1.708 +* Let j-th column of the current basis matrix B have to be replaced by
1.709 +* a new column a. This replacement is equivalent to removing the old
1.710 +* j-th column by fixing it at zero and introducing the new column as
1.711 +* follows:
1.712 +*
1.713 +* ( B F^| a )
1.714 +* ( B F^) ( | )
1.715 +* ( ) ---> ( G^ H^| 0 )
1.716 +* ( G^ H^) (-------+---)
1.717 +* ( e'j 0 | 0 )
1.718 +*
1.719 +* where ej is a unit vector with 1 in j-th position which used to fix
1.720 +* the old j-th column of B (at zero). Then using the main equality we
1.721 +* have:
1.722 +*
1.723 +* ( B F^| a ) ( B0 F | f )
1.724 +* ( | ) ( P 0 ) ( | ) ( Q 0 )
1.725 +* ( G^ H^| 0 ) = ( ) ( G H | g ) ( ) =
1.726 +* (-------+---) ( 0 1 ) (-------+---) ( 0 1 )
1.727 +* ( e'j 0 | 0 ) ( v' w'| 0 )
1.728 +*
1.729 +* [ ( B0 F )| ( f ) ] [ ( B0 F ) | ( f ) ]
1.730 +* [ P ( )| P ( ) ] ( Q 0 ) [ P ( ) Q| P ( ) ]
1.731 +* = [ ( G H )| ( g ) ] ( ) = [ ( G H ) | ( g ) ]
1.732 +* [------------+-------- ] ( 0 1 ) [-------------+---------]
1.733 +* [ ( v' w')| 0 ] [ ( v' w') Q| 0 ]
1.734 +*
1.735 +* where:
1.736 +*
1.737 +* ( a ) ( f ) ( f ) ( a )
1.738 +* ( ) = P ( ) => ( ) = P' * ( )
1.739 +* ( 0 ) ( g ) ( g ) ( 0 )
1.740 +*
1.741 +* ( ej ) ( v ) ( v ) ( ej )
1.742 +* ( e'j 0 ) = ( v' w' ) Q => ( ) = Q' ( ) => ( ) = Q ( )
1.743 +* ( 0 ) ( w ) ( w ) ( 0 )
1.744 +*
1.745 +* On the other hand:
1.746 +*
1.747 +* ( B0| F f )
1.748 +* ( P 0 ) (---+------) ( Q 0 ) ( B0 new F )
1.749 +* ( ) ( G | H g ) ( ) = new P ( ) new Q
1.750 +* ( 0 1 ) ( | ) ( 0 1 ) ( new G new H )
1.751 +* ( v'| w' 0 )
1.752 +*
1.753 +* where:
1.754 +* ( G ) ( H g )
1.755 +* new F = ( F f ), new G = ( ), new H = ( ),
1.756 +* ( v') ( w' 0 )
1.757 +*
1.758 +* ( P 0 ) ( Q 0 )
1.759 +* new P = ( ) , new Q = ( ) .
1.760 +* ( 0 1 ) ( 0 1 )
1.761 +*
1.762 +* The factorization structure for the new augmented matrix remains the
1.763 +* same, therefore:
1.764 +*
1.765 +* ( B0 new F ) ( L0 0 ) ( U0 new R )
1.766 +* new P ( ) new Q = ( ) ( )
1.767 +* ( new G new H ) ( new S I ) ( 0 new C )
1.768 +*
1.769 +* where:
1.770 +*
1.771 +* new F = L0 * new R =>
1.772 +*
1.773 +* new R = inv(L0) * new F = inv(L0) * (F f) = ( R inv(L0)*f )
1.774 +*
1.775 +* new G = new S * U0 =>
1.776 +*
1.777 +* ( G ) ( S )
1.778 +* new S = new G * inv(U0) = ( ) * inv(U0) = ( )
1.779 +* ( v') ( v'*inv(U0) )
1.780 +*
1.781 +* new H = new S * new R + new C =>
1.782 +*
1.783 +* new C = new H - new S * new R =
1.784 +*
1.785 +* ( H g ) ( S )
1.786 +* = ( ) - ( ) * ( R inv(L0)*f ) =
1.787 +* ( w' 0 ) ( v'*inv(U0) )
1.788 +*
1.789 +* ( H - S*R g - S*inv(L0)*f ) ( C x )
1.790 +* = ( ) = ( )
1.791 +* ( w'- v'*inv(U0)*R -v'*inv(U0)*inv(L0)*f) ( y' z )
1.792 +*
1.793 +* Note that new C is resulted by expanding old C with new column x,
1.794 +* row y', and diagonal element z, where:
1.795 +*
1.796 +* x = g - S * inv(L0) * f = g - S * (new column of R)
1.797 +*
1.798 +* y = w - R'* inv(U'0)* v = w - R'* (new row of S)
1.799 +*
1.800 +* z = - (new row of S) * (new column of R)
1.801 +*
1.802 +* Finally, to replace old B by new B we have to permute j-th and last
1.803 +* (just added) columns of the matrix
1.804 +*
1.805 +* ( B F^| a )
1.806 +* ( | )
1.807 +* ( G^ H^| 0 )
1.808 +* (-------+---)
1.809 +* ( e'j 0 | 0 )
1.810 +*
1.811 +* and to keep the main equality do the same for matrix Q. */
1.812 +
1.813 +int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[],
1.814 + const double val[])
1.815 +{ int m0 = lpf->m0;
1.816 + int m = lpf->m;
1.817 +#if _GLPLPF_DEBUG
1.818 + double *B = lpf->B;
1.819 +#endif
1.820 + int n = lpf->n;
1.821 + int *R_ptr = lpf->R_ptr;
1.822 + int *R_len = lpf->R_len;
1.823 + int *S_ptr = lpf->S_ptr;
1.824 + int *S_len = lpf->S_len;
1.825 + int *P_row = lpf->P_row;
1.826 + int *P_col = lpf->P_col;
1.827 + int *Q_row = lpf->Q_row;
1.828 + int *Q_col = lpf->Q_col;
1.829 + int v_ptr = lpf->v_ptr;
1.830 + int *v_ind = lpf->v_ind;
1.831 + double *v_val = lpf->v_val;
1.832 + double *a = lpf->work2; /* new column */
1.833 + double *fg = lpf->work1, *f = fg, *g = fg + m0;
1.834 + double *vw = lpf->work2, *v = vw, *w = vw + m0;
1.835 + double *x = g, *y = w, z;
1.836 + int i, ii, k, ret;
1.837 + xassert(bh == bh);
1.838 + if (!lpf->valid)
1.839 + xfault("lpf_update_it: the factorization is not valid\n");
1.840 + if (!(1 <= j && j <= m))
1.841 + xfault("lpf_update_it: j = %d; column number out of range\n",
1.842 + j);
1.843 + xassert(0 <= m && m <= m0 + n);
1.844 + /* check if the basis factorization can be expanded */
1.845 + if (n == lpf->n_max)
1.846 + { lpf->valid = 0;
1.847 + ret = LPF_ELIMIT;
1.848 + goto done;
1.849 + }
1.850 + /* convert new j-th column of B to dense format */
1.851 + for (i = 1; i <= m; i++)
1.852 + a[i] = 0.0;
1.853 + for (k = 1; k <= len; k++)
1.854 + { i = ind[k];
1.855 + if (!(1 <= i && i <= m))
1.856 + xfault("lpf_update_it: ind[%d] = %d; row number out of rang"
1.857 + "e\n", k, i);
1.858 + if (a[i] != 0.0)
1.859 + xfault("lpf_update_it: ind[%d] = %d; duplicate row index no"
1.860 + "t allowed\n", k, i);
1.861 + if (val[k] == 0.0)
1.862 + xfault("lpf_update_it: val[%d] = %g; zero element not allow"
1.863 + "ed\n", k, val[k]);
1.864 + a[i] = val[k];
1.865 + }
1.866 +#if _GLPLPF_DEBUG
1.867 + /* change column in the basis matrix for debugging */
1.868 + for (i = 1; i <= m; i++)
1.869 + B[(i - 1) * m + j] = a[i];
1.870 +#endif
1.871 + /* (f g) := inv(P) * (a 0) */
1.872 + for (i = 1; i <= m0+n; i++)
1.873 + fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0);
1.874 + /* (v w) := Q * (ej 0) */
1.875 + for (i = 1; i <= m0+n; i++) vw[i] = 0.0;
1.876 + vw[Q_col[j]] = 1.0;
1.877 + /* f1 := inv(L0) * f (new column of R) */
1.878 + luf_f_solve(lpf->luf, 0, f);
1.879 + /* v1 := inv(U'0) * v (new row of S) */
1.880 + luf_v_solve(lpf->luf, 1, v);
1.881 + /* we need at most 2 * m0 available locations in the SVA to store
1.882 + new column of matrix R and new row of matrix S */
1.883 + if (lpf->v_size < v_ptr + m0 + m0)
1.884 + { enlarge_sva(lpf, v_ptr + m0 + m0);
1.885 + v_ind = lpf->v_ind;
1.886 + v_val = lpf->v_val;
1.887 + }
1.888 + /* store new column of R */
1.889 + R_ptr[n+1] = v_ptr;
1.890 + for (i = 1; i <= m0; i++)
1.891 + { if (f[i] != 0.0)
1.892 + v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++;
1.893 + }
1.894 + R_len[n+1] = v_ptr - lpf->v_ptr;
1.895 + lpf->v_ptr = v_ptr;
1.896 + /* store new row of S */
1.897 + S_ptr[n+1] = v_ptr;
1.898 + for (i = 1; i <= m0; i++)
1.899 + { if (v[i] != 0.0)
1.900 + v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++;
1.901 + }
1.902 + S_len[n+1] = v_ptr - lpf->v_ptr;
1.903 + lpf->v_ptr = v_ptr;
1.904 + /* x := g - S * f1 (new column of C) */
1.905 + s_prod(lpf, x, -1.0, f);
1.906 + /* y := w - R' * v1 (new row of C) */
1.907 + rt_prod(lpf, y, -1.0, v);
1.908 + /* z := - v1 * f1 (new diagonal element of C) */
1.909 + z = 0.0;
1.910 + for (i = 1; i <= m0; i++) z -= v[i] * f[i];
1.911 + /* update factorization of new matrix C */
1.912 + switch (scf_update_exp(lpf->scf, x, y, z))
1.913 + { case 0:
1.914 + break;
1.915 + case SCF_ESING:
1.916 + lpf->valid = 0;
1.917 + ret = LPF_ESING;
1.918 + goto done;
1.919 + case SCF_ELIMIT:
1.920 + xassert(lpf != lpf);
1.921 + default:
1.922 + xassert(lpf != lpf);
1.923 + }
1.924 + /* expand matrix P */
1.925 + P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1;
1.926 + /* expand matrix Q */
1.927 + Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1;
1.928 + /* permute j-th and last (just added) column of matrix Q */
1.929 + i = Q_col[j], ii = Q_col[m0+n+1];
1.930 + Q_row[i] = m0+n+1, Q_col[m0+n+1] = i;
1.931 + Q_row[ii] = j, Q_col[j] = ii;
1.932 + /* increase the number of additional rows and columns */
1.933 + lpf->n++;
1.934 + xassert(lpf->n <= lpf->n_max);
1.935 + /* the factorization has been successfully updated */
1.936 + ret = 0;
1.937 +done: /* return to the calling program */
1.938 + return ret;
1.939 +}
1.940 +
1.941 +/***********************************************************************
1.942 +* NAME
1.943 +*
1.944 +* lpf_delete_it - delete LP basis factorization
1.945 +*
1.946 +* SYNOPSIS
1.947 +*
1.948 +* #include "glplpf.h"
1.949 +* void lpf_delete_it(LPF *lpf)
1.950 +*
1.951 +* DESCRIPTION
1.952 +*
1.953 +* The routine lpf_delete_it deletes LP basis factorization specified
1.954 +* by the parameter lpf and frees all memory allocated to this program
1.955 +* object. */
1.956 +
1.957 +void lpf_delete_it(LPF *lpf)
1.958 +{ luf_delete_it(lpf->luf);
1.959 +#if _GLPLPF_DEBUG
1.960 + if (lpf->B != NULL) xfree(lpf->B);
1.961 +#else
1.962 + xassert(lpf->B == NULL);
1.963 +#endif
1.964 + if (lpf->R_ptr != NULL) xfree(lpf->R_ptr);
1.965 + if (lpf->R_len != NULL) xfree(lpf->R_len);
1.966 + if (lpf->S_ptr != NULL) xfree(lpf->S_ptr);
1.967 + if (lpf->S_len != NULL) xfree(lpf->S_len);
1.968 + if (lpf->scf != NULL) scf_delete_it(lpf->scf);
1.969 + if (lpf->P_row != NULL) xfree(lpf->P_row);
1.970 + if (lpf->P_col != NULL) xfree(lpf->P_col);
1.971 + if (lpf->Q_row != NULL) xfree(lpf->Q_row);
1.972 + if (lpf->Q_col != NULL) xfree(lpf->Q_col);
1.973 + if (lpf->v_ind != NULL) xfree(lpf->v_ind);
1.974 + if (lpf->v_val != NULL) xfree(lpf->v_val);
1.975 + if (lpf->work1 != NULL) xfree(lpf->work1);
1.976 + if (lpf->work2 != NULL) xfree(lpf->work2);
1.977 + xfree(lpf);
1.978 + return;
1.979 +}
1.980 +
1.981 +/* eof */