diff -r d59bea55db9b -r c445c931472f src/glplpf.c --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/glplpf.c Mon Dec 06 13:09:21 2010 +0100 @@ -0,0 +1,978 @@ +/* glplpf.c (LP basis factorization, Schur complement version) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics, +* Moscow Aviation Institute, Moscow, Russia. All rights reserved. +* E-mail: . +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see . +***********************************************************************/ + +#include "glplpf.h" +#include "glpenv.h" +#define xfault xerror + +#define _GLPLPF_DEBUG 0 + +/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ + +#define M_MAX 100000000 /* = 100*10^6 */ +/* maximal order of the basis matrix */ + +/*********************************************************************** +* NAME +* +* lpf_create_it - create LP basis factorization +* +* SYNOPSIS +* +* #include "glplpf.h" +* LPF *lpf_create_it(void); +* +* DESCRIPTION +* +* The routine lpf_create_it creates a program object, which represents +* a factorization of LP basis. +* +* RETURNS +* +* The routine lpf_create_it returns a pointer to the object created. */ + +LPF *lpf_create_it(void) +{ LPF *lpf; +#if _GLPLPF_DEBUG + xprintf("lpf_create_it: warning: debug mode enabled\n"); +#endif + lpf = xmalloc(sizeof(LPF)); + lpf->valid = 0; + lpf->m0_max = lpf->m0 = 0; + lpf->luf = luf_create_it(); + lpf->m = 0; + lpf->B = NULL; + lpf->n_max = 50; + lpf->n = 0; + lpf->R_ptr = lpf->R_len = NULL; + lpf->S_ptr = lpf->S_len = NULL; + lpf->scf = NULL; + lpf->P_row = lpf->P_col = NULL; + lpf->Q_row = lpf->Q_col = NULL; + lpf->v_size = 1000; + lpf->v_ptr = 0; + lpf->v_ind = NULL; + lpf->v_val = NULL; + lpf->work1 = lpf->work2 = NULL; + return lpf; +} + +/*********************************************************************** +* NAME +* +* lpf_factorize - compute LP basis factorization +* +* SYNOPSIS +* +* #include "glplpf.h" +* int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col) +* (void *info, int j, int ind[], double val[]), void *info); +* +* DESCRIPTION +* +* The routine lpf_factorize computes the factorization of the basis +* matrix B specified by the routine col. +* +* The parameter lpf specified the basis factorization data structure +* created with the routine lpf_create_it. +* +* The parameter m specifies the order of B, m > 0. +* +* The array bh specifies the basis header: bh[j], 1 <= j <= m, is the +* number of j-th column of B in some original matrix. The array bh is +* optional and can be specified as NULL. +* +* The formal routine col specifies the matrix B to be factorized. To +* obtain j-th column of A the routine lpf_factorize calls the routine +* col with the parameter j (1 <= j <= n). In response the routine col +* should store row indices and numerical values of non-zero elements +* of j-th column of B to locations ind[1,...,len] and val[1,...,len], +* respectively, where len is the number of non-zeros in j-th column +* returned on exit. Neither zero nor duplicate elements are allowed. +* +* The parameter info is a transit pointer passed to the routine col. +* +* RETURNS +* +* 0 The factorization has been successfully computed. +* +* LPF_ESING +* The specified matrix is singular within the working precision. +* +* LPF_ECOND +* The specified matrix is ill-conditioned. +* +* For more details see comments to the routine luf_factorize. */ + +int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col) + (void *info, int j, int ind[], double val[]), void *info) +{ int k, ret; +#if _GLPLPF_DEBUG + int i, j, len, *ind; + double *B, *val; +#endif + xassert(bh == bh); + if (m < 1) + xfault("lpf_factorize: m = %d; invalid parameter\n", m); + if (m > M_MAX) + xfault("lpf_factorize: m = %d; matrix too big\n", m); + lpf->m0 = lpf->m = m; + /* invalidate the factorization */ + lpf->valid = 0; + /* allocate/reallocate arrays, if necessary */ + if (lpf->R_ptr == NULL) + lpf->R_ptr = xcalloc(1+lpf->n_max, sizeof(int)); + if (lpf->R_len == NULL) + lpf->R_len = xcalloc(1+lpf->n_max, sizeof(int)); + if (lpf->S_ptr == NULL) + lpf->S_ptr = xcalloc(1+lpf->n_max, sizeof(int)); + if (lpf->S_len == NULL) + lpf->S_len = xcalloc(1+lpf->n_max, sizeof(int)); + if (lpf->scf == NULL) + lpf->scf = scf_create_it(lpf->n_max); + if (lpf->v_ind == NULL) + lpf->v_ind = xcalloc(1+lpf->v_size, sizeof(int)); + if (lpf->v_val == NULL) + lpf->v_val = xcalloc(1+lpf->v_size, sizeof(double)); + if (lpf->m0_max < m) + { if (lpf->P_row != NULL) xfree(lpf->P_row); + if (lpf->P_col != NULL) xfree(lpf->P_col); + if (lpf->Q_row != NULL) xfree(lpf->Q_row); + if (lpf->Q_col != NULL) xfree(lpf->Q_col); + if (lpf->work1 != NULL) xfree(lpf->work1); + if (lpf->work2 != NULL) xfree(lpf->work2); + lpf->m0_max = m + 100; + lpf->P_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); + lpf->P_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); + lpf->Q_row = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); + lpf->Q_col = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(int)); + lpf->work1 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double)); + lpf->work2 = xcalloc(1+lpf->m0_max+lpf->n_max, sizeof(double)); + } + /* try to factorize the basis matrix */ + switch (luf_factorize(lpf->luf, m, col, info)) + { case 0: + break; + case LUF_ESING: + ret = LPF_ESING; + goto done; + case LUF_ECOND: + ret = LPF_ECOND; + goto done; + default: + xassert(lpf != lpf); + } + /* the basis matrix has been successfully factorized */ + lpf->valid = 1; +#if _GLPLPF_DEBUG + /* store the basis matrix for debugging */ + if (lpf->B != NULL) xfree(lpf->B); + xassert(m <= 32767); + lpf->B = B = xcalloc(1+m*m, sizeof(double)); + ind = xcalloc(1+m, sizeof(int)); + val = xcalloc(1+m, sizeof(double)); + for (k = 1; k <= m * m; k++) + B[k] = 0.0; + for (j = 1; j <= m; j++) + { len = col(info, j, ind, val); + xassert(0 <= len && len <= m); + for (k = 1; k <= len; k++) + { i = ind[k]; + xassert(1 <= i && i <= m); + xassert(B[(i - 1) * m + j] == 0.0); + xassert(val[k] != 0.0); + B[(i - 1) * m + j] = val[k]; + } + } + xfree(ind); + xfree(val); +#endif + /* B = B0, so there are no additional rows/columns */ + lpf->n = 0; + /* reset the Schur complement factorization */ + scf_reset_it(lpf->scf); + /* P := Q := I */ + for (k = 1; k <= m; k++) + { lpf->P_row[k] = lpf->P_col[k] = k; + lpf->Q_row[k] = lpf->Q_col[k] = k; + } + /* make all SVA locations free */ + lpf->v_ptr = 1; + ret = 0; +done: /* return to the calling program */ + return ret; +} + +/*********************************************************************** +* The routine r_prod computes the product y := y + alpha * R * x, +* where x is a n-vector, alpha is a scalar, y is a m0-vector. +* +* Since matrix R is available by columns, the product is computed as +* a linear combination: +* +* y := y + alpha * (R[1] * x[1] + ... + R[n] * x[n]), +* +* where R[j] is j-th column of R. */ + +static void r_prod(LPF *lpf, double y[], double a, const double x[]) +{ int n = lpf->n; + int *R_ptr = lpf->R_ptr; + int *R_len = lpf->R_len; + int *v_ind = lpf->v_ind; + double *v_val = lpf->v_val; + int j, beg, end, ptr; + double t; + for (j = 1; j <= n; j++) + { if (x[j] == 0.0) continue; + /* y := y + alpha * R[j] * x[j] */ + t = a * x[j]; + beg = R_ptr[j]; + end = beg + R_len[j]; + for (ptr = beg; ptr < end; ptr++) + y[v_ind[ptr]] += t * v_val[ptr]; + } + return; +} + +/*********************************************************************** +* The routine rt_prod computes the product y := y + alpha * R' * x, +* where R' is a matrix transposed to R, x is a m0-vector, alpha is a +* scalar, y is a n-vector. +* +* Since matrix R is available by columns, the product components are +* computed as inner products: +* +* y[j] := y[j] + alpha * (j-th column of R) * x +* +* for j = 1, 2, ..., n. */ + +static void rt_prod(LPF *lpf, double y[], double a, const double x[]) +{ int n = lpf->n; + int *R_ptr = lpf->R_ptr; + int *R_len = lpf->R_len; + int *v_ind = lpf->v_ind; + double *v_val = lpf->v_val; + int j, beg, end, ptr; + double t; + for (j = 1; j <= n; j++) + { /* t := (j-th column of R) * x */ + t = 0.0; + beg = R_ptr[j]; + end = beg + R_len[j]; + for (ptr = beg; ptr < end; ptr++) + t += v_val[ptr] * x[v_ind[ptr]]; + /* y[j] := y[j] + alpha * t */ + y[j] += a * t; + } + return; +} + +/*********************************************************************** +* The routine s_prod computes the product y := y + alpha * S * x, +* where x is a m0-vector, alpha is a scalar, y is a n-vector. +* +* Since matrix S is available by rows, the product components are +* computed as inner products: +* +* y[i] = y[i] + alpha * (i-th row of S) * x +* +* for i = 1, 2, ..., n. */ + +static void s_prod(LPF *lpf, double y[], double a, const double x[]) +{ int n = lpf->n; + int *S_ptr = lpf->S_ptr; + int *S_len = lpf->S_len; + int *v_ind = lpf->v_ind; + double *v_val = lpf->v_val; + int i, beg, end, ptr; + double t; + for (i = 1; i <= n; i++) + { /* t := (i-th row of S) * x */ + t = 0.0; + beg = S_ptr[i]; + end = beg + S_len[i]; + for (ptr = beg; ptr < end; ptr++) + t += v_val[ptr] * x[v_ind[ptr]]; + /* y[i] := y[i] + alpha * t */ + y[i] += a * t; + } + return; +} + +/*********************************************************************** +* The routine st_prod computes the product y := y + alpha * S' * x, +* where S' is a matrix transposed to S, x is a n-vector, alpha is a +* scalar, y is m0-vector. +* +* Since matrix R is available by rows, the product is computed as a +* linear combination: +* +* y := y + alpha * (S'[1] * x[1] + ... + S'[n] * x[n]), +* +* where S'[i] is i-th row of S. */ + +static void st_prod(LPF *lpf, double y[], double a, const double x[]) +{ int n = lpf->n; + int *S_ptr = lpf->S_ptr; + int *S_len = lpf->S_len; + int *v_ind = lpf->v_ind; + double *v_val = lpf->v_val; + int i, beg, end, ptr; + double t; + for (i = 1; i <= n; i++) + { if (x[i] == 0.0) continue; + /* y := y + alpha * S'[i] * x[i] */ + t = a * x[i]; + beg = S_ptr[i]; + end = beg + S_len[i]; + for (ptr = beg; ptr < end; ptr++) + y[v_ind[ptr]] += t * v_val[ptr]; + } + return; +} + +#if _GLPLPF_DEBUG +/*********************************************************************** +* The routine check_error computes the maximal relative error between +* left- and right-hand sides for the system B * x = b (if tr is zero) +* or B' * x = b (if tr is non-zero), where B' is a matrix transposed +* to B. (This routine is intended for debugging only.) */ + +static void check_error(LPF *lpf, int tr, const double x[], + const double b[]) +{ int m = lpf->m; + double *B = lpf->B; + int i, j; + double d, dmax = 0.0, s, t, tmax; + for (i = 1; i <= m; i++) + { s = 0.0; + tmax = 1.0; + for (j = 1; j <= m; j++) + { if (!tr) + t = B[m * (i - 1) + j] * x[j]; + else + t = B[m * (j - 1) + i] * x[j]; + if (tmax < fabs(t)) tmax = fabs(t); + s += t; + } + d = fabs(s - b[i]) / tmax; + if (dmax < d) dmax = d; + } + if (dmax > 1e-8) + xprintf("%s: dmax = %g; relative error too large\n", + !tr ? "lpf_ftran" : "lpf_btran", dmax); + return; +} +#endif + +/*********************************************************************** +* NAME +* +* lpf_ftran - perform forward transformation (solve system B*x = b) +* +* SYNOPSIS +* +* #include "glplpf.h" +* void lpf_ftran(LPF *lpf, double x[]); +* +* DESCRIPTION +* +* The routine lpf_ftran performs forward transformation, i.e. solves +* the system B*x = b, where B is the basis matrix, x is the vector of +* unknowns to be computed, b is the vector of right-hand sides. +* +* On entry elements of the vector b should be stored in dense format +* in locations x[1], ..., x[m], where m is the number of rows. On exit +* the routine stores elements of the vector x in the same locations. +* +* BACKGROUND +* +* Solution of the system B * x = b can be obtained by solving the +* following augmented system: +* +* ( B F^) ( x ) ( b ) +* ( ) ( ) = ( ) +* ( G^ H^) ( y ) ( 0 ) +* +* which, using the main equality, can be written as follows: +* +* ( L0 0 ) ( U0 R ) ( x ) ( b ) +* P ( ) ( ) Q ( ) = ( ) +* ( S I ) ( 0 C ) ( y ) ( 0 ) +* +* therefore, +* +* ( x ) ( U0 R )-1 ( L0 0 )-1 ( b ) +* ( ) = Q' ( ) ( ) P' ( ) +* ( y ) ( 0 C ) ( S I ) ( 0 ) +* +* Thus, computing the solution includes the following steps: +* +* 1. Compute +* +* ( f ) ( b ) +* ( ) = P' ( ) +* ( g ) ( 0 ) +* +* 2. Solve the system +* +* ( f1 ) ( L0 0 )-1 ( f ) ( L0 0 ) ( f1 ) ( f ) +* ( ) = ( ) ( ) => ( ) ( ) = ( ) +* ( g1 ) ( S I ) ( g ) ( S I ) ( g1 ) ( g ) +* +* from which it follows that: +* +* { L0 * f1 = f f1 = inv(L0) * f +* { => +* { S * f1 + g1 = g g1 = g - S * f1 +* +* 3. Solve the system +* +* ( f2 ) ( U0 R )-1 ( f1 ) ( U0 R ) ( f2 ) ( f1 ) +* ( ) = ( ) ( ) => ( ) ( ) = ( ) +* ( g2 ) ( 0 C ) ( g1 ) ( 0 C ) ( g2 ) ( g1 ) +* +* from which it follows that: +* +* { U0 * f2 + R * g2 = f1 f2 = inv(U0) * (f1 - R * g2) +* { => +* { C * g2 = g1 g2 = inv(C) * g1 +* +* 4. Compute +* +* ( x ) ( f2 ) +* ( ) = Q' ( ) +* ( y ) ( g2 ) */ + +void lpf_ftran(LPF *lpf, double x[]) +{ int m0 = lpf->m0; + int m = lpf->m; + int n = lpf->n; + int *P_col = lpf->P_col; + int *Q_col = lpf->Q_col; + double *fg = lpf->work1; + double *f = fg; + double *g = fg + m0; + int i, ii; +#if _GLPLPF_DEBUG + double *b; +#endif + if (!lpf->valid) + xfault("lpf_ftran: the factorization is not valid\n"); + xassert(0 <= m && m <= m0 + n); +#if _GLPLPF_DEBUG + /* save the right-hand side vector */ + b = xcalloc(1+m, sizeof(double)); + for (i = 1; i <= m; i++) b[i] = x[i]; +#endif + /* (f g) := inv(P) * (b 0) */ + for (i = 1; i <= m0 + n; i++) + fg[i] = ((ii = P_col[i]) <= m ? x[ii] : 0.0); + /* f1 := inv(L0) * f */ + luf_f_solve(lpf->luf, 0, f); + /* g1 := g - S * f1 */ + s_prod(lpf, g, -1.0, f); + /* g2 := inv(C) * g1 */ + scf_solve_it(lpf->scf, 0, g); + /* f2 := inv(U0) * (f1 - R * g2) */ + r_prod(lpf, f, -1.0, g); + luf_v_solve(lpf->luf, 0, f); + /* (x y) := inv(Q) * (f2 g2) */ + for (i = 1; i <= m; i++) + x[i] = fg[Q_col[i]]; +#if _GLPLPF_DEBUG + /* check relative error in solution */ + check_error(lpf, 0, x, b); + xfree(b); +#endif + return; +} + +/*********************************************************************** +* NAME +* +* lpf_btran - perform backward transformation (solve system B'*x = b) +* +* SYNOPSIS +* +* #include "glplpf.h" +* void lpf_btran(LPF *lpf, double x[]); +* +* DESCRIPTION +* +* The routine lpf_btran performs backward transformation, i.e. solves +* the system B'*x = b, where B' is a matrix transposed to the basis +* matrix B, x is the vector of unknowns to be computed, b is the vector +* of right-hand sides. +* +* On entry elements of the vector b should be stored in dense format +* in locations x[1], ..., x[m], where m is the number of rows. On exit +* the routine stores elements of the vector x in the same locations. +* +* BACKGROUND +* +* Solution of the system B' * x = b, where B' is a matrix transposed +* to B, can be obtained by solving the following augmented system: +* +* ( B F^)T ( x ) ( b ) +* ( ) ( ) = ( ) +* ( G^ H^) ( y ) ( 0 ) +* +* which, using the main equality, can be written as follows: +* +* T ( U0 R )T ( L0 0 )T T ( x ) ( b ) +* Q ( ) ( ) P ( ) = ( ) +* ( 0 C ) ( S I ) ( y ) ( 0 ) +* +* or, equivalently, as follows: +* +* ( U'0 0 ) ( L'0 S') ( x ) ( b ) +* Q' ( ) ( ) P' ( ) = ( ) +* ( R' C') ( 0 I ) ( y ) ( 0 ) +* +* therefore, +* +* ( x ) ( L'0 S')-1 ( U'0 0 )-1 ( b ) +* ( ) = P ( ) ( ) Q ( ) +* ( y ) ( 0 I ) ( R' C') ( 0 ) +* +* Thus, computing the solution includes the following steps: +* +* 1. Compute +* +* ( f ) ( b ) +* ( ) = Q ( ) +* ( g ) ( 0 ) +* +* 2. Solve the system +* +* ( f1 ) ( U'0 0 )-1 ( f ) ( U'0 0 ) ( f1 ) ( f ) +* ( ) = ( ) ( ) => ( ) ( ) = ( ) +* ( g1 ) ( R' C') ( g ) ( R' C') ( g1 ) ( g ) +* +* from which it follows that: +* +* { U'0 * f1 = f f1 = inv(U'0) * f +* { => +* { R' * f1 + C' * g1 = g g1 = inv(C') * (g - R' * f1) +* +* 3. Solve the system +* +* ( f2 ) ( L'0 S')-1 ( f1 ) ( L'0 S') ( f2 ) ( f1 ) +* ( ) = ( ) ( ) => ( ) ( ) = ( ) +* ( g2 ) ( 0 I ) ( g1 ) ( 0 I ) ( g2 ) ( g1 ) +* +* from which it follows that: +* +* { L'0 * f2 + S' * g2 = f1 +* { => f2 = inv(L'0) * ( f1 - S' * g2) +* { g2 = g1 +* +* 4. Compute +* +* ( x ) ( f2 ) +* ( ) = P ( ) +* ( y ) ( g2 ) */ + +void lpf_btran(LPF *lpf, double x[]) +{ int m0 = lpf->m0; + int m = lpf->m; + int n = lpf->n; + int *P_row = lpf->P_row; + int *Q_row = lpf->Q_row; + double *fg = lpf->work1; + double *f = fg; + double *g = fg + m0; + int i, ii; +#if _GLPLPF_DEBUG + double *b; +#endif + if (!lpf->valid) + xfault("lpf_btran: the factorization is not valid\n"); + xassert(0 <= m && m <= m0 + n); +#if _GLPLPF_DEBUG + /* save the right-hand side vector */ + b = xcalloc(1+m, sizeof(double)); + for (i = 1; i <= m; i++) b[i] = x[i]; +#endif + /* (f g) := Q * (b 0) */ + for (i = 1; i <= m0 + n; i++) + fg[i] = ((ii = Q_row[i]) <= m ? x[ii] : 0.0); + /* f1 := inv(U'0) * f */ + luf_v_solve(lpf->luf, 1, f); + /* g1 := inv(C') * (g - R' * f1) */ + rt_prod(lpf, g, -1.0, f); + scf_solve_it(lpf->scf, 1, g); + /* g2 := g1 */ + g = g; + /* f2 := inv(L'0) * (f1 - S' * g2) */ + st_prod(lpf, f, -1.0, g); + luf_f_solve(lpf->luf, 1, f); + /* (x y) := P * (f2 g2) */ + for (i = 1; i <= m; i++) + x[i] = fg[P_row[i]]; +#if _GLPLPF_DEBUG + /* check relative error in solution */ + check_error(lpf, 1, x, b); + xfree(b); +#endif + return; +} + +/*********************************************************************** +* The routine enlarge_sva enlarges the Sparse Vector Area to new_size +* locations by reallocating the arrays v_ind and v_val. */ + +static void enlarge_sva(LPF *lpf, int new_size) +{ int v_size = lpf->v_size; + int used = lpf->v_ptr - 1; + int *v_ind = lpf->v_ind; + double *v_val = lpf->v_val; + xassert(v_size < new_size); + while (v_size < new_size) v_size += v_size; + lpf->v_size = v_size; + lpf->v_ind = xcalloc(1+v_size, sizeof(int)); + lpf->v_val = xcalloc(1+v_size, sizeof(double)); + xassert(used >= 0); + memcpy(&lpf->v_ind[1], &v_ind[1], used * sizeof(int)); + memcpy(&lpf->v_val[1], &v_val[1], used * sizeof(double)); + xfree(v_ind); + xfree(v_val); + return; +} + +/*********************************************************************** +* NAME +* +* lpf_update_it - update LP basis factorization +* +* SYNOPSIS +* +* #include "glplpf.h" +* int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[], +* const double val[]); +* +* DESCRIPTION +* +* The routine lpf_update_it updates the factorization of the basis +* matrix B after replacing its j-th column by a new vector. +* +* The parameter j specifies the number of column of B, which has been +* replaced, 1 <= j <= m, where m is the order of B. +* +* The parameter bh specifies the basis header entry for the new column +* of B, which is the number of the new column in some original matrix. +* This parameter is optional and can be specified as 0. +* +* Row indices and numerical values of non-zero elements of the new +* column of B should be placed in locations ind[1], ..., ind[len] and +* val[1], ..., val[len], resp., where len is the number of non-zeros +* in the column. Neither zero nor duplicate elements are allowed. +* +* RETURNS +* +* 0 The factorization has been successfully updated. +* +* LPF_ESING +* New basis B is singular within the working precision. +* +* LPF_ELIMIT +* Maximal number of additional rows and columns has been reached. +* +* BACKGROUND +* +* Let j-th column of the current basis matrix B have to be replaced by +* a new column a. This replacement is equivalent to removing the old +* j-th column by fixing it at zero and introducing the new column as +* follows: +* +* ( B F^| a ) +* ( B F^) ( | ) +* ( ) ---> ( G^ H^| 0 ) +* ( G^ H^) (-------+---) +* ( e'j 0 | 0 ) +* +* where ej is a unit vector with 1 in j-th position which used to fix +* the old j-th column of B (at zero). Then using the main equality we +* have: +* +* ( B F^| a ) ( B0 F | f ) +* ( | ) ( P 0 ) ( | ) ( Q 0 ) +* ( G^ H^| 0 ) = ( ) ( G H | g ) ( ) = +* (-------+---) ( 0 1 ) (-------+---) ( 0 1 ) +* ( e'j 0 | 0 ) ( v' w'| 0 ) +* +* [ ( B0 F )| ( f ) ] [ ( B0 F ) | ( f ) ] +* [ P ( )| P ( ) ] ( Q 0 ) [ P ( ) Q| P ( ) ] +* = [ ( G H )| ( g ) ] ( ) = [ ( G H ) | ( g ) ] +* [------------+-------- ] ( 0 1 ) [-------------+---------] +* [ ( v' w')| 0 ] [ ( v' w') Q| 0 ] +* +* where: +* +* ( a ) ( f ) ( f ) ( a ) +* ( ) = P ( ) => ( ) = P' * ( ) +* ( 0 ) ( g ) ( g ) ( 0 ) +* +* ( ej ) ( v ) ( v ) ( ej ) +* ( e'j 0 ) = ( v' w' ) Q => ( ) = Q' ( ) => ( ) = Q ( ) +* ( 0 ) ( w ) ( w ) ( 0 ) +* +* On the other hand: +* +* ( B0| F f ) +* ( P 0 ) (---+------) ( Q 0 ) ( B0 new F ) +* ( ) ( G | H g ) ( ) = new P ( ) new Q +* ( 0 1 ) ( | ) ( 0 1 ) ( new G new H ) +* ( v'| w' 0 ) +* +* where: +* ( G ) ( H g ) +* new F = ( F f ), new G = ( ), new H = ( ), +* ( v') ( w' 0 ) +* +* ( P 0 ) ( Q 0 ) +* new P = ( ) , new Q = ( ) . +* ( 0 1 ) ( 0 1 ) +* +* The factorization structure for the new augmented matrix remains the +* same, therefore: +* +* ( B0 new F ) ( L0 0 ) ( U0 new R ) +* new P ( ) new Q = ( ) ( ) +* ( new G new H ) ( new S I ) ( 0 new C ) +* +* where: +* +* new F = L0 * new R => +* +* new R = inv(L0) * new F = inv(L0) * (F f) = ( R inv(L0)*f ) +* +* new G = new S * U0 => +* +* ( G ) ( S ) +* new S = new G * inv(U0) = ( ) * inv(U0) = ( ) +* ( v') ( v'*inv(U0) ) +* +* new H = new S * new R + new C => +* +* new C = new H - new S * new R = +* +* ( H g ) ( S ) +* = ( ) - ( ) * ( R inv(L0)*f ) = +* ( w' 0 ) ( v'*inv(U0) ) +* +* ( H - S*R g - S*inv(L0)*f ) ( C x ) +* = ( ) = ( ) +* ( w'- v'*inv(U0)*R -v'*inv(U0)*inv(L0)*f) ( y' z ) +* +* Note that new C is resulted by expanding old C with new column x, +* row y', and diagonal element z, where: +* +* x = g - S * inv(L0) * f = g - S * (new column of R) +* +* y = w - R'* inv(U'0)* v = w - R'* (new row of S) +* +* z = - (new row of S) * (new column of R) +* +* Finally, to replace old B by new B we have to permute j-th and last +* (just added) columns of the matrix +* +* ( B F^| a ) +* ( | ) +* ( G^ H^| 0 ) +* (-------+---) +* ( e'j 0 | 0 ) +* +* and to keep the main equality do the same for matrix Q. */ + +int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[], + const double val[]) +{ int m0 = lpf->m0; + int m = lpf->m; +#if _GLPLPF_DEBUG + double *B = lpf->B; +#endif + int n = lpf->n; + int *R_ptr = lpf->R_ptr; + int *R_len = lpf->R_len; + int *S_ptr = lpf->S_ptr; + int *S_len = lpf->S_len; + int *P_row = lpf->P_row; + int *P_col = lpf->P_col; + int *Q_row = lpf->Q_row; + int *Q_col = lpf->Q_col; + int v_ptr = lpf->v_ptr; + int *v_ind = lpf->v_ind; + double *v_val = lpf->v_val; + double *a = lpf->work2; /* new column */ + double *fg = lpf->work1, *f = fg, *g = fg + m0; + double *vw = lpf->work2, *v = vw, *w = vw + m0; + double *x = g, *y = w, z; + int i, ii, k, ret; + xassert(bh == bh); + if (!lpf->valid) + xfault("lpf_update_it: the factorization is not valid\n"); + if (!(1 <= j && j <= m)) + xfault("lpf_update_it: j = %d; column number out of range\n", + j); + xassert(0 <= m && m <= m0 + n); + /* check if the basis factorization can be expanded */ + if (n == lpf->n_max) + { lpf->valid = 0; + ret = LPF_ELIMIT; + goto done; + } + /* convert new j-th column of B to dense format */ + for (i = 1; i <= m; i++) + a[i] = 0.0; + for (k = 1; k <= len; k++) + { i = ind[k]; + if (!(1 <= i && i <= m)) + xfault("lpf_update_it: ind[%d] = %d; row number out of rang" + "e\n", k, i); + if (a[i] != 0.0) + xfault("lpf_update_it: ind[%d] = %d; duplicate row index no" + "t allowed\n", k, i); + if (val[k] == 0.0) + xfault("lpf_update_it: val[%d] = %g; zero element not allow" + "ed\n", k, val[k]); + a[i] = val[k]; + } +#if _GLPLPF_DEBUG + /* change column in the basis matrix for debugging */ + for (i = 1; i <= m; i++) + B[(i - 1) * m + j] = a[i]; +#endif + /* (f g) := inv(P) * (a 0) */ + for (i = 1; i <= m0+n; i++) + fg[i] = ((ii = P_col[i]) <= m ? a[ii] : 0.0); + /* (v w) := Q * (ej 0) */ + for (i = 1; i <= m0+n; i++) vw[i] = 0.0; + vw[Q_col[j]] = 1.0; + /* f1 := inv(L0) * f (new column of R) */ + luf_f_solve(lpf->luf, 0, f); + /* v1 := inv(U'0) * v (new row of S) */ + luf_v_solve(lpf->luf, 1, v); + /* we need at most 2 * m0 available locations in the SVA to store + new column of matrix R and new row of matrix S */ + if (lpf->v_size < v_ptr + m0 + m0) + { enlarge_sva(lpf, v_ptr + m0 + m0); + v_ind = lpf->v_ind; + v_val = lpf->v_val; + } + /* store new column of R */ + R_ptr[n+1] = v_ptr; + for (i = 1; i <= m0; i++) + { if (f[i] != 0.0) + v_ind[v_ptr] = i, v_val[v_ptr] = f[i], v_ptr++; + } + R_len[n+1] = v_ptr - lpf->v_ptr; + lpf->v_ptr = v_ptr; + /* store new row of S */ + S_ptr[n+1] = v_ptr; + for (i = 1; i <= m0; i++) + { if (v[i] != 0.0) + v_ind[v_ptr] = i, v_val[v_ptr] = v[i], v_ptr++; + } + S_len[n+1] = v_ptr - lpf->v_ptr; + lpf->v_ptr = v_ptr; + /* x := g - S * f1 (new column of C) */ + s_prod(lpf, x, -1.0, f); + /* y := w - R' * v1 (new row of C) */ + rt_prod(lpf, y, -1.0, v); + /* z := - v1 * f1 (new diagonal element of C) */ + z = 0.0; + for (i = 1; i <= m0; i++) z -= v[i] * f[i]; + /* update factorization of new matrix C */ + switch (scf_update_exp(lpf->scf, x, y, z)) + { case 0: + break; + case SCF_ESING: + lpf->valid = 0; + ret = LPF_ESING; + goto done; + case SCF_ELIMIT: + xassert(lpf != lpf); + default: + xassert(lpf != lpf); + } + /* expand matrix P */ + P_row[m0+n+1] = P_col[m0+n+1] = m0+n+1; + /* expand matrix Q */ + Q_row[m0+n+1] = Q_col[m0+n+1] = m0+n+1; + /* permute j-th and last (just added) column of matrix Q */ + i = Q_col[j], ii = Q_col[m0+n+1]; + Q_row[i] = m0+n+1, Q_col[m0+n+1] = i; + Q_row[ii] = j, Q_col[j] = ii; + /* increase the number of additional rows and columns */ + lpf->n++; + xassert(lpf->n <= lpf->n_max); + /* the factorization has been successfully updated */ + ret = 0; +done: /* return to the calling program */ + return ret; +} + +/*********************************************************************** +* NAME +* +* lpf_delete_it - delete LP basis factorization +* +* SYNOPSIS +* +* #include "glplpf.h" +* void lpf_delete_it(LPF *lpf) +* +* DESCRIPTION +* +* The routine lpf_delete_it deletes LP basis factorization specified +* by the parameter lpf and frees all memory allocated to this program +* object. */ + +void lpf_delete_it(LPF *lpf) +{ luf_delete_it(lpf->luf); +#if _GLPLPF_DEBUG + if (lpf->B != NULL) xfree(lpf->B); +#else + xassert(lpf->B == NULL); +#endif + if (lpf->R_ptr != NULL) xfree(lpf->R_ptr); + if (lpf->R_len != NULL) xfree(lpf->R_len); + if (lpf->S_ptr != NULL) xfree(lpf->S_ptr); + if (lpf->S_len != NULL) xfree(lpf->S_len); + if (lpf->scf != NULL) scf_delete_it(lpf->scf); + if (lpf->P_row != NULL) xfree(lpf->P_row); + if (lpf->P_col != NULL) xfree(lpf->P_col); + if (lpf->Q_row != NULL) xfree(lpf->Q_row); + if (lpf->Q_col != NULL) xfree(lpf->Q_col); + if (lpf->v_ind != NULL) xfree(lpf->v_ind); + if (lpf->v_val != NULL) xfree(lpf->v_val); + if (lpf->work1 != NULL) xfree(lpf->work1); + if (lpf->work2 != NULL) xfree(lpf->work2); + xfree(lpf); + return; +} + +/* eof */