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