You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
1097 lines
33 KiB
1097 lines
33 KiB
/* 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, 2011, 2013 Andrew Makhorin, Department for Applied
|
|
* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
|
|
* reserved. E-mail: <mao@gnu.org>.
|
|
*
|
|
* 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 <http://www.gnu.org/licenses/>.
|
|
***********************************************************************/
|
|
|
|
#include "env.h"
|
|
#include "glplpf.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;
|
|
#if 0 /* 06/VI-2013 */
|
|
lpf->luf = luf_create_it();
|
|
#else
|
|
lpf->lufint = lufint_create();
|
|
#endif
|
|
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;
|
|
#if 0 /* 11/VIII-2013 */
|
|
lpf->scf = NULL;
|
|
#else
|
|
lpf->ifu.n_max = 0;
|
|
lpf->ifu.n = 0;
|
|
lpf->ifu.f = NULL;
|
|
lpf->ifu.u = NULL;
|
|
lpf->t_opt = 0;
|
|
#endif
|
|
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 0 /* 11/VIII-2013 */
|
|
if (lpf->scf == NULL)
|
|
lpf->scf = scf_create_it(lpf->n_max);
|
|
#else
|
|
if (lpf->ifu.n_max == 0)
|
|
{ int n_max = lpf->n_max;
|
|
lpf->ifu.n_max = n_max;
|
|
lpf->ifu.n = 0;
|
|
xassert(n_max > 0);
|
|
xassert(lpf->ifu.f == NULL);
|
|
lpf->ifu.f = talloc(n_max * n_max, double);
|
|
xassert(lpf->ifu.u == NULL);
|
|
lpf->ifu.u = talloc(n_max * n_max, double);
|
|
lpf->t_opt = 0;
|
|
}
|
|
#endif
|
|
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 */
|
|
#if 0 /* 06/VI-2013 */
|
|
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);
|
|
}
|
|
#else
|
|
if (lufint_factorize(lpf->lufint, m, col, info) != 0)
|
|
{ ret = LPF_ESING;
|
|
goto done;
|
|
}
|
|
#endif
|
|
/* 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 */
|
|
#if 0 /* 11/VIII-2013 */
|
|
scf_reset_it(lpf->scf);
|
|
#else
|
|
lpf->ifu.n = 0;
|
|
#endif
|
|
/* 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 */
|
|
#if 0 /* 06/VI-2013 */
|
|
luf_f_solve(lpf->luf, 0, f);
|
|
#else
|
|
luf_f_solve(lpf->lufint->luf, f);
|
|
#endif
|
|
/* g1 := g - S * f1 */
|
|
s_prod(lpf, g, -1.0, f);
|
|
/* g2 := inv(C) * g1 */
|
|
#if 0 /* 11/VIII-2013 */
|
|
scf_solve_it(lpf->scf, 0, g);
|
|
#else
|
|
ifu_a_solve(&lpf->ifu, g, lpf->work2);
|
|
#endif
|
|
/* f2 := inv(U0) * (f1 - R * g2) */
|
|
r_prod(lpf, f, -1.0, g);
|
|
#if 0 /* 06/VI-2013 */
|
|
luf_v_solve(lpf->luf, 0, f);
|
|
#else
|
|
{ double *work = lpf->lufint->sgf->work;
|
|
luf_v_solve(lpf->lufint->luf, f, work);
|
|
memcpy(&f[1], &work[1], m0 * sizeof(double));
|
|
}
|
|
#endif
|
|
/* (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 */
|
|
#if 0 /* 06/VI-2013 */
|
|
luf_v_solve(lpf->luf, 1, f);
|
|
#else
|
|
{ double *work = lpf->lufint->sgf->work;
|
|
luf_vt_solve(lpf->lufint->luf, f, work);
|
|
memcpy(&f[1], &work[1], m0 * sizeof(double));
|
|
}
|
|
#endif
|
|
/* g1 := inv(C') * (g - R' * f1) */
|
|
rt_prod(lpf, g, -1.0, f);
|
|
#if 0 /* 11/VIII-2013 */
|
|
scf_solve_it(lpf->scf, 1, g);
|
|
#else
|
|
ifu_at_solve(&lpf->ifu, g, lpf->work2);
|
|
#endif
|
|
/* g2 := g1 */
|
|
g = g;
|
|
/* f2 := inv(L'0) * (f1 - S' * g2) */
|
|
st_prod(lpf, f, -1.0, g);
|
|
#if 0 /* 06/VI-2013 */
|
|
luf_f_solve(lpf->luf, 1, f);
|
|
#else
|
|
luf_ft_solve(lpf->lufint->luf, f);
|
|
#endif
|
|
/* (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) */
|
|
#if 0 /* 06/VI-2013 */
|
|
luf_f_solve(lpf->luf, 0, f);
|
|
#else
|
|
luf_f_solve(lpf->lufint->luf, f);
|
|
#endif
|
|
/* v1 := inv(U'0) * v (new row of S) */
|
|
#if 0 /* 06/VI-2013 */
|
|
luf_v_solve(lpf->luf, 1, v);
|
|
#else
|
|
{ double *work = lpf->lufint->sgf->work;
|
|
luf_vt_solve(lpf->lufint->luf, v, work);
|
|
memcpy(&v[1], &work[1], m0 * sizeof(double));
|
|
}
|
|
#endif
|
|
/* 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 */
|
|
#if 0 /* 11/VIII-2013 */
|
|
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);
|
|
}
|
|
#else
|
|
if (lpf->t_opt == SCF_TBG)
|
|
{ if (ifu_bg_update(&lpf->ifu, x, y, z) != 0)
|
|
#if 0
|
|
{ xprintf("Warning: insufficient accuracy (Bartels-Golub upda"
|
|
"te)\n");
|
|
#else
|
|
{
|
|
#endif
|
|
lpf->valid = 0;
|
|
ret = LPF_ESING;
|
|
goto done;
|
|
}
|
|
}
|
|
else if (lpf->t_opt == SCF_TGR)
|
|
{ if (ifu_gr_update(&lpf->ifu, x, y, z) != 0)
|
|
#if 0
|
|
{ xprintf("Warning: insufficient accuracy (Givens rotations u"
|
|
"pdate)\n");
|
|
#else
|
|
{
|
|
#endif
|
|
lpf->valid = 0;
|
|
ret = LPF_ESING;
|
|
goto done;
|
|
}
|
|
}
|
|
else
|
|
xassert(lpf != lpf);
|
|
#endif
|
|
/* 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)
|
|
#if 0 /* 06/VI-2013 */
|
|
{ luf_delete_it(lpf->luf);
|
|
#else
|
|
{ lufint_delete(lpf->lufint);
|
|
#endif
|
|
#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 0 /* 11/VIII-2013 */
|
|
if (lpf->scf != NULL) scf_delete_it(lpf->scf);
|
|
#else
|
|
if (lpf->ifu.f != NULL) tfree(lpf->ifu.f);
|
|
if (lpf->ifu.u != NULL) tfree(lpf->ifu.u);
|
|
#endif
|
|
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 */
|