diff options
Diffstat (limited to 'glpk-5.0/src/bflib')
| -rw-r--r-- | glpk-5.0/src/bflib/btf.c | 567 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/btf.h | 205 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/btfint.c | 405 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/btfint.h | 71 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/fhv.c | 584 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/fhv.h | 112 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/fhvint.c | 166 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/fhvint.h | 76 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/ifu.c | 390 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/ifu.h | 97 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/luf.c | 711 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/luf.h | 225 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/lufint.c | 180 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/lufint.h | 71 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/scf.c | 521 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/scf.h | 209 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/scfint.c | 253 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/scfint.h | 87 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/sgf.c | 1441 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/sgf.h | 201 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/sva.c | 570 | ||||
| -rw-r--r-- | glpk-5.0/src/bflib/sva.h | 159 |
22 files changed, 7301 insertions, 0 deletions
diff --git a/glpk-5.0/src/bflib/btf.c b/glpk-5.0/src/bflib/btf.c new file mode 100644 index 0000000..d0a7905 --- /dev/null +++ b/glpk-5.0/src/bflib/btf.c @@ -0,0 +1,567 @@ +/* btf.c (sparse block triangular LU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "btf.h" +#include "env.h" +#include "luf.h" +#include "mc13d.h" +#include "mc21a.h" + +/*********************************************************************** +* btf_store_a_cols - store pattern of matrix A in column-wise format +* +* This routine stores the pattern (that is, only indices of non-zero +* elements) of the original matrix A in column-wise format. +* +* On exit the routine returns the number of non-zeros in matrix A. */ + +int btf_store_a_cols(BTF *btf, int (*col)(void *info, int j, int ind[], + double val[]), void *info, int ind[], double val[]) +{ int n = btf->n; + SVA *sva = btf->sva; + int *sv_ind = sva->ind; + int ac_ref = btf->ac_ref; + int *ac_ptr = &sva->ptr[ac_ref-1]; + int *ac_len = &sva->len[ac_ref-1]; + int j, len, ptr, nnz; + nnz = 0; + for (j = 1; j <= n; j++) + { /* get j-th column */ + len = col(info, j, ind, val); + xassert(0 <= len && len <= n); + /* reserve locations for j-th column */ + if (len > 0) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + } + sva_reserve_cap(sva, ac_ref+(j-1), len); + } + /* store pattern of j-th column */ + ptr = ac_ptr[j]; + memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int)); + ac_len[j] = len; + nnz += len; + } + return nnz; +} + +/*********************************************************************** +* btf_make_blocks - permutations to block triangular form +* +* This routine analyzes the pattern of the original matrix A and +* determines permutation matrices P and Q such that A = P * A~* Q, +* where A~ is an upper block triangular matrix. +* +* On exit the routine returns symbolic rank of matrix A. */ + +int btf_make_blocks(BTF *btf) +{ int n = btf->n; + SVA *sva = btf->sva; + int *sv_ind = sva->ind; + int *pp_ind = btf->pp_ind; + int *pp_inv = btf->pp_inv; + int *qq_ind = btf->qq_ind; + int *qq_inv = btf->qq_inv; + int *beg = btf->beg; + int ac_ref = btf->ac_ref; + int *ac_ptr = &sva->ptr[ac_ref-1]; + int *ac_len = &sva->len[ac_ref-1]; + int i, j, rank, *iperm, *pr, *arp, *cv, *out, *ip, *lenr, *lowl, + *numb, *prev; + /* determine column permutation matrix M such that matrix A * M + * has zero-free diagonal */ + iperm = qq_inv; /* matrix M */ + pr = btf->p1_ind; /* working array */ + arp = btf->p1_inv; /* working array */ + cv = btf->q1_ind; /* working array */ + out = btf->q1_inv; /* working array */ + rank = mc21a(n, sv_ind, ac_ptr, ac_len, iperm, pr, arp, cv, out); + xassert(0 <= rank && rank <= n); + if (rank < n) + { /* A is structurally singular (rank is its symbolic rank) */ + goto done; + } + /* build pattern of matrix A * M */ + ip = pp_ind; /* working array */ + lenr = qq_ind; /* working array */ + for (j = 1; j <= n; j++) + { ip[j] = ac_ptr[iperm[j]]; + lenr[j] = ac_len[iperm[j]]; + } + /* determine symmetric permutation matrix S such that matrix + * S * (A * M) * S' = A~ is upper block triangular */ + lowl = btf->p1_ind; /* working array */ + numb = btf->p1_inv; /* working array */ + prev = btf->q1_ind; /* working array */ + btf->num = + mc13d(n, sv_ind, ip, lenr, pp_inv, beg, lowl, numb, prev); + xassert(beg[1] == 1); + beg[btf->num+1] = n+1; + /* A * M = S' * A~ * S ==> A = S' * A~ * (S * M') */ + /* determine permutation matrix P = S' */ + for (j = 1; j <= n; j++) + pp_ind[pp_inv[j]] = j; + /* determine permutation matrix Q = S * M' = P' * M' */ + for (i = 1; i <= n; i++) + qq_ind[i] = iperm[pp_inv[i]]; + for (i = 1; i <= n; i++) + qq_inv[qq_ind[i]] = i; +done: return rank; +} + +/*********************************************************************** +* btf_check_blocks - check structure of matrix A~ +* +* This routine checks that structure of upper block triangular matrix +* A~ is correct. +* +* NOTE: For testing/debugging only. */ + +void btf_check_blocks(BTF *btf) +{ int n = btf->n; + SVA *sva = btf->sva; + int *sv_ind = sva->ind; + int *pp_ind = btf->pp_ind; + int *pp_inv = btf->pp_inv; + int *qq_ind = btf->qq_ind; + int *qq_inv = btf->qq_inv; + int num = btf->num; + int *beg = btf->beg; + int ac_ref = btf->ac_ref; + int *ac_ptr = &sva->ptr[ac_ref-1]; + int *ac_len = &sva->len[ac_ref-1]; + int i, ii, j, jj, k, size, ptr, end, diag; + xassert(n > 0); + /* check permutation matrices P and Q */ + for (k = 1; k <= n; k++) + { xassert(1 <= pp_ind[k] && pp_ind[k] <= n); + xassert(pp_inv[pp_ind[k]] == k); + xassert(1 <= qq_ind[k] && qq_ind[k] <= n); + xassert(qq_inv[qq_ind[k]] == k); + } + /* check that matrix A~ is upper block triangular with non-zero + * diagonal */ + xassert(1 <= num && num <= n); + xassert(beg[1] == 1); + xassert(beg[num+1] == n+1); + /* walk thru blocks of A~ */ + for (k = 1; k <= num; k++) + { /* determine size of k-th block */ + size = beg[k+1] - beg[k]; + xassert(size >= 1); + /* walk thru columns of k-th block */ + for (jj = beg[k]; jj < beg[k+1]; jj++) + { diag = 0; + /* jj-th column of A~ = j-th column of A */ + j = qq_ind[jj]; + /* walk thru elements of j-th column of A */ + ptr = ac_ptr[j]; + end = ptr + ac_len[j]; + for (; ptr < end; ptr++) + { /* determine row index of a[i,j] */ + i = sv_ind[ptr]; + /* i-th row of A = ii-th row of A~ */ + ii = pp_ind[i]; + /* a~[ii,jj] should not be below k-th block */ + xassert(ii < beg[k+1]); + if (ii == jj) + { /* non-zero diagonal element of A~ encountered */ + diag = 1; + } + } + xassert(diag); + } + } + return; +} + +/*********************************************************************** +* btf_build_a_rows - build matrix A in row-wise format +* +* This routine builds the row-wise representation of matrix A in the +* right part of SVA using its column-wise representation. +* +* The working array len should have at least 1+n elements (len[0] is +* not used). */ + +void btf_build_a_rows(BTF *btf, int len[/*1+n*/]) +{ int n = btf->n; + SVA *sva = btf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int ar_ref = btf->ar_ref; + int *ar_ptr = &sva->ptr[ar_ref-1]; + int *ar_len = &sva->len[ar_ref-1]; + int ac_ref = btf->ac_ref; + int *ac_ptr = &sva->ptr[ac_ref-1]; + int *ac_len = &sva->len[ac_ref-1]; + int i, j, end, nnz, ptr, ptr1; + /* calculate the number of non-zeros in each row of matrix A and + * the total number of non-zeros */ + nnz = 0; + for (i = 1; i <= n; i++) + len[i] = 0; + for (j = 1; j <= n; j++) + { nnz += ac_len[j]; + for (end = (ptr = ac_ptr[j]) + ac_len[j]; ptr < end; ptr++) + len[sv_ind[ptr]]++; + } + /* we need at least nnz free locations in SVA */ + if (sva->r_ptr - sva->m_ptr < nnz) + { sva_more_space(sva, nnz); + sv_ind = sva->ind; + sv_val = sva->val; + } + /* reserve locations for rows of matrix A */ + for (i = 1; i <= n; i++) + { if (len[i] > 0) + sva_reserve_cap(sva, ar_ref-1+i, len[i]); + ar_len[i] = len[i]; + } + /* walk thru columns of matrix A and build its rows */ + for (j = 1; j <= n; j++) + { for (end = (ptr = ac_ptr[j]) + ac_len[j]; ptr < end; ptr++) + { i = sv_ind[ptr]; + sv_ind[ptr1 = ar_ptr[i] + (--len[i])] = j; + sv_val[ptr1] = sv_val[ptr]; + } + } + return; +} + +/*********************************************************************** +* btf_a_solve - solve system A * x = b +* +* This routine solves the system A * x = b, where A is the original +* matrix. +* +* On entry the array b should contain elements of the right-hand size +* vector b in locations b[1], ..., b[n], where n is the order of the +* matrix A. On exit the array x will contain elements of the solution +* vector in locations x[1], ..., x[n]. Note that the array b will be +* clobbered on exit. +* +* The routine also uses locations [1], ..., [max_size] of two working +* arrays w1 and w2, where max_size is the maximal size of diagonal +* blocks in BT-factorization (max_size <= n). */ + +void btf_a_solve(BTF *btf, double b[/*1+n*/], double x[/*1+n*/], + double w1[/*1+n*/], double w2[/*1+n*/]) +{ SVA *sva = btf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int *pp_inv = btf->pp_inv; + int *qq_ind = btf->qq_ind; + int num = btf->num; + int *beg = btf->beg; + int ac_ref = btf->ac_ref; + int *ac_ptr = &sva->ptr[ac_ref-1]; + int *ac_len = &sva->len[ac_ref-1]; + double *bb = w1; + double *xx = w2; + LUF luf; + int i, j, jj, k, beg_k, flag; + double t; + for (k = num; k >= 1; k--) + { /* determine order of diagonal block A~[k,k] */ + luf.n = beg[k+1] - (beg_k = beg[k]); + if (luf.n == 1) + { /* trivial case */ + /* solve system A~[k,k] * X[k] = B[k] */ + t = x[qq_ind[beg_k]] = + b[pp_inv[beg_k]] / btf->vr_piv[beg_k]; + /* substitute X[k] into other equations */ + if (t != 0.0) + { int ptr = ac_ptr[qq_ind[beg_k]]; + int end = ptr + ac_len[qq_ind[beg_k]]; + for (; ptr < end; ptr++) + b[sv_ind[ptr]] -= sv_val[ptr] * t; + } + } + else + { /* general case */ + /* construct B[k] */ + flag = 0; + for (i = 1; i <= luf.n; i++) + { if ((bb[i] = b[pp_inv[i + (beg_k-1)]]) != 0.0) + flag = 1; + } + /* solve system A~[k,k] * X[k] = B[k] */ + if (!flag) + { /* B[k] = 0, so X[k] = 0 */ + for (j = 1; j <= luf.n; j++) + x[qq_ind[j + (beg_k-1)]] = 0.0; + continue; + } + luf.sva = sva; + luf.fr_ref = btf->fr_ref + (beg_k-1); + luf.fc_ref = btf->fc_ref + (beg_k-1); + luf.vr_ref = btf->vr_ref + (beg_k-1); + luf.vr_piv = btf->vr_piv + (beg_k-1); + luf.vc_ref = btf->vc_ref + (beg_k-1); + luf.pp_ind = btf->p1_ind + (beg_k-1); + luf.pp_inv = btf->p1_inv + (beg_k-1); + luf.qq_ind = btf->q1_ind + (beg_k-1); + luf.qq_inv = btf->q1_inv + (beg_k-1); + luf_f_solve(&luf, bb); + luf_v_solve(&luf, bb, xx); + /* store X[k] and substitute it into other equations */ + for (j = 1; j <= luf.n; j++) + { jj = j + (beg_k-1); + t = x[qq_ind[jj]] = xx[j]; + if (t != 0.0) + { int ptr = ac_ptr[qq_ind[jj]]; + int end = ptr + ac_len[qq_ind[jj]]; + for (; ptr < end; ptr++) + b[sv_ind[ptr]] -= sv_val[ptr] * t; + } + } + } + } + return; +} + +/*********************************************************************** +* btf_at_solve - solve system A'* x = b +* +* This routine solves the system A'* x = b, where A' is a matrix +* transposed to the original matrix A. +* +* On entry the array b should contain elements of the right-hand size +* vector b in locations b[1], ..., b[n], where n is the order of the +* matrix A. On exit the array x will contain elements of the solution +* vector in locations x[1], ..., x[n]. Note that the array b will be +* clobbered on exit. +* +* The routine also uses locations [1], ..., [max_size] of two working +* arrays w1 and w2, where max_size is the maximal size of diagonal +* blocks in BT-factorization (max_size <= n). */ + +void btf_at_solve(BTF *btf, double b[/*1+n*/], double x[/*1+n*/], + double w1[/*1+n*/], double w2[/*1+n*/]) +{ SVA *sva = btf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int *pp_inv = btf->pp_inv; + int *qq_ind = btf->qq_ind; + int num = btf->num; + int *beg = btf->beg; + int ar_ref = btf->ar_ref; + int *ar_ptr = &sva->ptr[ar_ref-1]; + int *ar_len = &sva->len[ar_ref-1]; + double *bb = w1; + double *xx = w2; + LUF luf; + int i, j, jj, k, beg_k, flag; + double t; + for (k = 1; k <= num; k++) + { /* determine order of diagonal block A~[k,k] */ + luf.n = beg[k+1] - (beg_k = beg[k]); + if (luf.n == 1) + { /* trivial case */ + /* solve system A~'[k,k] * X[k] = B[k] */ + t = x[pp_inv[beg_k]] = + b[qq_ind[beg_k]] / btf->vr_piv[beg_k]; + /* substitute X[k] into other equations */ + if (t != 0.0) + { int ptr = ar_ptr[pp_inv[beg_k]]; + int end = ptr + ar_len[pp_inv[beg_k]]; + for (; ptr < end; ptr++) + b[sv_ind[ptr]] -= sv_val[ptr] * t; + } + } + else + { /* general case */ + /* construct B[k] */ + flag = 0; + for (i = 1; i <= luf.n; i++) + { if ((bb[i] = b[qq_ind[i + (beg_k-1)]]) != 0.0) + flag = 1; + } + /* solve system A~'[k,k] * X[k] = B[k] */ + if (!flag) + { /* B[k] = 0, so X[k] = 0 */ + for (j = 1; j <= luf.n; j++) + x[pp_inv[j + (beg_k-1)]] = 0.0; + continue; + } + luf.sva = sva; + luf.fr_ref = btf->fr_ref + (beg_k-1); + luf.fc_ref = btf->fc_ref + (beg_k-1); + luf.vr_ref = btf->vr_ref + (beg_k-1); + luf.vr_piv = btf->vr_piv + (beg_k-1); + luf.vc_ref = btf->vc_ref + (beg_k-1); + luf.pp_ind = btf->p1_ind + (beg_k-1); + luf.pp_inv = btf->p1_inv + (beg_k-1); + luf.qq_ind = btf->q1_ind + (beg_k-1); + luf.qq_inv = btf->q1_inv + (beg_k-1); + luf_vt_solve(&luf, bb, xx); + luf_ft_solve(&luf, xx); + /* store X[k] and substitute it into other equations */ + for (j = 1; j <= luf.n; j++) + { jj = j + (beg_k-1); + t = x[pp_inv[jj]] = xx[j]; + if (t != 0.0) + { int ptr = ar_ptr[pp_inv[jj]]; + int end = ptr + ar_len[pp_inv[jj]]; + for (; ptr < end; ptr++) + b[sv_ind[ptr]] -= sv_val[ptr] * t; + } + } + } + } + return; +} + +/*********************************************************************** +* btf_at_solve1 - solve system A'* y = e' to cause growth in y +* +* This routine is a special version of btf_at_solve. It solves the +* system A'* y = e' = e + delta e, where A' is a matrix transposed to +* the original matrix A, e is the specified right-hand side vector, +* and delta e is a vector of +1 and -1 chosen to cause growth in the +* solution vector y. +* +* On entry the array e should contain elements of the right-hand size +* vector e in locations e[1], ..., e[n], where n is the order of the +* matrix A. On exit the array y will contain elements of the solution +* vector in locations y[1], ..., y[n]. Note that the array e will be +* clobbered on exit. +* +* The routine also uses locations [1], ..., [max_size] of two working +* arrays w1 and w2, where max_size is the maximal size of diagonal +* blocks in BT-factorization (max_size <= n). */ + +void btf_at_solve1(BTF *btf, double e[/*1+n*/], double y[/*1+n*/], + double w1[/*1+n*/], double w2[/*1+n*/]) +{ SVA *sva = btf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int *pp_inv = btf->pp_inv; + int *qq_ind = btf->qq_ind; + int num = btf->num; + int *beg = btf->beg; + int ar_ref = btf->ar_ref; + int *ar_ptr = &sva->ptr[ar_ref-1]; + int *ar_len = &sva->len[ar_ref-1]; + double *ee = w1; + double *yy = w2; + LUF luf; + int i, j, jj, k, beg_k, ptr, end; + double e_k, y_k; + for (k = 1; k <= num; k++) + { /* determine order of diagonal block A~[k,k] */ + luf.n = beg[k+1] - (beg_k = beg[k]); + if (luf.n == 1) + { /* trivial case */ + /* determine E'[k] = E[k] + delta E[k] */ + e_k = e[qq_ind[beg_k]]; + e_k = (e_k >= 0.0 ? e_k + 1.0 : e_k - 1.0); + /* solve system A~'[k,k] * Y[k] = E[k] */ + y_k = y[pp_inv[beg_k]] = e_k / btf->vr_piv[beg_k]; + /* substitute Y[k] into other equations */ + ptr = ar_ptr[pp_inv[beg_k]]; + end = ptr + ar_len[pp_inv[beg_k]]; + for (; ptr < end; ptr++) + e[sv_ind[ptr]] -= sv_val[ptr] * y_k; + } + else + { /* general case */ + /* construct E[k] */ + for (i = 1; i <= luf.n; i++) + ee[i] = e[qq_ind[i + (beg_k-1)]]; + /* solve system A~'[k,k] * Y[k] = E[k] + delta E[k] */ + luf.sva = sva; + luf.fr_ref = btf->fr_ref + (beg_k-1); + luf.fc_ref = btf->fc_ref + (beg_k-1); + luf.vr_ref = btf->vr_ref + (beg_k-1); + luf.vr_piv = btf->vr_piv + (beg_k-1); + luf.vc_ref = btf->vc_ref + (beg_k-1); + luf.pp_ind = btf->p1_ind + (beg_k-1); + luf.pp_inv = btf->p1_inv + (beg_k-1); + luf.qq_ind = btf->q1_ind + (beg_k-1); + luf.qq_inv = btf->q1_inv + (beg_k-1); + luf_vt_solve1(&luf, ee, yy); + luf_ft_solve(&luf, yy); + /* store Y[k] and substitute it into other equations */ + for (j = 1; j <= luf.n; j++) + { jj = j + (beg_k-1); + y_k = y[pp_inv[jj]] = yy[j]; + ptr = ar_ptr[pp_inv[jj]]; + end = ptr + ar_len[pp_inv[jj]]; + for (; ptr < end; ptr++) + e[sv_ind[ptr]] -= sv_val[ptr] * y_k; + } + } + } + return; +} + +/*********************************************************************** +* btf_estimate_norm - estimate 1-norm of inv(A) +* +* This routine estimates 1-norm of inv(A) by one step of inverse +* iteration for the small singular vector as described in [1]. This +* involves solving two systems of equations: +* +* A'* y = e, +* +* A * z = y, +* +* where A' is a matrix transposed to A, and e is a vector of +1 and -1 +* chosen to cause growth in y. Then +* +* estimate 1-norm of inv(A) = (1-norm of z) / (1-norm of y) +* +* REFERENCES +* +* 1. G.E.Forsythe, M.A.Malcolm, C.B.Moler. Computer Methods for +* Mathematical Computations. Prentice-Hall, Englewood Cliffs, N.J., +* pp. 30-62 (subroutines DECOMP and SOLVE). */ + +double btf_estimate_norm(BTF *btf, double w1[/*1+n*/], double + w2[/*1+n*/], double w3[/*1+n*/], double w4[/*1+n*/]) +{ int n = btf->n; + double *e = w1; + double *y = w2; + double *z = w1; + int i; + double y_norm, z_norm; + /* compute y = inv(A') * e to cause growth in y */ + for (i = 1; i <= n; i++) + e[i] = 0.0; + btf_at_solve1(btf, e, y, w3, w4); + /* compute 1-norm of y = sum |y[i]| */ + y_norm = 0.0; + for (i = 1; i <= n; i++) + y_norm += (y[i] >= 0.0 ? +y[i] : -y[i]); + /* compute z = inv(A) * y */ + btf_a_solve(btf, y, z, w3, w4); + /* compute 1-norm of z = sum |z[i]| */ + z_norm = 0.0; + for (i = 1; i <= n; i++) + z_norm += (z[i] >= 0.0 ? +z[i] : -z[i]); + /* estimate 1-norm of inv(A) = (1-norm of z) / (1-norm of y) */ + return z_norm / y_norm; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/btf.h b/glpk-5.0/src/bflib/btf.h new file mode 100644 index 0000000..16e8c00 --- /dev/null +++ b/glpk-5.0/src/bflib/btf.h @@ -0,0 +1,205 @@ +/* btf.h (sparse block triangular LU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef BTF_H +#define BTF_H + +#include "sva.h" + +/*********************************************************************** +* The structure BTF describes BT-factorization, which is sparse block +* triangular LU-factorization. +* +* The BT-factorization has the following format: +* +* A = P * A~ * Q, (1) +* +* where A is a given (unsymmetric) square matrix, A~ is an upper block +* triangular matrix (see below), P and Q are permutation matrices. All +* the matrices have the same order n. +* +* The matrix A~, which is a permuted version of the original matrix A, +* has the following structure: +* +* A~[1,1] A~[1,2] ... A~[1,num-1] A~[1,num] +* +* A~[2,2] ... A~[2,num-1] A~[2,num] +* +* . . . . . . . . . (2) +* +* A~[num-1,num-1] A~[num-1,num] +* +* A~[num,num] +* +* where A~[i,j] is a submatrix called a "block," num is the number of +* blocks. Each diagonal block A~[k,k] is a non-singular square matrix, +* and each subdiagonal block A~[i,j], i > j, is a zero submatrix, thus +* A~ is an upper block triangular matrix. +* +* Permutation matrices P and Q are stored in ordinary arrays in both +* row- and column-like formats. +* +* The original matrix A is stored in both row- and column-wise sparse +* formats in the associated sparse vector area (SVA). Should note that +* elements of all diagonal blocks A~[k,k] in matrix A are set to zero +* (i.e. removed), so only elements of non-diagonal blocks are stored. +* +* Each diagonal block A~[k,k], 1 <= k <= num, is stored in the form of +* LU-factorization (see the module LUF). */ + +typedef struct BTF BTF; + +struct BTF +{ /* sparse block triangular LU-factorization */ + int n; + /* order of matrices A, A~, P, Q */ + SVA *sva; + /* associated sparse vector area used to store rows and columns + * of matrix A as well as sparse vectors for LU-factorizations of + * all diagonal blocks A~[k,k] */ + /*--------------------------------------------------------------*/ + /* matrix P */ + int *pp_ind; /* int pp_ind[1+n]; */ + /* pp_ind[i] = j means that P[i,j] = 1 */ + int *pp_inv; /* int pp_inv[1+n]; */ + /* pp_inv[j] = i means that P[i,j] = 1 */ + /* if i-th row of matrix A is i'-th row of matrix A~, then + * pp_ind[i] = i' and pp_inv[i'] = i */ + /*--------------------------------------------------------------*/ + /* matrix Q */ + int *qq_ind; /* int qq_ind[1+n]; */ + /* qq_ind[i] = j means that Q[i,j] = 1 */ + int *qq_inv; /* int qq_inv[1+n]; */ + /* qq_inv[j] = i means that Q[i,j] = 1 */ + /* if j-th column of matrix A is j'-th column of matrix A~, then + * qq_ind[j'] = j and qq_inv[j] = j' */ + /*--------------------------------------------------------------*/ + /* block triangular structure of matrix A~ */ + int num; + /* number of diagonal blocks, 1 <= num <= n */ + int *beg; /* int beg[1+num+1]; */ + /* beg[0] is not used; + * beg[k], 1 <= k <= num, is index of first row/column of k-th + * block of matrix A~; + * beg[num+1] is always n+1; + * note that order (size) of k-th diagonal block can be computed + * as beg[k+1] - beg[k] */ + /*--------------------------------------------------------------*/ + /* original matrix A in row-wise format */ + /* NOTE: elements of all diagonal blocks A~[k,k] are removed */ + int ar_ref; + /* reference number of sparse vector in SVA, which is the first + * row of matrix A */ +#if 0 + 0 + int *ar_ptr = &sva->ptr[ar_ref-1]; + /* ar_ptr[0] is not used; + * ar_ptr[i], 1 <= i <= n, is pointer to i-th row in SVA */ + int *ar_len = &sva->ptr[ar_ref-1]; + /* ar_len[0] is not used; + * ar_len[i], 1 <= i <= n, is length of i-th row */ +#endif + /*--------------------------------------------------------------*/ + /* original matrix A in column-wise format */ + /* NOTE: elements of all diagonal blocks A~[k,k] are removed */ + int ac_ref; + /* reference number of sparse vector in SVA, which is the first + * column of matrix A */ +#if 0 + 0 + int *ac_ptr = &sva->ptr[ac_ref-1]; + /* ac_ptr[0] is not used; + * ac_ptr[j], 1 <= j <= n, is pointer to j-th column in SVA */ + int *ac_len = &sva->ptr[ac_ref-1]; + /* ac_len[0] is not used; + * ac_len[j], 1 <= j <= n, is length of j-th column */ +#endif + /*--------------------------------------------------------------*/ + /* LU-factorizations of diagonal blocks A~[k,k] */ + /* to decrease overhead expenses similar arrays for all LUFs are + * packed into a single array; for example, elements fr_ptr[1], + * ..., fr_ptr[n1], where n1 = beg[2] - beg[1], are related to + * LUF for first diagonal block A~[1,1], elements fr_ptr[n1+1], + * ..., fr_ptr[n1+n2], where n2 = beg[3] - beg[2], are related to + * LUF for second diagonal block A~[2,2], etc.; in other words, + * elements related to LUF for k-th diagonal block A~[k,k] have + * indices beg[k], beg[k]+1, ..., beg[k+1]-1 */ + /* for details about LUF see description of the LUF module */ + int fr_ref; + /* reference number of sparse vector in SVA, which is the first + row of matrix F for first diagonal block A~[1,1] */ + int fc_ref; + /* reference number of sparse vector in SVA, which is the first + column of matrix F for first diagonal block A~[1,1] */ + int vr_ref; + /* reference number of sparse vector in SVA, which is the first + row of matrix V for first diagonal block A~[1,1] */ + double *vr_piv; /* double vr_piv[1+n]; */ + /* vr_piv[0] is not used; + vr_piv[1,...,n] are pivot elements for all diagonal blocks */ + int vc_ref; + /* reference number of sparse vector in SVA, which is the first + column of matrix V for first diagonal block A~[1,1] */ + int *p1_ind; /* int p1_ind[1+n]; */ + int *p1_inv; /* int p1_inv[1+n]; */ + int *q1_ind; /* int q1_ind[1+n]; */ + int *q1_inv; /* int q1_inv[1+n]; */ + /* permutation matrices P and Q for all diagonal blocks */ +}; + +#define btf_store_a_cols _glp_btf_store_a_cols +int btf_store_a_cols(BTF *btf, int (*col)(void *info, int j, int ind[], + double val[]), void *info, int ind[], double val[]); +/* store pattern of matrix A in column-wise format */ + +#define btf_make_blocks _glp_btf_make_blocks +int btf_make_blocks(BTF *btf); +/* permutations to block triangular form */ + +#define btf_check_blocks _glp_btf_check_blocks +void btf_check_blocks(BTF *btf); +/* check structure of matrix A~ */ + +#define btf_build_a_rows _glp_btf_build_a_rows +void btf_build_a_rows(BTF *btf, int len[/*1+n*/]); +/* build matrix A in row-wise format */ + +#define btf_a_solve _glp_btf_a_solve +void btf_a_solve(BTF *btf, double b[/*1+n*/], double x[/*1+n*/], + double w1[/*1+n*/], double w2[/*1+n*/]); +/* solve system A * x = b */ + +#define btf_at_solve _glp_btf_at_solve +void btf_at_solve(BTF *btf, double b[/*1+n*/], double x[/*1+n*/], + double w1[/*1+n*/], double w2[/*1+n*/]); +/* solve system A'* x = b */ + +#define btf_at_solve1 _glp_btf_at_solve1 +void btf_at_solve1(BTF *btf, double e[/*1+n*/], double y[/*1+n*/], + double w1[/*1+n*/], double w2[/*1+n*/]); +/* solve system A'* y = e' to cause growth in y */ + +#define btf_estimate_norm _glp_btf_estimate_norm +double btf_estimate_norm(BTF *btf, double w1[/*1+n*/], double + w2[/*1+n*/], double w3[/*1+n*/], double w4[/*1+n*/]); +/* estimate 1-norm of inv(A) */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/btfint.c b/glpk-5.0/src/bflib/btfint.c new file mode 100644 index 0000000..da23aa0 --- /dev/null +++ b/glpk-5.0/src/bflib/btfint.c @@ -0,0 +1,405 @@ +/* btfint.c (interface to BT-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "btfint.h" + +BTFINT *btfint_create(void) +{ /* create interface to BT-factorization */ + BTFINT *fi; + fi = talloc(1, BTFINT); + fi->n_max = 0; + fi->valid = 0; + fi->sva = NULL; + fi->btf = NULL; + fi->sgf = NULL; + fi->sva_n_max = fi->sva_size = 0; + fi->delta_n0 = fi->delta_n = 0; + fi->sgf_piv_tol = 0.10; + fi->sgf_piv_lim = 4; + fi->sgf_suhl = 1; + fi->sgf_eps_tol = DBL_EPSILON; + return fi; +} + +static void factorize_triv(BTFINT *fi, int k, int (*col)(void *info, + int j, int ind[], double val[]), void *info) +{ /* compute LU-factorization of diagonal block A~[k,k] and store + * corresponding columns of matrix A except elements of A~[k,k] + * (trivial case when the block has unity size) */ + SVA *sva = fi->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + BTF *btf = fi->btf; + int *pp_inv = btf->pp_inv; + int *qq_ind = btf->qq_ind; + int *beg = btf->beg; + int ac_ref = btf->ac_ref; + int *ac_ptr = &sva->ptr[ac_ref-1]; + int *ac_len = &sva->len[ac_ref-1]; + SGF *sgf = fi->sgf; + int *ind = (int *)sgf->vr_max; /* working array */ + double *val = sgf->work; /* working array */ + int i, j, t, len, ptr, beg_k; + /* diagonal block A~[k,k] has the only element in matrix A~, + * which is a~[beg[k],beg[k]] = a[i,j] */ + beg_k = beg[k]; + i = pp_inv[beg_k]; + j = qq_ind[beg_k]; + /* get j-th column of A */ + len = col(info, j, ind, val); + /* find element a[i,j] = a~[beg[k],beg[k]] in j-th column */ + for (t = 1; t <= len; t++) + { if (ind[t] == i) + break; + } + xassert(t <= len); + /* compute LU-factorization of diagonal block A~[k,k], where + * F = (1), V = (a[i,j]), P = Q = (1) (see the module LUF) */ +#if 1 /* FIXME */ + xassert(val[t] != 0.0); +#endif + btf->vr_piv[beg_k] = val[t]; + btf->p1_ind[beg_k] = btf->p1_inv[beg_k] = 1; + btf->q1_ind[beg_k] = btf->q1_inv[beg_k] = 1; + /* remove element a[i,j] = a~[beg[k],beg[k]] from j-th column */ + memmove(&ind[t], &ind[t+1], (len-t) * sizeof(int)); + memmove(&val[t], &val[t+1], (len-t) * sizeof(double)); + len--; + /* and store resulting j-th column of A into BTF */ + if (len > 0) + { /* reserve locations for j-th column of A */ + if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_reserve_cap(sva, ac_ref+(j-1), len); + /* store j-th column of A (except elements of A~[k,k]) */ + ptr = ac_ptr[j]; + memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int)); + memcpy(&sv_val[ptr], &val[1], len * sizeof(double)); + ac_len[j] = len; + } + return; +} + +static int factorize_block(BTFINT *fi, int k, int (*col)(void *info, + int j, int ind[], double val[]), void *info) +{ /* compute LU-factorization of diagonal block A~[k,k] and store + * corresponding columns of matrix A except elements of A~[k,k] + * (general case) */ + SVA *sva = fi->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + BTF *btf = fi->btf; + int *pp_ind = btf->pp_ind; + int *qq_ind = btf->qq_ind; + int *beg = btf->beg; + int ac_ref = btf->ac_ref; + int *ac_ptr = &sva->ptr[ac_ref-1]; + int *ac_len = &sva->len[ac_ref-1]; + SGF *sgf = fi->sgf; + int *ind = (int *)sgf->vr_max; /* working array */ + double *val = sgf->work; /* working array */ + LUF luf; + int *vc_ptr, *vc_len, *vc_cap; + int i, ii, j, jj, t, len, cnt, ptr, beg_k; + /* construct fake LUF for LU-factorization of A~[k,k] */ + sgf->luf = &luf; + luf.n = beg[k+1] - (beg_k = beg[k]); + luf.sva = sva; + luf.fr_ref = btf->fr_ref + (beg_k-1); + luf.fc_ref = btf->fc_ref + (beg_k-1); + luf.vr_ref = btf->vr_ref + (beg_k-1); + luf.vr_piv = btf->vr_piv + (beg_k-1); + luf.vc_ref = btf->vc_ref + (beg_k-1); + luf.pp_ind = btf->p1_ind + (beg_k-1); + luf.pp_inv = btf->p1_inv + (beg_k-1); + luf.qq_ind = btf->q1_ind + (beg_k-1); + luf.qq_inv = btf->q1_inv + (beg_k-1); + /* process columns of k-th block of matrix A~ */ + vc_ptr = &sva->ptr[luf.vc_ref-1]; + vc_len = &sva->len[luf.vc_ref-1]; + vc_cap = &sva->cap[luf.vc_ref-1]; + for (jj = 1; jj <= luf.n; jj++) + { /* jj-th column of A~ = j-th column of A */ + j = qq_ind[jj + (beg_k-1)]; + /* get j-th column of A */ + len = col(info, j, ind, val); + /* move elements of diagonal block A~[k,k] to the beginning of + * the column list */ + cnt = 0; + for (t = 1; t <= len; t++) + { /* i = row index of element a[i,j] */ + i = ind[t]; + /* i-th row of A = ii-th row of A~ */ + ii = pp_ind[i]; + if (ii >= beg_k) + { /* a~[ii,jj] = a[i,j] is in diagonal block A~[k,k] */ + double temp; + cnt++; + ind[t] = ind[cnt]; + ind[cnt] = ii - (beg_k-1); /* local index */ + temp = val[t], val[t] = val[cnt], val[cnt] = temp; + } + } + /* first cnt elements in the column list give jj-th column of + * diagonal block A~[k,k], which is initial matrix V in LUF */ + /* enlarge capacity of jj-th column of V = A~[k,k] */ + if (vc_cap[jj] < cnt) + { if (sva->r_ptr - sva->m_ptr < cnt) + { sva_more_space(sva, cnt); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, luf.vc_ref+(jj-1), cnt, 0); + } + /* store jj-th column of V = A~[k,k] */ + ptr = vc_ptr[jj]; + memcpy(&sv_ind[ptr], &ind[1], cnt * sizeof(int)); + memcpy(&sv_val[ptr], &val[1], cnt * sizeof(double)); + vc_len[jj] = cnt; + /* other (len-cnt) elements in the column list are stored in + * j-th column of the original matrix A */ + len -= cnt; + if (len > 0) + { /* reserve locations for j-th column of A */ + if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_reserve_cap(sva, ac_ref-1+j, len); + /* store j-th column of A (except elements of A~[k,k]) */ + ptr = ac_ptr[j]; + memcpy(&sv_ind[ptr], &ind[cnt+1], len * sizeof(int)); + memcpy(&sv_val[ptr], &val[cnt+1], len * sizeof(double)); + ac_len[j] = len; + } + } + /* compute LU-factorization of diagonal block A~[k,k]; may note + * that A~[k,k] is irreducible (strongly connected), so singleton + * phase will have no effect */ + k = sgf_factorize(sgf, 0 /* disable singleton phase */); + /* now left (dynamic) part of SVA should be empty (wichtig!) */ + xassert(sva->m_ptr == 1); + return k; +} + +int btfint_factorize(BTFINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info) +{ /* compute BT-factorization of specified matrix A */ + SVA *sva; + BTF *btf; + SGF *sgf; + int k, rank; + xassert(n > 0); + fi->valid = 0; + /* create sparse vector area (SVA), if necessary */ + sva = fi->sva; + if (sva == NULL) + { int sva_n_max = fi->sva_n_max; + int sva_size = fi->sva_size; + if (sva_n_max == 0) + sva_n_max = 6 * n; + if (sva_size == 0) + sva_size = 10 * n; + sva = fi->sva = sva_create_area(sva_n_max, sva_size); + } + /* allocate/reallocate underlying objects, if necessary */ + if (fi->n_max < n) + { int n_max = fi->n_max; + if (n_max == 0) + n_max = fi->n_max = n + fi->delta_n0; + else + n_max = fi->n_max = n + fi->delta_n; + xassert(n_max >= n); + /* allocate/reallocate block triangular factorization (BTF) */ + btf = fi->btf; + if (btf == NULL) + { btf = fi->btf = talloc(1, BTF); + memset(btf, 0, sizeof(BTF)); + btf->sva = sva; + } + else + { tfree(btf->pp_ind); + tfree(btf->pp_inv); + tfree(btf->qq_ind); + tfree(btf->qq_inv); + tfree(btf->beg); + tfree(btf->vr_piv); + tfree(btf->p1_ind); + tfree(btf->p1_inv); + tfree(btf->q1_ind); + tfree(btf->q1_inv); + } + btf->pp_ind = talloc(1+n_max, int); + btf->pp_inv = talloc(1+n_max, int); + btf->qq_ind = talloc(1+n_max, int); + btf->qq_inv = talloc(1+n_max, int); + btf->beg = talloc(1+n_max+1, int); + btf->vr_piv = talloc(1+n_max, double); + btf->p1_ind = talloc(1+n_max, int); + btf->p1_inv = talloc(1+n_max, int); + btf->q1_ind = talloc(1+n_max, int); + btf->q1_inv = talloc(1+n_max, int); + /* allocate/reallocate factorizer workspace (SGF) */ + /* (note that for SGF we could use the size of largest block + * rather than n_max) */ + sgf = fi->sgf; + sgf = fi->sgf; + if (sgf == NULL) + { sgf = fi->sgf = talloc(1, SGF); + memset(sgf, 0, sizeof(SGF)); + } + else + { tfree(sgf->rs_head); + tfree(sgf->rs_prev); + tfree(sgf->rs_next); + tfree(sgf->cs_head); + tfree(sgf->cs_prev); + tfree(sgf->cs_next); + tfree(sgf->vr_max); + tfree(sgf->flag); + tfree(sgf->work); + } + sgf->rs_head = talloc(1+n_max, int); + sgf->rs_prev = talloc(1+n_max, int); + sgf->rs_next = talloc(1+n_max, int); + sgf->cs_head = talloc(1+n_max, int); + sgf->cs_prev = talloc(1+n_max, int); + sgf->cs_next = talloc(1+n_max, int); + sgf->vr_max = talloc(1+n_max, double); + sgf->flag = talloc(1+n_max, char); + sgf->work = talloc(1+n_max, double); + } + btf = fi->btf; + btf->n = n; + sgf = fi->sgf; +#if 1 /* FIXME */ + /* initialize SVA */ + sva->n = 0; + sva->m_ptr = 1; + sva->r_ptr = sva->size + 1; + sva->head = sva->tail = 0; +#endif + /* store pattern of original matrix A in column-wise format */ + btf->ac_ref = sva_alloc_vecs(btf->sva, btf->n); + btf_store_a_cols(btf, col, info, btf->pp_ind, btf->vr_piv); +#ifdef GLP_DEBUG + sva_check_area(sva); +#endif + /* analyze pattern of original matrix A and determine permutation + * matrices P and Q such that A = P * A~* Q, where A~ is an upper + * block triangular matrix */ + rank = btf_make_blocks(btf); + if (rank != n) + { /* original matrix A is structurally singular */ + return 1; + } +#ifdef GLP_DEBUG + btf_check_blocks(btf); +#endif +#if 1 /* FIXME */ + /* initialize SVA */ + sva->n = 0; + sva->m_ptr = 1; + sva->r_ptr = sva->size + 1; + sva->head = sva->tail = 0; +#endif + /* allocate sparse vectors in SVA */ + btf->ar_ref = sva_alloc_vecs(btf->sva, btf->n); + btf->ac_ref = sva_alloc_vecs(btf->sva, btf->n); + btf->fr_ref = sva_alloc_vecs(btf->sva, btf->n); + btf->fc_ref = sva_alloc_vecs(btf->sva, btf->n); + btf->vr_ref = sva_alloc_vecs(btf->sva, btf->n); + btf->vc_ref = sva_alloc_vecs(btf->sva, btf->n); + /* setup factorizer control parameters */ + sgf->updat = 0; /* wichtig! */ + sgf->piv_tol = fi->sgf_piv_tol; + sgf->piv_lim = fi->sgf_piv_lim; + sgf->suhl = fi->sgf_suhl; + sgf->eps_tol = fi->sgf_eps_tol; + /* compute LU-factorizations of diagonal blocks A~[k,k] and also + * store corresponding columns of matrix A except elements of all + * blocks A~[k,k] */ + for (k = 1; k <= btf->num; k++) + { if (btf->beg[k+1] - btf->beg[k] == 1) + { /* trivial case (A~[k,k] has unity order) */ + factorize_triv(fi, k, col, info); + } + else + { /* general case */ + if (factorize_block(fi, k, col, info) != 0) + return 2; /* factorization of A~[k,k] failed */ + } + } +#ifdef GLP_DEBUG + sva_check_area(sva); +#endif + /* build row-wise representation of matrix A */ + btf_build_a_rows(fi->btf, fi->sgf->rs_head); +#ifdef GLP_DEBUG + sva_check_area(sva); +#endif + /* BT-factorization has been successfully computed */ + fi->valid = 1; + return 0; +} + +void btfint_delete(BTFINT *fi) +{ /* delete interface to BT-factorization */ + SVA *sva = fi->sva; + BTF *btf = fi->btf; + SGF *sgf = fi->sgf; + if (sva != NULL) + sva_delete_area(sva); + if (btf != NULL) + { tfree(btf->pp_ind); + tfree(btf->pp_inv); + tfree(btf->qq_ind); + tfree(btf->qq_inv); + tfree(btf->beg); + tfree(btf->vr_piv); + tfree(btf->p1_ind); + tfree(btf->p1_inv); + tfree(btf->q1_ind); + tfree(btf->q1_inv); + tfree(btf); + } + if (sgf != NULL) + { tfree(sgf->rs_head); + tfree(sgf->rs_prev); + tfree(sgf->rs_next); + tfree(sgf->cs_head); + tfree(sgf->cs_prev); + tfree(sgf->cs_next); + tfree(sgf->vr_max); + tfree(sgf->flag); + tfree(sgf->work); + tfree(sgf); + } + tfree(fi); + return; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/btfint.h b/glpk-5.0/src/bflib/btfint.h new file mode 100644 index 0000000..5c7f043 --- /dev/null +++ b/glpk-5.0/src/bflib/btfint.h @@ -0,0 +1,71 @@ +/* btfint.h (interface to BT-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef BTFINT_H +#define BTFINT_H + +#include "btf.h" +#include "sgf.h" + +typedef struct BTFINT BTFINT; + +struct BTFINT +{ /* interface to BT-factorization */ + int n_max; + /* maximal value of n (increased automatically) */ + int valid; + /* factorization is valid only if this flag is set */ + SVA *sva; + /* sparse vector area (SVA) */ + BTF *btf; + /* sparse block triangular LU-factorization */ + SGF *sgf; + /* sparse Gaussian factorizer workspace */ + /*--------------------------------------------------------------*/ + /* control parameters */ + int sva_n_max, sva_size; + /* parameters passed to sva_create_area */ + int delta_n0, delta_n; + /* if n_max = 0, set n_max = n + delta_n0 + * if n_max < n, set n_max = n + delta_n */ + double sgf_piv_tol; + int sgf_piv_lim; + int sgf_suhl; + double sgf_eps_tol; + /* factorizer control parameters */ +}; + +#define btfint_create _glp_btfint_create +BTFINT *btfint_create(void); +/* create interface to BT-factorization */ + +#define btfint_factorize _glp_btfint_factorize +int btfint_factorize(BTFINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info); +/* compute BT-factorization of specified matrix A */ + +#define btfint_delete _glp_btfint_delete +void btfint_delete(BTFINT *fi); +/* delete interface to BT-factorization */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/fhv.c b/glpk-5.0/src/bflib/fhv.c new file mode 100644 index 0000000..61d7be2 --- /dev/null +++ b/glpk-5.0/src/bflib/fhv.c @@ -0,0 +1,584 @@ +/* fhv.c (sparse updatable FHV-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "fhv.h" + +/*********************************************************************** +* fhv_ft_update - update FHV-factorization (Forrest-Tomlin) +* +* This routine updates FHV-factorization of the original matrix A +* after replacing its j-th column by a new one. The routine is based +* on the method proposed by Forrest and Tomlin [1]. +* +* The parameter q specifies the number of column of A, which has been +* replaced, 1 <= q <= n, where n is the order of A. +* +* Row indices and numerical values of non-zero elements of the new +* j-th column of A should be placed in locations aq_ind[1], ..., +* aq_ind[aq_len] and aq_val[1], ..., aq_val[aq_len], respectively, +* where aq_len is the number of non-zeros. Neither zero nor duplicate +* elements are allowed. +* +* The working arrays ind, val, and work should have at least 1+n +* elements (0-th elements are not used). +* +* RETURNS +* +* 0 The factorization has been successfully updated. +* +* 1 New matrix U = P'* V * Q' is upper triangular with zero diagonal +* element u[s,s]. (Elimination was not performed.) +* +* 2 New matrix U = P'* V * Q' is upper triangular, and its diagonal +* element u[s,s] or u[t,t] is too small in magnitude. (Elimination +* was not performed.) +* +* 3 The same as 2, but after performing elimination. +* +* 4 The factorization has not been updated, because maximal number of +* updates has been reached. +* +* 5 Accuracy test failed for the updated factorization. +* +* BACKGROUND +* +* The routine is based on the updating method proposed by Forrest and +* Tomlin [1]. +* +* Let q-th column of the original matrix A have been replaced by new +* column A[q]. Then, to keep the equality A = F * H * V, q-th column +* of matrix V should be replaced by column V[q] = inv(F * H) * A[q]. +* From the standpoint of matrix U = P'* V * Q' such replacement is +* equivalent to replacement of s-th column of matrix U, where s is +* determined from q by permutation matrix Q. Thus, matrix U loses its +* upper triangular form and becomes the following: +* +* 1 s t n +* 1 x x * x x x x x x +* . x * x x x x x x +* s . . * x x x x x x +* . . * x x x x x x +* . . * . x x x x x +* . . * . . x x x x +* t . . * . . . x x x +* . . . . . . . x x +* n . . . . . . . . x +* +* where t is largest row index of a non-zero element in s-th column. +* +* The routine makes matrix U upper triangular as follows. First, it +* moves rows and columns s+1, ..., t by one position to the left and +* upwards, resp., and moves s-th row and s-th column to position t. +* Due to such symmetric permutations matrix U becomes the following +* (note that all diagonal elements remain on the diagonal, and element +* u[s,s] becomes u[t,t]): +* +* 1 s t n +* 1 x x x x x x * x x +* . x x x x x * x x +* s . . x x x x * x x +* . . . x x x * x x +* . . . . x x * x x +* . . . . . x * x x +* t . . x x x x * x x +* . . . . . . . x x +* n . . . . . . . . x +* +* Then the routine performs gaussian elimination to eliminate +* subdiagonal elements u[t,s], ..., u[t,t-1] using diagonal elements +* u[s,s], ..., u[t-1,t-1] as pivots. During the elimination process +* the routine permutes neither rows nor columns, so only t-th row is +* changed. Should note that actually all operations are performed on +* matrix V = P * U * Q, since matrix U is not stored. +* +* To keep the equality A = F * H * V, the routine appends new row-like +* factor H[k] to matrix H, and every time it applies elementary +* gaussian transformation to eliminate u[t,j'] = v[p,j] using pivot +* u[j',j'] = v[i,j], it also adds new element f[p,j] = v[p,j] / v[i,j] +* (gaussian multiplier) to factor H[k], which initially is a unity +* matrix. At the end of elimination process the row-like factor H[k] +* may look as follows: +* +* 1 n 1 s t n +* 1 1 . . . . . . . . 1 1 . . . . . . . . +* . 1 . . . . . . . . 1 . . . . . . . +* . . 1 . . . . . . s . . 1 . . . . . . +* p . x x 1 . x . x . . . . 1 . . . . . +* . . . . 1 . . . . . . . . 1 . . . . +* . . . . . 1 . . . . . . . . 1 . . . +* . . . . . . 1 . . t . . x x x x 1 . . +* . . . . . . . 1 . . . . . . . . 1 . +* n . . . . . . . . 1 n . . . . . . . . 1 +* +* H[k] inv(P) * H[k] * P +* +* If, however, s = t, no elimination is needed, in which case no new +* row-like factor is created. +* +* REFERENCES +* +* 1. J.J.H.Forrest and J.A.Tomlin, "Updated triangular factors of the +* basis to maintain sparsity in the product form simplex method," +* Math. Prog. 2 (1972), pp. 263-78. */ + +int fhv_ft_update(FHV *fhv, int q, int aq_len, const int aq_ind[], + const double aq_val[], int ind[/*1+n*/], double val[/*1+n*/], + double work[/*1+n*/]) +{ LUF *luf = fhv->luf; + int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int *vr_cap = &sva->cap[vr_ref-1]; + double *vr_piv = luf->vr_piv; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *vc_cap = &sva->cap[vc_ref-1]; + int *pp_ind = luf->pp_ind; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int *qq_inv = luf->qq_inv; + int *hh_ind = fhv->hh_ind; + int hh_ref = fhv->hh_ref; + int *hh_ptr = &sva->ptr[hh_ref-1]; + int *hh_len = &sva->len[hh_ref-1]; +#if 1 /* FIXME */ + const double eps_tol = DBL_EPSILON; + const double vpq_tol = 1e-5; + const double err_tol = 1e-10; +#endif + int end, i, i_end, i_ptr, j, j_end, j_ptr, k, len, nnz, p, p_end, + p_ptr, ptr, q_end, q_ptr, s, t; + double f, vpq, temp; + /*--------------------------------------------------------------*/ + /* replace current q-th column of matrix V by new one */ + /*--------------------------------------------------------------*/ + xassert(1 <= q && q <= n); + /* convert new q-th column of matrix A to dense format */ + for (i = 1; i <= n; i++) + val[i] = 0.0; + xassert(0 <= aq_len && aq_len <= n); + for (k = 1; k <= aq_len; k++) + { i = aq_ind[k]; + xassert(1 <= i && i <= n); + xassert(val[i] == 0.0); + xassert(aq_val[k] != 0.0); + val[i] = aq_val[k]; + } + /* compute new q-th column of matrix V: + * new V[q] = inv(F * H) * (new A[q]) */ + luf->pp_ind = fhv->p0_ind; + luf->pp_inv = fhv->p0_inv; + luf_f_solve(luf, val); + luf->pp_ind = pp_ind; + luf->pp_inv = pp_inv; + fhv_h_solve(fhv, val); + /* q-th column of V = s-th column of U */ + s = qq_inv[q]; + /* determine row number of element v[p,q] that corresponds to + * diagonal element u[s,s] */ + p = pp_inv[s]; + /* convert new q-th column of V to sparse format; + * element v[p,q] = u[s,s] is not included in the element list + * and stored separately */ + vpq = 0.0; + len = 0; + for (i = 1; i <= n; i++) + { temp = val[i]; +#if 1 /* FIXME */ + if (-eps_tol < temp && temp < +eps_tol) +#endif + /* nop */; + else if (i == p) + vpq = temp; + else + { ind[++len] = i; + val[len] = temp; + } + } + /* clear q-th column of matrix V */ + for (q_end = (q_ptr = vc_ptr[q]) + vc_len[q]; + q_ptr < q_end; q_ptr++) + { /* get row index of v[i,q] */ + i = sv_ind[q_ptr]; + /* find and remove v[i,q] from i-th row */ + for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + sv_ind[i_ptr] != q; i_ptr++) + /* nop */; + xassert(i_ptr < i_end); + sv_ind[i_ptr] = sv_ind[i_end-1]; + sv_val[i_ptr] = sv_val[i_end-1]; + vr_len[i]--; + } + /* now q-th column of matrix V is empty */ + vc_len[q] = 0; + /* put new q-th column of V (except element v[p,q] = u[s,s]) in + * column-wise format */ + if (len > 0) + { if (vc_cap[q] < len) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, vc_ref-1+q, len, 0); + } + ptr = vc_ptr[q]; + memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int)); + memcpy(&sv_val[ptr], &val[1], len * sizeof(double)); + vc_len[q] = len; + } + /* put new q-th column of V (except element v[p,q] = u[s,s]) in + * row-wise format, and determine largest row number t such that + * u[s,t] != 0 */ + t = (vpq == 0.0 ? 0 : s); + for (k = 1; k <= len; k++) + { /* get row index of v[i,q] */ + i = ind[k]; + /* put v[i,q] to i-th row */ + if (vr_cap[i] == vr_len[i]) + { /* reserve extra locations in i-th row to reduce further + * relocations of that row */ +#if 1 /* FIXME */ + int need = vr_len[i] + 5; +#endif + if (sva->r_ptr - sva->m_ptr < need) + { sva_more_space(sva, need); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, vr_ref-1+i, need, 0); + } + sv_ind[ptr = vr_ptr[i] + (vr_len[i]++)] = q; + sv_val[ptr] = val[k]; + /* v[i,q] is non-zero; increase t */ + if (t < pp_ind[i]) + t = pp_ind[i]; + } + /*--------------------------------------------------------------*/ + /* check if matrix U is already upper triangular */ + /*--------------------------------------------------------------*/ + /* check if there is a spike in s-th column of matrix U, which + * is q-th column of matrix V */ + if (s >= t) + { /* no spike; matrix U is already upper triangular */ + /* store its diagonal element u[s,s] = v[p,q] */ + vr_piv[p] = vpq; + if (s > t) + { /* matrix U is structurally singular, because its diagonal + * element u[s,s] = v[p,q] is exact zero */ + xassert(vpq == 0.0); + return 1; + } +#if 1 /* FIXME */ + else if (-vpq_tol < vpq && vpq < +vpq_tol) +#endif + { /* matrix U is not well conditioned, because its diagonal + * element u[s,s] = v[p,q] is too small in magnitude */ + return 2; + } + else + { /* normal case */ + return 0; + } + } + /*--------------------------------------------------------------*/ + /* perform implicit symmetric permutations of rows and columns */ + /* of matrix U */ + /*--------------------------------------------------------------*/ + /* currently v[p,q] = u[s,s] */ + xassert(p == pp_inv[s] && q == qq_ind[s]); + for (k = s; k < t; k++) + { pp_ind[pp_inv[k] = pp_inv[k+1]] = k; + qq_inv[qq_ind[k] = qq_ind[k+1]] = k; + } + /* now v[p,q] = u[t,t] */ + pp_ind[pp_inv[t] = p] = qq_inv[qq_ind[t] = q] = t; + /*--------------------------------------------------------------*/ + /* check if matrix U is already upper triangular */ + /*--------------------------------------------------------------*/ + /* check if there is a spike in t-th row of matrix U, which is + * p-th row of matrix V */ + for (p_end = (p_ptr = vr_ptr[p]) + vr_len[p]; + p_ptr < p_end; p_ptr++) + { if (qq_inv[sv_ind[p_ptr]] < t) + break; /* spike detected */ + } + if (p_ptr == p_end) + { /* no spike; matrix U is already upper triangular */ + /* store its diagonal element u[t,t] = v[p,q] */ + vr_piv[p] = vpq; +#if 1 /* FIXME */ + if (-vpq_tol < vpq && vpq < +vpq_tol) +#endif + { /* matrix U is not well conditioned, because its diagonal + * element u[t,t] = v[p,q] is too small in magnitude */ + return 2; + } + else + { /* normal case */ + return 0; + } + } + /*--------------------------------------------------------------*/ + /* copy p-th row of matrix V, which is t-th row of matrix U, to */ + /* working array */ + /*--------------------------------------------------------------*/ + /* copy p-th row of matrix V, including element v[p,q] = u[t,t], + * to the working array in dense format and remove these elements + * from matrix V; since no pivoting is used, only this row will + * change during elimination */ + for (j = 1; j <= n; j++) + work[j] = 0.0; + work[q] = vpq; + for (p_end = (p_ptr = vr_ptr[p]) + vr_len[p]; + p_ptr < p_end; p_ptr++) + { /* get column index of v[p,j] and store this element to the + * working array */ + work[j = sv_ind[p_ptr]] = sv_val[p_ptr]; + /* find and remove v[p,j] from j-th column */ + for (j_end = (j_ptr = vc_ptr[j]) + vc_len[j]; + sv_ind[j_ptr] != p; j_ptr++) + /* nop */; + xassert(j_ptr < j_end); + sv_ind[j_ptr] = sv_ind[j_end-1]; + sv_val[j_ptr] = sv_val[j_end-1]; + vc_len[j]--; + } + /* now p-th row of matrix V is temporarily empty */ + vr_len[p] = 0; + /*--------------------------------------------------------------*/ + /* perform gaussian elimination */ + /*--------------------------------------------------------------*/ + /* transform p-th row of matrix V stored in working array, which + * is t-th row of matrix U, to eliminate subdiagonal elements + * u[t,s], ..., u[t,t-1]; corresponding gaussian multipliers will + * form non-trivial row of new row-like factor */ + nnz = 0; /* number of non-zero gaussian multipliers */ + for (k = s; k < t; k++) + { /* diagonal element u[k,k] = v[i,j] is used as pivot */ + i = pp_inv[k], j = qq_ind[k]; + /* take subdiagonal element u[t,k] = v[p,j] */ + temp = work[j]; +#if 1 /* FIXME */ + if (-eps_tol < temp && temp < +eps_tol) + continue; +#endif + /* compute and save gaussian multiplier: + * f := u[t,k] / u[k,k] = v[p,j] / v[i,j] */ + ind[++nnz] = i; + val[nnz] = f = work[j] / vr_piv[i]; + /* gaussian transformation to eliminate u[t,k] = v[p,j]: + * (p-th row of V) := (p-th row of V) - f * (i-th row of V) */ + for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + i_ptr < i_end; i_ptr++) + work[sv_ind[i_ptr]] -= f * sv_val[i_ptr]; + } + /* now matrix U is again upper triangular */ +#if 1 /* FIXME */ + if (-vpq_tol < work[q] && work[q] < +vpq_tol) +#endif + { /* however, its new diagonal element u[t,t] = v[p,q] is too + * small in magnitude */ + return 3; + } + /*--------------------------------------------------------------*/ + /* create new row-like factor H[k] and add to eta file H */ + /*--------------------------------------------------------------*/ + /* (nnz = 0 means that all subdiagonal elements were too small + * in magnitude) */ + if (nnz > 0) + { if (fhv->nfs == fhv->nfs_max) + { /* maximal number of row-like factors has been reached */ + return 4; + } + k = ++(fhv->nfs); + hh_ind[k] = p; + /* store non-trivial row of H[k] in right (dynamic) part of + * SVA (diagonal unity element is not stored) */ + if (sva->r_ptr - sva->m_ptr < nnz) + { sva_more_space(sva, nnz); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_reserve_cap(sva, fhv->hh_ref-1+k, nnz); + ptr = hh_ptr[k]; + memcpy(&sv_ind[ptr], &ind[1], nnz * sizeof(int)); + memcpy(&sv_val[ptr], &val[1], nnz * sizeof(double)); + hh_len[k] = nnz; + } + /*--------------------------------------------------------------*/ + /* copy transformed p-th row of matrix V, which is t-th row of */ + /* matrix U, from working array back to matrix V */ + /*--------------------------------------------------------------*/ + /* copy elements of transformed p-th row of matrix V, which are + * non-diagonal elements u[t,t+1], ..., u[t,n] of matrix U, from + * working array to corresponding columns of matrix V (note that + * diagonal element u[t,t] = v[p,q] not copied); also transform + * p-th row of matrix V to sparse format */ + len = 0; + for (k = t+1; k <= n; k++) + { /* j-th column of V = k-th column of U */ + j = qq_ind[k]; + /* take non-diagonal element v[p,j] = u[t,k] */ + temp = work[j]; +#if 1 /* FIXME */ + if (-eps_tol < temp && temp < +eps_tol) + continue; +#endif + /* add v[p,j] to j-th column of matrix V */ + if (vc_cap[j] == vc_len[j]) + { /* reserve extra locations in j-th column to reduce further + * relocations of that column */ +#if 1 /* FIXME */ + int need = vc_len[j] + 5; +#endif + if (sva->r_ptr - sva->m_ptr < need) + { sva_more_space(sva, need); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, vc_ref-1+j, need, 0); + } + sv_ind[ptr = vc_ptr[j] + (vc_len[j]++)] = p; + sv_val[ptr] = temp; + /* store element v[p,j] = u[t,k] to working sparse vector */ + ind[++len] = j; + val[len] = temp; + } + /* copy elements from working sparse vector to p-th row of matrix + * V (this row is currently empty) */ + if (vr_cap[p] < len) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, vr_ref-1+p, len, 0); + } + ptr = vr_ptr[p]; + memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int)); + memcpy(&sv_val[ptr], &val[1], len * sizeof(double)); + vr_len[p] = len; + /* store new diagonal element u[t,t] = v[p,q] */ + vr_piv[p] = work[q]; + /*--------------------------------------------------------------*/ + /* perform accuracy test (only if new H[k] was added) */ + /*--------------------------------------------------------------*/ + if (nnz > 0) + { /* copy p-th (non-trivial) row of row-like factor H[k] (except + * unity diagonal element) to working array in dense format */ + for (j = 1; j <= n; j++) + work[j] = 0.0; + k = fhv->nfs; + for (end = (ptr = hh_ptr[k]) + hh_len[k]; ptr < end; ptr++) + work[sv_ind[ptr]] = sv_val[ptr]; + /* compute inner product of p-th (non-trivial) row of matrix + * H[k] and q-th column of matrix V */ + temp = vr_piv[p]; /* 1 * v[p,q] */ + ptr = vc_ptr[q]; + end = ptr + vc_len[q]; + for (; ptr < end; ptr++) + temp += work[sv_ind[ptr]] * sv_val[ptr]; + /* inner product should be equal to element v[p,q] *before* + * matrix V was transformed */ + /* compute relative error */ + temp = fabs(vpq - temp) / (1.0 + fabs(vpq)); +#if 1 /* FIXME */ + if (temp > err_tol) +#endif + { /* relative error is too large */ + return 5; + } + } + /* factorization has been successfully updated */ + return 0; +} + +/*********************************************************************** +* fhv_h_solve - solve system H * x = b +* +* This routine solves the system H * x = b, where the matrix H is the +* middle factor of the sparse updatable FHV-factorization. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix H. On exit this array will contain elements of the solution +* vector x in the same locations. */ + +void fhv_h_solve(FHV *fhv, double x[/*1+n*/]) +{ SVA *sva = fhv->luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int nfs = fhv->nfs; + int *hh_ind = fhv->hh_ind; + int hh_ref = fhv->hh_ref; + int *hh_ptr = &sva->ptr[hh_ref-1]; + int *hh_len = &sva->len[hh_ref-1]; + int i, k, end, ptr; + double x_i; + for (k = 1; k <= nfs; k++) + { x_i = x[i = hh_ind[k]]; + for (end = (ptr = hh_ptr[k]) + hh_len[k]; ptr < end; ptr++) + x_i -= sv_val[ptr] * x[sv_ind[ptr]]; + x[i] = x_i; + } + return; +} + +/*********************************************************************** +* fhv_ht_solve - solve system H' * x = b +* +* This routine solves the system H' * x = b, where H' is a matrix +* transposed to the matrix H, which is the middle factor of the sparse +* updatable FHV-factorization. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix H. On exit this array will contain elements of the solution +* vector x in the same locations. */ + +void fhv_ht_solve(FHV *fhv, double x[/*1+n*/]) +{ SVA *sva = fhv->luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int nfs = fhv->nfs; + int *hh_ind = fhv->hh_ind; + int hh_ref = fhv->hh_ref; + int *hh_ptr = &sva->ptr[hh_ref-1]; + int *hh_len = &sva->len[hh_ref-1]; + int k, end, ptr; + double x_j; + for (k = nfs; k >= 1; k--) + { if ((x_j = x[hh_ind[k]]) == 0.0) + continue; + for (end = (ptr = hh_ptr[k]) + hh_len[k]; ptr < end; ptr++) + x[sv_ind[ptr]] -= sv_val[ptr] * x_j; + } + return; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/fhv.h b/glpk-5.0/src/bflib/fhv.h new file mode 100644 index 0000000..7b9d9a0 --- /dev/null +++ b/glpk-5.0/src/bflib/fhv.h @@ -0,0 +1,112 @@ +/* fhv.h (sparse updatable FHV-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef FHV_H +#define FHV_H + +#include "luf.h" + +/*********************************************************************** +* The structure FHV describes sparse updatable FHV-factorization. +* +* The FHV-factorization has the following format: +* +* A = F * H * V, (1) +* +* F = P0 * L * P0', (2) +* +* H = H[1] * H[2] * ... * H[nfs], (3) +* +* V = P * U * Q, (4) +* +* where: A is a given (unsymmetric) square matrix; F, H, V are matrix +* factors actually computed; L is a lower triangular matrix with unity +* diagonal; U is an upper tringular matrix; H[k], k = 1, 2, ..., nfs, +* is a row-like factor, which differs from unity matrix only in one +* row called a non-trivial row; P0, P, Q are permutation matrices; and +* P0' is a matrix transposed to P0. +* +* Matrices F, V, P, Q are stored in the underlying LUF object. +* +* Non-trivial rows of factors H[k] are stored as sparse vectors in the +* right (static) part of the sparse vector area (SVA). Note that unity +* diagonal elements of non-trivial rows are not stored. +* +* Matrix P0 is stored in the same way as matrix P. +* +* Matrices L and U are completely defined by matrices F, V, P, and Q, +* and therefore not stored explicitly. */ + +typedef struct FHV FHV; + +struct FHV +{ /* FHV-factorization */ + LUF *luf; + /* LU-factorization (contains matrices F, V, P, Q) */ + /*--------------------------------------------------------------*/ + /* matrix H in the form of eta file */ + int nfs_max; + /* maximal number of row-like factors (this limits the number of + * updates of the factorization) */ + int nfs; + /* current number of row-like factors, 0 <= nfs <= nfs_max */ + int *hh_ind; /* int hh_ind[1+nfs_max]; */ + /* hh_ind[0] is not used; + * hh_ind[k], 1 <= k <= nfs, is number of non-trivial row of + * factor H[k] */ + int hh_ref; + /* reference number of sparse vector in SVA, which is non-trivial + * row of factor H[1] */ +#if 0 + 0 + int *hh_ptr = &sva->ptr[hh_ref-1]; + /* hh_ptr[0] is not used; + * hh_ptr[k], 1 <= k <= nfs, is pointer to non-trivial row of + * factor H[k] */ + int *hh_len = &sva->len[hh_ref-1]; + /* hh_len[0] is not used; + * hh_len[k], 1 <= k <= nfs, is number of non-zero elements in + * non-trivial row of factor H[k] */ +#endif + /*--------------------------------------------------------------*/ + /* matrix P0 */ + int *p0_ind; /* int p0_ind[1+n]; */ + /* p0_ind[i] = j means that P0[i,j] = 1 */ + int *p0_inv; /* int p0_inv[1+n]; */ + /* p0_inv[j] = i means that P0[i,j] = 1 */ +}; + +#define fhv_ft_update _glp_fhv_ft_update +int fhv_ft_update(FHV *fhv, int q, int aq_len, const int aq_ind[], + const double aq_val[], int ind[/*1+n*/], double val[/*1+n*/], + double work[/*1+n*/]); +/* update FHV-factorization (Forrest-Tomlin) */ + +#define fhv_h_solve _glp_fhv_h_solve +void fhv_h_solve(FHV *fhv, double x[/*1+n*/]); +/* solve system H * x = b */ + +#define fhv_ht_solve _glp_fhv_ht_solve +void fhv_ht_solve(FHV *fhv, double x[/*1+n*/]); +/* solve system H' * x = b */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/fhvint.c b/glpk-5.0/src/bflib/fhvint.c new file mode 100644 index 0000000..b09071a --- /dev/null +++ b/glpk-5.0/src/bflib/fhvint.c @@ -0,0 +1,166 @@ +/* fhvint.c (interface to FHV-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "fhvint.h" + +FHVINT *fhvint_create(void) +{ /* create interface to FHV-factorization */ + FHVINT *fi; + fi = talloc(1, FHVINT); + memset(fi, 0, sizeof(FHVINT)); + fi->lufi = lufint_create(); + return fi; +} + +int fhvint_factorize(FHVINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info) +{ /* compute FHV-factorization of specified matrix A */ + int nfs_max, old_n_max, n_max, k, ret; + xassert(n > 0); + fi->valid = 0; + /* get required value of nfs_max */ + nfs_max = fi->nfs_max; + if (nfs_max == 0) + nfs_max = 100; + xassert(nfs_max > 0); + /* compute factorization of specified matrix A */ + old_n_max = fi->lufi->n_max; + fi->lufi->sva_n_max = 4 * n + nfs_max; + fi->lufi->sgf_updat = 1; + ret = lufint_factorize(fi->lufi, n, col, info); + n_max = fi->lufi->n_max; + /* allocate/reallocate arrays, if necessary */ + if (fi->fhv.nfs_max != nfs_max) + { if (fi->fhv.hh_ind != NULL) + tfree(fi->fhv.hh_ind); + fi->fhv.hh_ind = talloc(1+nfs_max, int); + } + if (old_n_max < n_max) + { if (fi->fhv.p0_ind != NULL) + tfree(fi->fhv.p0_ind); + if (fi->fhv.p0_inv != NULL) + tfree(fi->fhv.p0_inv); + fi->fhv.p0_ind = talloc(1+n_max, int); + fi->fhv.p0_inv = talloc(1+n_max, int); + } + /* initialize FHV-factorization */ + fi->fhv.luf = fi->lufi->luf; + fi->fhv.nfs_max = nfs_max; + /* H := I */ + fi->fhv.nfs = 0; + fi->fhv.hh_ref = sva_alloc_vecs(fi->lufi->sva, nfs_max); + /* P0 := P */ + for (k = 1; k <= n; k++) + { fi->fhv.p0_ind[k] = fi->fhv.luf->pp_ind[k]; + fi->fhv.p0_inv[k] = fi->fhv.luf->pp_inv[k]; + } + /* set validation flag */ + if (ret == 0) + fi->valid = 1; + return ret; +} + +int fhvint_update(FHVINT *fi, int j, int len, const int ind[], + const double val[]) +{ /* update FHV-factorization after replacing j-th column of A */ + SGF *sgf = fi->lufi->sgf; + int *ind1 = sgf->rs_next; + double *val1 = sgf->vr_max; + double *work = sgf->work; + int ret; + xassert(fi->valid); + ret = fhv_ft_update(&fi->fhv, j, len, ind, val, ind1, val1, work); + if (ret != 0) + fi->valid = 0; + return ret; +} + +void fhvint_ftran(FHVINT *fi, double x[]) +{ /* solve system A * x = b */ + FHV *fhv = &fi->fhv; + LUF *luf = fhv->luf; + int n = luf->n; + int *pp_ind = luf->pp_ind; + int *pp_inv = luf->pp_inv; + SGF *sgf = fi->lufi->sgf; + double *work = sgf->work; + xassert(fi->valid); + /* A = F * H * V */ + /* x = inv(A) * b = inv(V) * inv(H) * inv(F) * b */ + luf->pp_ind = fhv->p0_ind; + luf->pp_inv = fhv->p0_inv; + luf_f_solve(luf, x); + luf->pp_ind = pp_ind; + luf->pp_inv = pp_inv; + fhv_h_solve(fhv, x); + luf_v_solve(luf, x, work); + memcpy(&x[1], &work[1], n * sizeof(double)); + return; +} + +void fhvint_btran(FHVINT *fi, double x[]) +{ /* solve system A'* x = b */ + FHV *fhv = &fi->fhv; + LUF *luf = fhv->luf; + int n = luf->n; + int *pp_ind = luf->pp_ind; + int *pp_inv = luf->pp_inv; + SGF *sgf = fi->lufi->sgf; + double *work = sgf->work; + xassert(fi->valid); + /* A' = (F * H * V)' = V'* H'* F' */ + /* x = inv(A') * b = inv(F') * inv(H') * inv(V') * b */ + luf_vt_solve(luf, x, work); + fhv_ht_solve(fhv, work); + luf->pp_ind = fhv->p0_ind; + luf->pp_inv = fhv->p0_inv; + luf_ft_solve(luf, work); + luf->pp_ind = pp_ind; + luf->pp_inv = pp_inv; + memcpy(&x[1], &work[1], n * sizeof(double)); + return; +} + +double fhvint_estimate(FHVINT *fi) +{ /* estimate 1-norm of inv(A) */ + double norm; + xassert(fi->valid); + xassert(fi->fhv.nfs == 0); + norm = luf_estimate_norm(fi->fhv.luf, fi->lufi->sgf->vr_max, + fi->lufi->sgf->work); + return norm; +} + +void fhvint_delete(FHVINT *fi) +{ /* delete interface to FHV-factorization */ + lufint_delete(fi->lufi); + if (fi->fhv.hh_ind != NULL) + tfree(fi->fhv.hh_ind); + if (fi->fhv.p0_ind != NULL) + tfree(fi->fhv.p0_ind); + if (fi->fhv.p0_inv != NULL) + tfree(fi->fhv.p0_inv); + tfree(fi); + return; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/fhvint.h b/glpk-5.0/src/bflib/fhvint.h new file mode 100644 index 0000000..0be7945 --- /dev/null +++ b/glpk-5.0/src/bflib/fhvint.h @@ -0,0 +1,76 @@ +/* fhvint.h (interface to FHV-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef FHVINT_H +#define FHVINT_H + +#include "fhv.h" +#include "lufint.h" + +typedef struct FHVINT FHVINT; + +struct FHVINT +{ /* interface to FHV-factorization */ + int valid; + /* factorization is valid only if this flag is set */ + FHV fhv; + /* FHV-factorization */ + LUFINT *lufi; + /* interface to underlying LU-factorization */ + /*--------------------------------------------------------------*/ + /* control parameters */ + int nfs_max; + /* required maximal number of row-like factors */ +}; + +#define fhvint_create _glp_fhvint_create +FHVINT *fhvint_create(void); +/* create interface to FHV-factorization */ + +#define fhvint_factorize _glp_fhvint_factorize +int fhvint_factorize(FHVINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info); +/* compute FHV-factorization of specified matrix A */ + +#define fhvint_update _glp_fhvint_update +int fhvint_update(FHVINT *fi, int j, int len, const int ind[], + const double val[]); +/* update FHV-factorization after replacing j-th column of A */ + +#define fhvint_ftran _glp_fhvint_ftran +void fhvint_ftran(FHVINT *fi, double x[]); +/* solve system A * x = b */ + +#define fhvint_btran _glp_fhvint_btran +void fhvint_btran(FHVINT *fi, double x[]); +/* solve system A'* x = b */ + +#define fhvint_estimate _glp_fhvint_estimate +double fhvint_estimate(FHVINT *fi); +/* estimate 1-norm of inv(A) */ + +#define fhvint_delete _glp_fhvint_delete +void fhvint_delete(FHVINT *fi); +/* delete interface to FHV-factorization */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/ifu.c b/glpk-5.0/src/bflib/ifu.c new file mode 100644 index 0000000..5fcfa9a --- /dev/null +++ b/glpk-5.0/src/bflib/ifu.c @@ -0,0 +1,390 @@ +/* ifu.c (dense updatable IFU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "ifu.h" + +/*********************************************************************** +* ifu_expand - expand IFU-factorization +* +* This routine expands the IFU-factorization of the matrix A according +* to the following expansion of A: +* +* ( A c ) +* new A = ( ) +* ( r' d ) +* +* where c[1,...,n] is a new column, r[1,...,n] is a new row, and d is +* a new diagonal element. +* +* From the main equality F * A = U it follows that: +* +* ( F 0 ) ( A c ) ( FA Fc ) ( U Fc ) +* ( ) ( ) = ( ) = ( ), +* ( 0 1 ) ( r' d ) ( r' d ) ( r' d ) +* +* thus, +* +* ( F 0 ) ( U Fc ) +* new F = ( ), new U = ( ). +* ( 0 1 ) ( r' d ) +* +* Note that the resulting matrix U loses its upper triangular form due +* to row spike r', which should be eliminated. */ + +void ifu_expand(IFU *ifu, double c[/*1+n*/], double r[/*1+n*/], + double d) +{ /* non-optimized version */ + int n_max = ifu->n_max; + int n = ifu->n; + double *f_ = ifu->f; + double *u_ = ifu->u; + int i, j; + double t; +# define f(i,j) f_[(i)*n_max+(j)] +# define u(i,j) u_[(i)*n_max+(j)] + xassert(0 <= n && n < n_max); + /* adjust indexing */ + c++, r++; + /* set new zero column of matrix F */ + for (i = 0; i < n; i++) + f(i,n) = 0.0; + /* set new zero row of matrix F */ + for (j = 0; j < n; j++) + f(n,j) = 0.0; + /* set new unity diagonal element of matrix F */ + f(n,n) = 1.0; + /* set new column of matrix U to vector (old F) * c */ + for (i = 0; i < n; i++) + { /* u[i,n] := (i-th row of old F) * c */ + t = 0.0; + for (j = 0; j < n; j++) + t += f(i,j) * c[j]; + u(i,n) = t; + } + /* set new row of matrix U to vector r */ + for (j = 0; j < n; j++) + u(n,j) = r[j]; + /* set new diagonal element of matrix U to scalar d */ + u(n,n) = d; + /* increase factorization order */ + ifu->n++; +# undef f +# undef u + return; +} + +/*********************************************************************** +* ifu_bg_update - update IFU-factorization (Bartels-Golub) +* +* This routine updates IFU-factorization of the matrix A according to +* its expansion (see comments to the routine ifu_expand). The routine +* is based on the method proposed by Bartels and Golub [1]. +* +* RETURNS +* +* 0 The factorization has been successfully updated. +* +* 1 On some elimination step diagional element u[k,k] to be used as +* pivot is too small in magnitude. +* +* 2 Diagonal element u[n,n] is too small in magnitude (at the end of +* update). +* +* REFERENCES +* +* 1. R.H.Bartels, G.H.Golub, "The Simplex Method of Linear Programming +* Using LU-decomposition", Comm. ACM, 12, pp. 266-68, 1969. */ + +int ifu_bg_update(IFU *ifu, double c[/*1+n*/], double r[/*1+n*/], + double d) +{ /* non-optimized version */ + int n_max = ifu->n_max; + int n = ifu->n; + double *f_ = ifu->f; + double *u_ = ifu->u; +#if 1 /* FIXME */ + double tol = 1e-5; +#endif + int j, k; + double t; +# define f(i,j) f_[(i)*n_max+(j)] +# define u(i,j) u_[(i)*n_max+(j)] + /* expand factorization */ + ifu_expand(ifu, c, r, d); + /* NOTE: n keeps its old value */ + /* eliminate spike (non-zero subdiagonal elements) in last row of + * matrix U */ + for (k = 0; k < n; k++) + { /* if |u[k,k]| < |u[n,k]|, interchange k-th and n-th rows to + * provide |u[k,k]| >= |u[n,k]| for numeric stability */ + if (fabs(u(k,k)) < fabs(u(n,k))) + { /* interchange k-th and n-th rows of matrix U */ + for (j = k; j <= n; j++) + t = u(k,j), u(k,j) = u(n,j), u(n,j) = t; + /* interchange k-th and n-th rows of matrix F to keep the + * main equality F * A = U */ + for (j = 0; j <= n; j++) + t = f(k,j), f(k,j) = f(n,j), f(n,j) = t; + } + /* now |u[k,k]| >= |u[n,k]| */ + /* check if diagonal element u[k,k] can be used as pivot */ + if (fabs(u(k,k)) < tol) + { /* u[k,k] is too small in magnitude */ + return 1; + } + /* if u[n,k] = 0, elimination is not needed */ + if (u(n,k) == 0.0) + continue; + /* compute gaussian multiplier t = u[n,k] / u[k,k] */ + t = u(n,k) / u(k,k); + /* apply gaussian transformation to eliminate u[n,k] */ + /* (n-th row of U) := (n-th row of U) - t * (k-th row of U) */ + for (j = k+1; j <= n; j++) + u(n,j) -= t * u(k,j); + /* apply the same transformation to matrix F to keep the main + * equality F * A = U */ + for (j = 0; j <= n; j++) + f(n,j) -= t * f(k,j); + } + /* now matrix U is upper triangular */ + if (fabs(u(n,n)) < tol) + { /* u[n,n] is too small in magnitude */ + return 2; + } +# undef f +# undef u + return 0; +} + +/*********************************************************************** +* The routine givens computes the parameters of Givens plane rotation +* c = cos(teta) and s = sin(teta) such that: +* +* ( c -s ) ( a ) ( r ) +* ( ) ( ) = ( ) , +* ( s c ) ( b ) ( 0 ) +* +* where a and b are given scalars. +* +* REFERENCES +* +* G.H.Golub, C.F.Van Loan, "Matrix Computations", 2nd ed. */ + +static void givens(double a, double b, double *c, double *s) +{ /* non-optimized version */ + double t; + if (b == 0.0) + (*c) = 1.0, (*s) = 0.0; + else if (fabs(a) <= fabs(b)) + t = - a / b, (*s) = 1.0 / sqrt(1.0 + t * t), (*c) = (*s) * t; + else + t = - b / a, (*c) = 1.0 / sqrt(1.0 + t * t), (*s) = (*c) * t; + return; +} + +/*********************************************************************** +* ifu_gr_update - update IFU-factorization (Givens rotations) +* +* This routine updates IFU-factorization of the matrix A according to +* its expansion (see comments to the routine ifu_expand). The routine +* is based on Givens plane rotations [1]. +* +* RETURNS +* +* 0 The factorization has been successfully updated. +* +* 1 On some elimination step both elements u[k,k] and u[n,k] are too +* small in magnitude. +* +* 2 Diagonal element u[n,n] is too small in magnitude (at the end of +* update). +* +* REFERENCES +* +* 1. G.H.Golub, C.F.Van Loan, "Matrix Computations", 2nd ed. */ + +int ifu_gr_update(IFU *ifu, double c[/*1+n*/], double r[/*1+n*/], + double d) +{ /* non-optimized version */ + int n_max = ifu->n_max; + int n = ifu->n; + double *f_ = ifu->f; + double *u_ = ifu->u; +#if 1 /* FIXME */ + double tol = 1e-5; +#endif + int j, k; + double cs, sn; +# define f(i,j) f_[(i)*n_max+(j)] +# define u(i,j) u_[(i)*n_max+(j)] + /* expand factorization */ + ifu_expand(ifu, c, r, d); + /* NOTE: n keeps its old value */ + /* eliminate spike (non-zero subdiagonal elements) in last row of + * matrix U */ + for (k = 0; k < n; k++) + { /* check if elements u[k,k] and u[n,k] are eligible */ + if (fabs(u(k,k)) < tol && fabs(u(n,k)) < tol) + { /* both u[k,k] and u[n,k] are too small in magnitude */ + return 1; + } + /* if u[n,k] = 0, elimination is not needed */ + if (u(n,k) == 0.0) + continue; + /* compute parameters of Givens plane rotation */ + givens(u(k,k), u(n,k), &cs, &sn); + /* apply Givens rotation to k-th and n-th rows of matrix U to + * eliminate u[n,k] */ + for (j = k; j <= n; j++) + { double ukj = u(k,j), unj = u(n,j); + u(k,j) = cs * ukj - sn * unj; + u(n,j) = sn * ukj + cs * unj; + } + /* apply the same transformation to matrix F to keep the main + * equality F * A = U */ + for (j = 0; j <= n; j++) + { double fkj = f(k,j), fnj = f(n,j); + f(k,j) = cs * fkj - sn * fnj; + f(n,j) = sn * fkj + cs * fnj; + } + } + /* now matrix U is upper triangular */ + if (fabs(u(n,n)) < tol) + { /* u[n,n] is too small in magnitude */ + return 2; + } +# undef f +# undef u + return 0; +} + +/*********************************************************************** +* ifu_a_solve - solve system A * x = b +* +* This routine solves the system A * x = b, where the matrix A is +* specified by its IFU-factorization. +* +* Using the main equality F * A = U we have: +* +* A * x = b => F * A * x = F * b => U * x = F * b => +* +* x = inv(U) * F * b. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix A. On exit this array will contain elements of the solution +* vector x in the same locations. +* +* The working array w should have at least 1+n elements (0-th element +* is not used). */ + +void ifu_a_solve(IFU *ifu, double x[/*1+n*/], double w[/*1+n*/]) +{ /* non-optimized version */ + int n_max = ifu->n_max; + int n = ifu->n; + double *f_ = ifu->f; + double *u_ = ifu->u; + int i, j; + double t; +# define f(i,j) f_[(i)*n_max+(j)] +# define u(i,j) u_[(i)*n_max+(j)] + xassert(0 <= n && n <= n_max); + /* adjust indexing */ + x++, w++; + /* y := F * b */ + memcpy(w, x, n * sizeof(double)); + for (i = 0; i < n; i++) + { /* y[i] := (i-th row of F) * b */ + t = 0.0; + for (j = 0; j < n; j++) + t += f(i,j) * w[j]; + x[i] = t; + } + /* x := inv(U) * y */ + for (i = n-1; i >= 0; i--) + { t = x[i]; + for (j = i+1; j < n; j++) + t -= u(i,j) * x[j]; + x[i] = t / u(i,i); + } +# undef f +# undef u + return; +} + +/*********************************************************************** +* ifu_at_solve - solve system A'* x = b +* +* This routine solves the system A'* x = b, where A' is a matrix +* transposed to the matrix A, specified by its IFU-factorization. +* +* Using the main equality F * A = U, from which it follows that +* A'* F' = U', we have: +* +* A'* x = b => A'* F'* inv(F') * x = b => +* +* U'* inv(F') * x = b => inv(F') * x = inv(U') * b => +* +* x = F' * inv(U') * b. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix A. On exit this array will contain elements of the solution +* vector x in the same locations. +* +* The working array w should have at least 1+n elements (0-th element +* is not used). */ + +void ifu_at_solve(IFU *ifu, double x[/*1+n*/], double w[/*1+n*/]) +{ /* non-optimized version */ + int n_max = ifu->n_max; + int n = ifu->n; + double *f_ = ifu->f; + double *u_ = ifu->u; + int i, j; + double t; +# define f(i,j) f_[(i)*n_max+(j)] +# define u(i,j) u_[(i)*n_max+(j)] + xassert(0 <= n && n <= n_max); + /* adjust indexing */ + x++, w++; + /* y := inv(U') * b */ + for (i = 0; i < n; i++) + { t = (x[i] /= u(i,i)); + for (j = i+1; j < n; j++) + x[j] -= u(i,j) * t; + } + /* x := F'* y */ + for (j = 0; j < n; j++) + { /* x[j] := (j-th column of F) * y */ + t = 0.0; + for (i = 0; i < n; i++) + t += f(i,j) * x[i]; + w[j] = t; + } + memcpy(x, w, n * sizeof(double)); +# undef f +# undef u + return; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/ifu.h b/glpk-5.0/src/bflib/ifu.h new file mode 100644 index 0000000..8614696 --- /dev/null +++ b/glpk-5.0/src/bflib/ifu.h @@ -0,0 +1,97 @@ +/* ifu.h (dense updatable IFU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef IFU_H +#define IFU_H + +/*********************************************************************** +* The structure IFU describes dense updatable IFU-factorization. +* +* The IFU-factorization has the following format: +* +* A = inv(F) * U, (1) +* +* where A is a given (unsymmetric) nxn square matrix, F is a square +* matrix, U is an upper triangular matrix. Obviously, the equality (1) +* is equivalent to the following equality: +* +* F * A = U. (2) +* +* It is assumed that matrix A is small and dense, so matrices F and U +* are stored by rows in dense format as follows: +* +* 1 n n_max 1 n n_max +* 1 * * * * * * x x x x 1 * * * * * * x x x x +* * * * * * * x x x x ? * * * * * x x x x +* * * * * * * x x x x ? ? * * * * x x x x +* * * * * * * x x x x ? ? ? * * * x x x x +* * * * * * * x x x x ? ? ? ? * * x x x x +* n * * * * * * x x x x n ? ? ? ? ? * x x x x +* x x x x x x x x x x x x x x x x x x x x +* x x x x x x x x x x x x x x x x x x x x +* x x x x x x x x x x x x x x x x x x x x +* n_max x x x x x x x x x x n_max x x x x x x x x x x +* +* matrix F matrix U +* +* where '*' are matrix elements, '?' are unused locations, 'x' are +* reserved locations. */ + +typedef struct IFU IFU; + +struct IFU +{ /* IFU-factorization */ + int n_max; + /* maximal order of matrices A, F, U; n_max >= 1 */ + int n; + /* current order of matrices A, F, U; 0 <= n <= n_max */ + double *f; /* double f[n_max*n_max]; */ + /* matrix F stored by rows */ + double *u; /* double u[n_max*n_max]; */ + /* matrix U stored by rows */ +}; + +#define ifu_expand _glp_ifu_expand +void ifu_expand(IFU *ifu, double c[/*1+n*/], double r[/*1+n*/], + double d); +/* expand IFU-factorization */ + +#define ifu_bg_update _glp_ifu_bg_update +int ifu_bg_update(IFU *ifu, double c[/*1+n*/], double r[/*1+n*/], + double d); +/* update IFU-factorization (Bartels-Golub) */ + +#define ifu_gr_update _glp_ifu_gr_update +int ifu_gr_update(IFU *ifu, double c[/*1+n*/], double r[/*1+n*/], + double d); +/* update IFU-factorization (Givens rotations) */ + +#define ifu_a_solve _glp_ifu_a_solve +void ifu_a_solve(IFU *ifu, double x[/*1+n*/], double w[/*1+n*/]); +/* solve system A * x = b */ + +#define ifu_at_solve _glp_ifu_at_solve +void ifu_at_solve(IFU *ifu, double x[/*1+n*/], double w[/*1+n*/]); +/* solve system A'* x = b */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/luf.c b/glpk-5.0/src/bflib/luf.c new file mode 100644 index 0000000..3a1e628 --- /dev/null +++ b/glpk-5.0/src/bflib/luf.c @@ -0,0 +1,711 @@ +/* luf.c (sparse LU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "luf.h" + +/*********************************************************************** +* luf_store_v_cols - store matrix V = A in column-wise format +* +* This routine stores matrix V = A in column-wise format, where A is +* the original matrix to be factorized. +* +* On exit the routine returns the number of non-zeros in matrix V. */ + +int luf_store_v_cols(LUF *luf, int (*col)(void *info, int j, int ind[], + double val[]), void *info, int ind[], double val[]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *vc_cap = &sva->cap[vc_ref-1]; + int j, len, ptr, nnz; + nnz = 0; + for (j = 1; j <= n; j++) + { /* get j-th column */ + len = col(info, j, ind, val); + xassert(0 <= len && len <= n); + /* enlarge j-th column capacity */ + if (vc_cap[j] < len) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, vc_ref-1+j, len, 0); + } + /* store j-th column */ + ptr = vc_ptr[j]; + memcpy(&sv_ind[ptr], &ind[1], len * sizeof(int)); + memcpy(&sv_val[ptr], &val[1], len * sizeof(double)); + vc_len[j] = len; + nnz += len; + } + return nnz; +} + +/*********************************************************************** +* luf_check_all - check LU-factorization before k-th elimination step +* +* This routine checks that before performing k-th elimination step, +* 1 <= k <= n+1, all components of the LU-factorization are correct. +* +* In case of k = n+1, i.e. after last elimination step, it is assumed +* that rows of F and columns of V are *not* built yet. +* +* NOTE: For testing/debugging only. */ + +void luf_check_all(LUF *luf, int k) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fr_ref = luf->fr_ref; + int *fr_len = &sva->len[fr_ref-1]; + int fc_ref = luf->fc_ref; + int *fc_ptr = &sva->ptr[fc_ref-1]; + int *fc_len = &sva->len[fc_ref-1]; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *pp_ind = luf->pp_ind; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int *qq_inv = luf->qq_inv; + int i, ii, i_ptr, i_end, j, jj, j_ptr, j_end; + xassert(n > 0); + xassert(1 <= k && k <= n+1); + /* check permutation matrix P */ + for (i = 1; i <= n; i++) + { ii = pp_ind[i]; + xassert(1 <= ii && ii <= n); + xassert(pp_inv[ii] == i); + } + /* check permutation matrix Q */ + for (j = 1; j <= n; j++) + { jj = qq_inv[j]; + xassert(1 <= jj && jj <= n); + xassert(qq_ind[jj] == j); + } + /* check row-wise representation of matrix F */ + for (i = 1; i <= n; i++) + xassert(fr_len[i] == 0); + /* check column-wise representation of matrix F */ + for (j = 1; j <= n; j++) + { /* j-th column of F = jj-th column of L */ + jj = pp_ind[j]; + if (jj < k) + { j_ptr = fc_ptr[j]; + j_end = j_ptr + fc_len[j]; + for (; j_ptr < j_end; j_ptr++) + { i = sv_ind[j_ptr]; + xassert(1 <= i && i <= n); + ii = pp_ind[i]; /* f[i,j] = l[ii,jj] */ + xassert(ii > jj); + xassert(sv_val[j_ptr] != 0.0); + } + } + else /* jj >= k */ + xassert(fc_len[j] == 0); + } + /* check row-wise representation of matrix V */ + for (i = 1; i <= n; i++) + { /* i-th row of V = ii-th row of U */ + ii = pp_ind[i]; + i_ptr = vr_ptr[i]; + i_end = i_ptr + vr_len[i]; + for (; i_ptr < i_end; i_ptr++) + { j = sv_ind[i_ptr]; + xassert(1 <= j && j <= n); + jj = qq_inv[j]; /* v[i,j] = u[ii,jj] */ + if (ii < k) + xassert(jj > ii); + else /* ii >= k */ + { xassert(jj >= k); + /* find v[i,j] in j-th column */ + j_ptr = vc_ptr[j]; + j_end = j_ptr + vc_len[j]; + for (; sv_ind[j_ptr] != i; j_ptr++) + /* nop */; + xassert(j_ptr < j_end); + } + xassert(sv_val[i_ptr] != 0.0); + } + } + /* check column-wise representation of matrix V */ + for (j = 1; j <= n; j++) + { /* j-th column of V = jj-th column of U */ + jj = qq_inv[j]; + if (jj < k) + xassert(vc_len[j] == 0); + else /* jj >= k */ + { j_ptr = vc_ptr[j]; + j_end = j_ptr + vc_len[j]; + for (; j_ptr < j_end; j_ptr++) + { i = sv_ind[j_ptr]; + ii = pp_ind[i]; /* v[i,j] = u[ii,jj] */ + xassert(ii >= k); + /* find v[i,j] in i-th row */ + i_ptr = vr_ptr[i]; + i_end = i_ptr + vr_len[i]; + for (; sv_ind[i_ptr] != j; i_ptr++) + /* nop */; + xassert(i_ptr < i_end); + } + } + } + return; +} + +/*********************************************************************** +* luf_build_v_rows - build matrix V in row-wise format +* +* This routine builds the row-wise representation of matrix V in the +* left part of SVA using its column-wise representation. +* +* NOTE: On entry to the routine all rows of matrix V should have zero +* capacity. +* +* The working array len should have at least 1+n elements (len[0] is +* not used). */ + +void luf_build_v_rows(LUF *luf, int len[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int i, j, end, nnz, ptr, ptr1; + /* calculate the number of non-zeros in each row of matrix V and + * the total number of non-zeros */ + nnz = 0; + for (i = 1; i <= n; i++) + len[i] = 0; + for (j = 1; j <= n; j++) + { nnz += vc_len[j]; + for (end = (ptr = vc_ptr[j]) + vc_len[j]; ptr < end; ptr++) + len[sv_ind[ptr]]++; + } + /* we need at least nnz free locations in SVA */ + if (sva->r_ptr - sva->m_ptr < nnz) + { sva_more_space(sva, nnz); + sv_ind = sva->ind; + sv_val = sva->val; + } + /* reserve locations for rows of matrix V */ + for (i = 1; i <= n; i++) + { if (len[i] > 0) + sva_enlarge_cap(sva, vr_ref-1+i, len[i], 0); + vr_len[i] = len[i]; + } + /* walk thru column of matrix V and build its rows */ + for (j = 1; j <= n; j++) + { for (end = (ptr = vc_ptr[j]) + vc_len[j]; ptr < end; ptr++) + { i = sv_ind[ptr]; + sv_ind[ptr1 = vr_ptr[i] + (--len[i])] = j; + sv_val[ptr1] = sv_val[ptr]; + } + } + return; +} + +/*********************************************************************** +* luf_build_f_rows - build matrix F in row-wise format +* +* This routine builds the row-wise representation of matrix F in the +* right part of SVA using its column-wise representation. +* +* NOTE: On entry to the routine all rows of matrix F should have zero +* capacity. +* +* The working array len should have at least 1+n elements (len[0] is +* not used). */ + +void luf_build_f_rows(LUF *luf, int len[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fr_ref = luf->fr_ref; + int *fr_ptr = &sva->ptr[fr_ref-1]; + int *fr_len = &sva->len[fr_ref-1]; + int fc_ref = luf->fc_ref; + int *fc_ptr = &sva->ptr[fc_ref-1]; + int *fc_len = &sva->len[fc_ref-1]; + int i, j, end, nnz, ptr, ptr1; + /* calculate the number of non-zeros in each row of matrix F and + * the total number of non-zeros (except diagonal elements) */ + nnz = 0; + for (i = 1; i <= n; i++) + len[i] = 0; + for (j = 1; j <= n; j++) + { nnz += fc_len[j]; + for (end = (ptr = fc_ptr[j]) + fc_len[j]; ptr < end; ptr++) + len[sv_ind[ptr]]++; + } + /* we need at least nnz free locations in SVA */ + if (sva->r_ptr - sva->m_ptr < nnz) + { sva_more_space(sva, nnz); + sv_ind = sva->ind; + sv_val = sva->val; + } + /* reserve locations for rows of matrix F */ + for (i = 1; i <= n; i++) + { if (len[i] > 0) + sva_reserve_cap(sva, fr_ref-1+i, len[i]); + fr_len[i] = len[i]; + } + /* walk through columns of matrix F and build its rows */ + for (j = 1; j <= n; j++) + { for (end = (ptr = fc_ptr[j]) + fc_len[j]; ptr < end; ptr++) + { i = sv_ind[ptr]; + sv_ind[ptr1 = fr_ptr[i] + (--len[i])] = j; + sv_val[ptr1] = sv_val[ptr]; + } + } + return; +} + +/*********************************************************************** +* luf_build_v_cols - build matrix V in column-wise format +* +* This routine builds the column-wise representation of matrix V in +* the left (if the flag updat is set) or right (if the flag updat is +* clear) part of SVA using its row-wise representation. +* +* NOTE: On entry to the routine all columns of matrix V should have +* zero capacity. +* +* The working array len should have at least 1+n elements (len[0] is +* not used). */ + +void luf_build_v_cols(LUF *luf, int updat, int len[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int i, j, end, nnz, ptr, ptr1; + /* calculate the number of non-zeros in each column of matrix V + * and the total number of non-zeros (except pivot elements) */ + nnz = 0; + for (j = 1; j <= n; j++) + len[j] = 0; + for (i = 1; i <= n; i++) + { nnz += vr_len[i]; + for (end = (ptr = vr_ptr[i]) + vr_len[i]; ptr < end; ptr++) + len[sv_ind[ptr]]++; + } + /* we need at least nnz free locations in SVA */ + if (sva->r_ptr - sva->m_ptr < nnz) + { sva_more_space(sva, nnz); + sv_ind = sva->ind; + sv_val = sva->val; + } + /* reserve locations for columns of matrix V */ + for (j = 1; j <= n; j++) + { if (len[j] > 0) + { if (updat) + sva_enlarge_cap(sva, vc_ref-1+j, len[j], 0); + else + sva_reserve_cap(sva, vc_ref-1+j, len[j]); + } + vc_len[j] = len[j]; + } + /* walk through rows of matrix V and build its columns */ + for (i = 1; i <= n; i++) + { for (end = (ptr = vr_ptr[i]) + vr_len[i]; ptr < end; ptr++) + { j = sv_ind[ptr]; + sv_ind[ptr1 = vc_ptr[j] + (--len[j])] = i; + sv_val[ptr1] = sv_val[ptr]; + } + } + return; +} + +/*********************************************************************** +* luf_check_f_rc - check rows and columns of matrix F +* +* This routine checks that the row- and column-wise representations +* of matrix F are identical. +* +* NOTE: For testing/debugging only. */ + +void luf_check_f_rc(LUF *luf) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fr_ref = luf->fr_ref; + int *fr_ptr = &sva->ptr[fr_ref-1]; + int *fr_len = &sva->len[fr_ref-1]; + int fc_ref = luf->fc_ref; + int *fc_ptr = &sva->ptr[fc_ref-1]; + int *fc_len = &sva->len[fc_ref-1]; + int i, i_end, i_ptr, j, j_end, j_ptr; + /* walk thru rows of matrix F */ + for (i = 1; i <= n; i++) + { for (i_end = (i_ptr = fr_ptr[i]) + fr_len[i]; + i_ptr < i_end; i_ptr++) + { j = sv_ind[i_ptr]; + /* find element f[i,j] in j-th column of matrix F */ + for (j_end = (j_ptr = fc_ptr[j]) + fc_len[j]; + sv_ind[j_ptr] != i; j_ptr++) + /* nop */; + xassert(j_ptr < j_end); + xassert(sv_val[i_ptr] == sv_val[j_ptr]); + /* mark element f[i,j] */ + sv_ind[j_ptr] = -i; + } + } + /* walk thru column of matix F and check that all elements has + been marked */ + for (j = 1; j <= n; j++) + { for (j_end = (j_ptr = fc_ptr[j]) + fc_len[j]; + j_ptr < j_end; j_ptr++) + { xassert((i = sv_ind[j_ptr]) < 0); + /* unmark element f[i,j] */ + sv_ind[j_ptr] = -i; + } + } + return; +} + +/*********************************************************************** +* luf_check_v_rc - check rows and columns of matrix V +* +* This routine checks that the row- and column-wise representations +* of matrix V are identical. +* +* NOTE: For testing/debugging only. */ + +void luf_check_v_rc(LUF *luf) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int i, i_end, i_ptr, j, j_end, j_ptr; + /* walk thru rows of matrix V */ + for (i = 1; i <= n; i++) + { for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + i_ptr < i_end; i_ptr++) + { j = sv_ind[i_ptr]; + /* find element v[i,j] in j-th column of matrix V */ + for (j_end = (j_ptr = vc_ptr[j]) + vc_len[j]; + sv_ind[j_ptr] != i; j_ptr++) + /* nop */; + xassert(j_ptr < j_end); + xassert(sv_val[i_ptr] == sv_val[j_ptr]); + /* mark element v[i,j] */ + sv_ind[j_ptr] = -i; + } + } + /* walk thru column of matix V and check that all elements has + been marked */ + for (j = 1; j <= n; j++) + { for (j_end = (j_ptr = vc_ptr[j]) + vc_len[j]; + j_ptr < j_end; j_ptr++) + { xassert((i = sv_ind[j_ptr]) < 0); + /* unmark element v[i,j] */ + sv_ind[j_ptr] = -i; + } + } + return; +} + +/*********************************************************************** +* luf_f_solve - solve system F * x = b +* +* This routine solves the system F * x = b, where the matrix F is the +* left factor of the sparse LU-factorization. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix F. On exit this array will contain elements of the solution +* vector x in the same locations. */ + +void luf_f_solve(LUF *luf, double x[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fc_ref = luf->fc_ref; + int *fc_ptr = &sva->ptr[fc_ref-1]; + int *fc_len = &sva->len[fc_ref-1]; + int *pp_inv = luf->pp_inv; + int j, k, ptr, end; + double x_j; + for (k = 1; k <= n; k++) + { /* k-th column of L = j-th column of F */ + j = pp_inv[k]; + /* x[j] is already computed */ + /* walk thru j-th column of matrix F and substitute x[j] into + * other equations */ + if ((x_j = x[j]) != 0.0) + { for (end = (ptr = fc_ptr[j]) + fc_len[j]; ptr < end; ptr++) + x[sv_ind[ptr]] -= sv_val[ptr] * x_j; + } + } + return; +} + +/*********************************************************************** +* luf_ft_solve - solve system F' * x = b +* +* This routine solves the system F' * x = b, where F' is a matrix +* transposed to the matrix F, which is the left factor of the sparse +* LU-factorization. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix F. On exit this array will contain elements of the solution +* vector x in the same locations. */ + +void luf_ft_solve(LUF *luf, double x[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fr_ref = luf->fr_ref; + int *fr_ptr = &sva->ptr[fr_ref-1]; + int *fr_len = &sva->len[fr_ref-1]; + int *pp_inv = luf->pp_inv; + int i, k, ptr, end; + double x_i; + for (k = n; k >= 1; k--) + { /* k-th column of L' = i-th row of F */ + i = pp_inv[k]; + /* x[i] is already computed */ + /* walk thru i-th row of matrix F and substitute x[i] into + * other equations */ + if ((x_i = x[i]) != 0.0) + { for (end = (ptr = fr_ptr[i]) + fr_len[i]; ptr < end; ptr++) + x[sv_ind[ptr]] -= sv_val[ptr] * x_i; + } + } + return; +} + +/*********************************************************************** +* luf_v_solve - solve system V * x = b +* +* This routine solves the system V * x = b, where the matrix V is the +* right factor of the sparse LU-factorization. +* +* On entry the array b should contain elements of the right-hand side +* vector b in locations b[1], ..., b[n], where n is the order of the +* matrix V. On exit the array x will contain elements of the solution +* vector x in locations x[1], ..., x[n]. Note that the array b will be +* clobbered on exit. */ + +void luf_v_solve(LUF *luf, double b[/*1+n*/], double x[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + double *vr_piv = luf->vr_piv; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int i, j, k, ptr, end; + double x_j; + for (k = n; k >= 1; k--) + { /* k-th row of U = i-th row of V */ + /* k-th column of U = j-th column of V */ + i = pp_inv[k]; + j = qq_ind[k]; + /* compute x[j] = b[i] / u[k,k], where u[k,k] = v[i,j]; + * walk through j-th column of matrix V and substitute x[j] + * into other equations */ + if ((x_j = x[j] = b[i] / vr_piv[i]) != 0.0) + { for (end = (ptr = vc_ptr[j]) + vc_len[j]; ptr < end; ptr++) + b[sv_ind[ptr]] -= sv_val[ptr] * x_j; + } + } + return; +} + +/*********************************************************************** +* luf_vt_solve - solve system V' * x = b +* +* This routine solves the system V' * x = b, where V' is a matrix +* transposed to the matrix V, which is the right factor of the sparse +* LU-factorization. +* +* On entry the array b should contain elements of the right-hand side +* vector b in locations b[1], ..., b[n], where n is the order of the +* matrix V. On exit the array x will contain elements of the solution +* vector x in locations x[1], ..., x[n]. Note that the array b will be +* clobbered on exit. */ + +void luf_vt_solve(LUF *luf, double b[/*1+n*/], double x[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + double *vr_piv = luf->vr_piv; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int i, j, k, ptr, end; + double x_i; + for (k = 1; k <= n; k++) + { /* k-th row of U' = j-th column of V */ + /* k-th column of U' = i-th row of V */ + i = pp_inv[k]; + j = qq_ind[k]; + /* compute x[i] = b[j] / u'[k,k], where u'[k,k] = v[i,j]; + * walk through i-th row of matrix V and substitute x[i] into + * other equations */ + if ((x_i = x[i] = b[j] / vr_piv[i]) != 0.0) + { for (end = (ptr = vr_ptr[i]) + vr_len[i]; ptr < end; ptr++) + b[sv_ind[ptr]] -= sv_val[ptr] * x_i; + } + } + return; +} + +/*********************************************************************** +* luf_vt_solve1 - solve system V' * y = e' to cause growth in y +* +* This routine is a special version of luf_vt_solve. It solves the +* system V'* y = e' = e + delta e, where V' is a matrix transposed to +* the matrix V, e is the specified right-hand side vector, and delta e +* is a vector of +1 and -1 chosen to cause growth in the solution +* vector y. +* +* On entry the array e should contain elements of the right-hand side +* vector e in locations e[1], ..., e[n], where n is the order of the +* matrix V. On exit the array y will contain elements of the solution +* vector y in locations y[1], ..., y[n]. Note that the array e will be +* clobbered on exit. */ + +void luf_vt_solve1(LUF *luf, double e[/*1+n*/], double y[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + double *vr_piv = luf->vr_piv; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int i, j, k, ptr, end; + double e_j, y_i; + for (k = 1; k <= n; k++) + { /* k-th row of U' = j-th column of V */ + /* k-th column of U' = i-th row of V */ + i = pp_inv[k]; + j = qq_ind[k]; + /* determine e'[j] = e[j] + delta e[j] */ + e_j = (e[j] >= 0.0 ? e[j] + 1.0 : e[j] - 1.0); + /* compute y[i] = e'[j] / u'[k,k], where u'[k,k] = v[i,j] */ + y_i = y[i] = e_j / vr_piv[i]; + /* walk through i-th row of matrix V and substitute y[i] into + * other equations */ + for (end = (ptr = vr_ptr[i]) + vr_len[i]; ptr < end; ptr++) + e[sv_ind[ptr]] -= sv_val[ptr] * y_i; + } + return; +} + +/*********************************************************************** +* luf_estimate_norm - estimate 1-norm of inv(A) +* +* This routine estimates 1-norm of inv(A) by one step of inverse +* iteration for the small singular vector as described in [1]. This +* involves solving two systems of equations: +* +* A'* y = e, +* +* A * z = y, +* +* where A' is a matrix transposed to A, and e is a vector of +1 and -1 +* chosen to cause growth in y. Then +* +* estimate 1-norm of inv(A) = (1-norm of z) / (1-norm of y) +* +* REFERENCES +* +* 1. G.E.Forsythe, M.A.Malcolm, C.B.Moler. Computer Methods for +* Mathematical Computations. Prentice-Hall, Englewood Cliffs, N.J., +* pp. 30-62 (subroutines DECOMP and SOLVE). */ + +double luf_estimate_norm(LUF *luf, double w1[/*1+n*/], double + w2[/*1+n*/]) +{ int n = luf->n; + double *e = w1; + double *y = w2; + double *z = w1; + int i; + double y_norm, z_norm; + /* y = inv(A') * e = inv(F') * inv(V') * e */ + /* compute y' = inv(V') * e to cause growth in y' */ + for (i = 1; i <= n; i++) + e[i] = 0.0; + luf_vt_solve1(luf, e, y); + /* compute y = inv(F') * y' */ + luf_ft_solve(luf, y); + /* compute 1-norm of y = sum |y[i]| */ + y_norm = 0.0; + for (i = 1; i <= n; i++) + y_norm += (y[i] >= 0.0 ? +y[i] : -y[i]); + /* z = inv(A) * y = inv(V) * inv(F) * y */ + /* compute z' = inv(F) * y */ + luf_f_solve(luf, y); + /* compute z = inv(V) * z' */ + luf_v_solve(luf, y, z); + /* compute 1-norm of z = sum |z[i]| */ + z_norm = 0.0; + for (i = 1; i <= n; i++) + z_norm += (z[i] >= 0.0 ? +z[i] : -z[i]); + /* estimate 1-norm of inv(A) = (1-norm of z) / (1-norm of y) */ + return z_norm / y_norm; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/luf.h b/glpk-5.0/src/bflib/luf.h new file mode 100644 index 0000000..f1eb5e2 --- /dev/null +++ b/glpk-5.0/src/bflib/luf.h @@ -0,0 +1,225 @@ +/* luf.h (sparse LU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef LUF_H +#define LUF_H + +#include "sva.h" + +/*********************************************************************** +* The structure LUF describes sparse LU-factorization. +* +* The LU-factorization has the following format: +* +* A = F * V = P * L * U * Q, (1) +* +* F = P * L * P', (2) +* +* V = P * U * Q, (3) +* +* where A is a given (unsymmetric) square matrix, F and V are matrix +* factors actually computed, L is a lower triangular matrix with unity +* diagonal, U is an upper triangular matrix, P and Q are permutation +* matrices, P' is a matrix transposed to P. All the matrices have the +* same order n. +* +* Matrices F and V are stored in both row- and column-wise sparse +* formats in the associated sparse vector area (SVA). Unity diagonal +* elements of matrix F are not stored. Pivot elements of matrix V +* (which correspond to diagonal elements of matrix U) are stored in +* a separate ordinary array. +* +* Permutation matrices P and Q are stored in ordinary arrays in both +* row- and column-like formats. +* +* Matrices L and U are completely defined by matrices F, V, P, and Q, +* and therefore not stored explicitly. */ + +typedef struct LUF LUF; + +struct LUF +{ /* sparse LU-factorization */ + int n; + /* order of matrices A, F, V, P, Q */ + SVA *sva; + /* associated sparse vector area (SVA) used to store rows and + * columns of matrices F and V; note that different objects may + * share the same SVA */ + /*--------------------------------------------------------------*/ + /* matrix F in row-wise format */ + /* during the factorization process this object is not used */ + int fr_ref; + /* reference number of sparse vector in SVA, which is the first + * row of matrix F */ +#if 0 + 0 + int *fr_ptr = &sva->ptr[fr_ref-1]; + /* fr_ptr[0] is not used; + * fr_ptr[i], 1 <= i <= n, is pointer to i-th row in SVA */ + int *fr_len = &sva->len[fr_ref-1]; + /* fr_len[0] is not used; + * fr_len[i], 1 <= i <= n, is length of i-th row */ +#endif + /*--------------------------------------------------------------*/ + /* matrix F in column-wise format */ + /* during the factorization process this object is constructed + * by columns */ + int fc_ref; + /* reference number of sparse vector in SVA, which is the first + * column of matrix F */ +#if 0 + 0 + int *fc_ptr = &sva->ptr[fc_ref-1]; + /* fc_ptr[0] is not used; + * fc_ptr[j], 1 <= j <= n, is pointer to j-th column in SVA */ + int *fc_len = &sva->len[fc_ref-1]; + /* fc_len[0] is not used; + * fc_len[j], 1 <= j <= n, is length of j-th column */ +#endif + /*--------------------------------------------------------------*/ + /* matrix V in row-wise format */ + int vr_ref; + /* reference number of sparse vector in SVA, which is the first + * row of matrix V */ +#if 0 + 0 + int *vr_ptr = &sva->ptr[vr_ref-1]; + /* vr_ptr[0] is not used; + * vr_ptr[i], 1 <= i <= n, is pointer to i-th row in SVA */ + int *vr_len = &sva->len[vr_ref-1]; + /* vr_len[0] is not used; + * vr_len[i], 1 <= i <= n, is length of i-th row */ + int *vr_cap = &sva->cap[vr_ref-1]; + /* vr_cap[0] is not used; + * vr_cap[i], 1 <= i <= n, is capacity of i-th row */ +#endif + double *vr_piv; /* double vr_piv[1+n]; */ + /* vr_piv[0] is not used; + * vr_piv[i], 1 <= i <= n, is pivot element of i-th row */ + /*--------------------------------------------------------------*/ + /* matrix V in column-wise format */ + /* during the factorization process this object contains only the + * patterns (row indices) of columns of the active submatrix */ + int vc_ref; + /* reference number of sparse vector in SVA, which is the first + * column of matrix V */ +#if 0 + 0 + int *vc_ptr = &sva->ptr[vc_ref-1]; + /* vc_ptr[0] is not used; + * vc_ptr[j], 1 <= j <= n, is pointer to j-th column in SVA */ + int *vc_len = &sva->len[vc_ref-1]; + /* vc_len[0] is not used; + * vc_len[j], 1 <= j <= n, is length of j-th column */ + int *vc_cap = &sva->cap[vc_ref-1]; + /* vc_cap[0] is not used; + * vc_cap[j], 1 <= j <= n, is capacity of j-th column */ +#endif + /*--------------------------------------------------------------*/ + /* matrix P */ + int *pp_ind; /* int pp_ind[1+n]; */ + /* pp_ind[i] = j means that P[i,j] = 1 */ + int *pp_inv; /* int pp_inv[1+n]; */ + /* pp_inv[j] = i means that P[i,j] = 1 */ + /* if i-th row or column of matrix F is i'-th row or column of + * matrix L, or if i-th row of matrix V is i'-th row of matrix U, + * then pp_ind[i] = i' and pp_inv[i'] = i */ + /*--------------------------------------------------------------*/ + /* matrix Q */ + int *qq_ind; /* int qq_ind[1+n]; */ + /* qq_ind[i] = j means that Q[i,j] = 1 */ + int *qq_inv; /* int qq_inv[1+n]; */ + /* qq_inv[j] = i means that Q[i,j] = 1 */ + /* if j-th column of matrix V is j'-th column of matrix U, then + * qq_ind[j'] = j and qq_inv[j] = j' */ +}; + +#define luf_swap_u_rows(i1, i2) \ + do \ + { int j1, j2; \ + j1 = pp_inv[i1], j2 = pp_inv[i2]; \ + pp_ind[j1] = i2, pp_inv[i2] = j1; \ + pp_ind[j2] = i1, pp_inv[i1] = j2; \ + } while (0) +/* swap rows i1 and i2 of matrix U = P'* V * Q' */ + +#define luf_swap_u_cols(j1, j2) \ + do \ + { int i1, i2; \ + i1 = qq_ind[j1], i2 = qq_ind[j2]; \ + qq_ind[j1] = i2, qq_inv[i2] = j1; \ + qq_ind[j2] = i1, qq_inv[i1] = j2; \ + } while (0) +/* swap columns j1 and j2 of matrix U = P'* V * Q' */ + +#define luf_store_v_cols _glp_luf_store_v_cols +int luf_store_v_cols(LUF *luf, int (*col)(void *info, int j, int ind[], + double val[]), void *info, int ind[], double val[]); +/* store matrix V = A in column-wise format */ + +#define luf_check_all _glp_luf_check_all +void luf_check_all(LUF *luf, int k); +/* check LU-factorization before k-th elimination step */ + +#define luf_build_v_rows _glp_luf_build_v_rows +void luf_build_v_rows(LUF *luf, int len[/*1+n*/]); +/* build matrix V in row-wise format */ + +#define luf_build_f_rows _glp_luf_build_f_rows +void luf_build_f_rows(LUF *luf, int len[/*1+n*/]); +/* build matrix F in row-wise format */ + +#define luf_build_v_cols _glp_luf_build_v_cols +void luf_build_v_cols(LUF *luf, int updat, int len[/*1+n*/]); +/* build matrix V in column-wise format */ + +#define luf_check_f_rc _glp_luf_check_f_rc +void luf_check_f_rc(LUF *luf); +/* check rows and columns of matrix F */ + +#define luf_check_v_rc _glp_luf_check_v_rc +void luf_check_v_rc(LUF *luf); +/* check rows and columns of matrix V */ + +#define luf_f_solve _glp_luf_f_solve +void luf_f_solve(LUF *luf, double x[/*1+n*/]); +/* solve system F * x = b */ + +#define luf_ft_solve _glp_luf_ft_solve +void luf_ft_solve(LUF *luf, double x[/*1+n*/]); +/* solve system F' * x = b */ + +#define luf_v_solve _glp_luf_v_solve +void luf_v_solve(LUF *luf, double b[/*1+n*/], double x[/*1+n*/]); +/* solve system V * x = b */ + +#define luf_vt_solve _glp_luf_vt_solve +void luf_vt_solve(LUF *luf, double b[/*1+n*/], double x[/*1+n*/]); +/* solve system V' * x = b */ + +#define luf_vt_solve1 _glp_luf_vt_solve1 +void luf_vt_solve1(LUF *luf, double e[/*1+n*/], double y[/*1+n*/]); +/* solve system V' * y = e' to cause growth in y */ + +#define luf_estimate_norm _glp_luf_estimate_norm +double luf_estimate_norm(LUF *luf, double w1[/*1+n*/], double + w2[/*1+n*/]); +/* estimate 1-norm of inv(A) */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/lufint.c b/glpk-5.0/src/bflib/lufint.c new file mode 100644 index 0000000..6ab90ab --- /dev/null +++ b/glpk-5.0/src/bflib/lufint.c @@ -0,0 +1,180 @@ +/* lufint.c (interface to LU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "lufint.h" + +LUFINT *lufint_create(void) +{ /* create interface to LU-factorization */ + LUFINT *fi; + fi = talloc(1, LUFINT); + fi->n_max = 0; + fi->valid = 0; + fi->sva = NULL; + fi->luf = NULL; + fi->sgf = NULL; + fi->sva_n_max = fi->sva_size = 0; + fi->delta_n0 = fi->delta_n = 0; + fi->sgf_updat = 0; + fi->sgf_piv_tol = 0.10; + fi->sgf_piv_lim = 4; + fi->sgf_suhl = 1; + fi->sgf_eps_tol = DBL_EPSILON; + return fi; +} + +int lufint_factorize(LUFINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info) +{ /* compute LU-factorization of specified matrix A */ + SVA *sva; + LUF *luf; + SGF *sgf; + int k; + xassert(n > 0); + fi->valid = 0; + /* create sparse vector area (SVA), if necessary */ + sva = fi->sva; + if (sva == NULL) + { int sva_n_max = fi->sva_n_max; + int sva_size = fi->sva_size; + if (sva_n_max == 0) + sva_n_max = 4 * n; + if (sva_size == 0) + sva_size = 10 * n; + sva = fi->sva = sva_create_area(sva_n_max, sva_size); + } + /* allocate/reallocate underlying objects, if necessary */ + if (fi->n_max < n) + { int n_max = fi->n_max; + if (n_max == 0) + n_max = fi->n_max = n + fi->delta_n0; + else + n_max = fi->n_max = n + fi->delta_n; + xassert(n_max >= n); + /* allocate/reallocate LU-factorization (LUF) */ + luf = fi->luf; + if (luf == NULL) + { luf = fi->luf = talloc(1, LUF); + memset(luf, 0, sizeof(LUF)); + luf->sva = sva; + } + else + { tfree(luf->vr_piv); + tfree(luf->pp_ind); + tfree(luf->pp_inv); + tfree(luf->qq_ind); + tfree(luf->qq_inv); + } + luf->vr_piv = talloc(1+n_max, double); + luf->pp_ind = talloc(1+n_max, int); + luf->pp_inv = talloc(1+n_max, int); + luf->qq_ind = talloc(1+n_max, int); + luf->qq_inv = talloc(1+n_max, int); + /* allocate/reallocate factorizer workspace (SGF) */ + sgf = fi->sgf; + if (sgf == NULL) + { sgf = fi->sgf = talloc(1, SGF); + memset(sgf, 0, sizeof(SGF)); + sgf->luf = luf; + } + else + { tfree(sgf->rs_head); + tfree(sgf->rs_prev); + tfree(sgf->rs_next); + tfree(sgf->cs_head); + tfree(sgf->cs_prev); + tfree(sgf->cs_next); + tfree(sgf->vr_max); + tfree(sgf->flag); + tfree(sgf->work); + } + sgf->rs_head = talloc(1+n_max, int); + sgf->rs_prev = talloc(1+n_max, int); + sgf->rs_next = talloc(1+n_max, int); + sgf->cs_head = talloc(1+n_max, int); + sgf->cs_prev = talloc(1+n_max, int); + sgf->cs_next = talloc(1+n_max, int); + sgf->vr_max = talloc(1+n_max, double); + sgf->flag = talloc(1+n_max, char); + sgf->work = talloc(1+n_max, double); + } + luf = fi->luf; + sgf = fi->sgf; +#if 1 /* FIXME */ + /* initialize SVA */ + sva->n = 0; + sva->m_ptr = 1; + sva->r_ptr = sva->size + 1; + sva->head = sva->tail = 0; +#endif + /* allocate sparse vectors in SVA */ + luf->n = n; + luf->fr_ref = sva_alloc_vecs(sva, n); + luf->fc_ref = sva_alloc_vecs(sva, n); + luf->vr_ref = sva_alloc_vecs(sva, n); + luf->vc_ref = sva_alloc_vecs(sva, n); + /* store matrix V = A in column-wise format */ + luf_store_v_cols(luf, col, info, sgf->rs_prev, sgf->work); + /* setup factorizer control parameters */ + sgf->updat = fi->sgf_updat; + sgf->piv_tol = fi->sgf_piv_tol; + sgf->piv_lim = fi->sgf_piv_lim; + sgf->suhl = fi->sgf_suhl; + sgf->eps_tol = fi->sgf_eps_tol; + /* compute LU-factorization of specified matrix A */ + k = sgf_factorize(sgf, 1); + if (k == 0) + fi->valid = 1; + return k; +} + +void lufint_delete(LUFINT *fi) +{ /* delete interface to LU-factorization */ + SVA *sva = fi->sva; + LUF *luf = fi->luf; + SGF *sgf = fi->sgf; + if (sva != NULL) + sva_delete_area(sva); + if (luf != NULL) + { tfree(luf->vr_piv); + tfree(luf->pp_ind); + tfree(luf->pp_inv); + tfree(luf->qq_ind); + tfree(luf->qq_inv); + tfree(luf); + } + if (sgf != NULL) + { tfree(sgf->rs_head); + tfree(sgf->rs_prev); + tfree(sgf->rs_next); + tfree(sgf->cs_head); + tfree(sgf->cs_prev); + tfree(sgf->cs_next); + tfree(sgf->vr_max); + tfree(sgf->flag); + tfree(sgf->work); + tfree(sgf); + } + tfree(fi); + return; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/lufint.h b/glpk-5.0/src/bflib/lufint.h new file mode 100644 index 0000000..468f749 --- /dev/null +++ b/glpk-5.0/src/bflib/lufint.h @@ -0,0 +1,71 @@ +/* lufint.h (interface to LU-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef LUFINT_H +#define LUFINT_H + +#include "sgf.h" + +typedef struct LUFINT LUFINT; + +struct LUFINT +{ /* interface to LU-factorization */ + int n_max; + /* maximal value of n (increased automatically) */ + int valid; + /* factorization is valid only if this flag is set */ + SVA *sva; + /* sparse vector area (SVA) */ + LUF *luf; + /* sparse LU-factorization */ + SGF *sgf; + /* sparse Gaussian factorizer workspace */ + /*--------------------------------------------------------------*/ + /* control parameters */ + int sva_n_max, sva_size; + /* parameters passed to sva_create_area */ + int delta_n0, delta_n; + /* if n_max = 0, set n_max = n + delta_n0 + * if n_max < n, set n_max = n + delta_n */ + int sgf_updat; + double sgf_piv_tol; + int sgf_piv_lim; + int sgf_suhl; + double sgf_eps_tol; + /* factorizer control parameters */ +}; + +#define lufint_create _glp_lufint_create +LUFINT *lufint_create(void); +/* create interface to LU-factorization */ + +#define lufint_factorize _glp_lufint_factorize +int lufint_factorize(LUFINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info); +/* compute LU-factorization of specified matrix A */ + +#define lufint_delete _glp_lufint_delete +void lufint_delete(LUFINT *fi); +/* delete interface to LU-factorization */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/scf.c b/glpk-5.0/src/bflib/scf.c new file mode 100644 index 0000000..b4f3a2c --- /dev/null +++ b/glpk-5.0/src/bflib/scf.c @@ -0,0 +1,521 @@ +/* scf.c (sparse updatable Schur-complement-based factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "scf.h" + +/*********************************************************************** +* scf_r0_solve - solve system R0 * x = b or R0'* x = b +* +* This routine solves the system R0 * x = b (if tr is zero) or the +* system R0'* x = b (if tr is non-zero), where R0 is the left factor +* of the initial matrix A0 = R0 * S0. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n0], where n0 is the order of the +* matrix R0. On exit the array x will contain elements of the solution +* vector in the same locations. */ + +void scf_r0_solve(SCF *scf, int tr, double x[/*1+n0*/]) +{ switch (scf->type) + { case 1: + /* A0 = F0 * V0, so R0 = F0 */ + if (!tr) + luf_f_solve(scf->a0.luf, x); + else + luf_ft_solve(scf->a0.luf, x); + break; + case 2: + /* A0 = I * A0, so R0 = I */ + break; + default: + xassert(scf != scf); + } + return; +} + +/*********************************************************************** +* scf_s0_solve - solve system S0 * x = b or S0'* x = b +* +* This routine solves the system S0 * x = b (if tr is zero) or the +* system S0'* x = b (if tr is non-zero), where S0 is the right factor +* of the initial matrix A0 = R0 * S0. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n0], where n0 is the order of the +* matrix S0. On exit the array x will contain elements of the solution +* vector in the same locations. +* +* The routine uses locations [1], ..., [n0] of three working arrays +* w1, w2, and w3. (In case of type = 1 arrays w2 and w3 are not used +* and can be specified as NULL.) */ + +void scf_s0_solve(SCF *scf, int tr, double x[/*1+n0*/], + double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]) +{ int n0 = scf->n0; + switch (scf->type) + { case 1: + /* A0 = F0 * V0, so S0 = V0 */ + if (!tr) + luf_v_solve(scf->a0.luf, x, w1); + else + luf_vt_solve(scf->a0.luf, x, w1); + break; + case 2: + /* A0 = I * A0, so S0 = A0 */ + if (!tr) + btf_a_solve(scf->a0.btf, x, w1, w2, w3); + else + btf_at_solve(scf->a0.btf, x, w1, w2, w3); + break; + default: + xassert(scf != scf); + } + memcpy(&x[1], &w1[1], n0 * sizeof(double)); + return; +} + +/*********************************************************************** +* scf_r_prod - compute product y := y + alpha * R * x +* +* This routine computes the product y := y + alpha * R * x, where +* x is a n0-vector, alpha is a scalar, y is a nn-vector. +* +* Since matrix R is available by rows, the product components are +* computed as inner products: +* +* y[i] = y[i] + alpha * (i-th row of R) * x +* +* for i = 1, 2, ..., nn. */ + +void scf_r_prod(SCF *scf, double y[/*1+nn*/], double a, const double + x[/*1+n0*/]) +{ int nn = scf->nn; + SVA *sva = scf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int rr_ref = scf->rr_ref; + int *rr_ptr = &sva->ptr[rr_ref-1]; + int *rr_len = &sva->len[rr_ref-1]; + int i, ptr, end; + double t; + for (i = 1; i <= nn; i++) + { /* t := (i-th row of R) * x */ + t = 0.0; + for (end = (ptr = rr_ptr[i]) + rr_len[i]; ptr < end; ptr++) + t += sv_val[ptr] * x[sv_ind[ptr]]; + /* y[i] := y[i] + alpha * t */ + y[i] += a * t; + } + return; +} + +/*********************************************************************** +* scf_rt_prod - compute product y := y + alpha * R'* x +* +* This routine computes the product y := y + alpha * R'* x, where +* R' is a matrix transposed to R, x is a nn-vector, alpha is a scalar, +* y is a n0-vector. +* +* Since matrix R is available by rows, the product is computed as a +* linear combination: +* +* y := y + alpha * (R'[1] * x[1] + ... + R'[nn] * x[nn]), +* +* where R'[i] is i-th row of R. */ + +void scf_rt_prod(SCF *scf, double y[/*1+n0*/], double a, const double + x[/*1+nn*/]) +{ int nn = scf->nn; + SVA *sva = scf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int rr_ref = scf->rr_ref; + int *rr_ptr = &sva->ptr[rr_ref-1]; + int *rr_len = &sva->len[rr_ref-1]; + int i, ptr, end; + double t; + for (i = 1; i <= nn; i++) + { if (x[i] == 0.0) + continue; + /* y := y + alpha * R'[i] * x[i] */ + t = a * x[i]; + for (end = (ptr = rr_ptr[i]) + rr_len[i]; ptr < end; ptr++) + y[sv_ind[ptr]] += sv_val[ptr] * t; + } + return; +} + +/*********************************************************************** +* scf_s_prod - compute product y := y + alpha * S * x +* +* This routine computes the product y := y + alpha * S * x, where +* x is a nn-vector, alpha is a scalar, y is a n0 vector. +* +* Since matrix S is available by columns, the product is computed as +* a linear combination: +* +* y := y + alpha * (S[1] * x[1] + ... + S[nn] * x[nn]), +* +* where S[j] is j-th column of S. */ + +void scf_s_prod(SCF *scf, double y[/*1+n0*/], double a, const double + x[/*1+nn*/]) +{ int nn = scf->nn; + SVA *sva = scf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int ss_ref = scf->ss_ref; + int *ss_ptr = &sva->ptr[ss_ref-1]; + int *ss_len = &sva->len[ss_ref-1]; + int j, ptr, end; + double t; + for (j = 1; j <= nn; j++) + { if (x[j] == 0.0) + continue; + /* y := y + alpha * S[j] * x[j] */ + t = a * x[j]; + for (end = (ptr = ss_ptr[j]) + ss_len[j]; ptr < end; ptr++) + y[sv_ind[ptr]] += sv_val[ptr] * t; + } + return; +} + +/*********************************************************************** +* scf_st_prod - compute product y := y + alpha * S'* x +* +* This routine computes the product y := y + alpha * S'* x, where +* S' is a matrix transposed to S, x is a n0-vector, alpha is a scalar, +* y is a nn-vector. +* +* Since matrix S is available by columns, the product components are +* computed as inner products: +* +* y[j] := y[j] + alpha * (j-th column of S) * x +* +* for j = 1, 2, ..., nn. */ + +void scf_st_prod(SCF *scf, double y[/*1+nn*/], double a, const double + x[/*1+n0*/]) +{ int nn = scf->nn; + SVA *sva = scf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int ss_ref = scf->ss_ref; + int *ss_ptr = &sva->ptr[ss_ref-1]; + int *ss_len = &sva->len[ss_ref-1]; + int j, ptr, end; + double t; + for (j = 1; j <= nn; j++) + { /* t := (j-th column of S) * x */ + t = 0.0; + for (end = (ptr = ss_ptr[j]) + ss_len[j]; ptr < end; ptr++) + t += sv_val[ptr] * x[sv_ind[ptr]]; + /* y[j] := y[j] + alpha * t */ + y[j] += a * t; + } + return; +} + +/*********************************************************************** +* scf_a_solve - solve system A * x = b +* +* This routine solves the system A * x = b, where A is the current +* matrix. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix A. On exit the array x will contain elements of the solution +* vector in the same locations. +* +* For details see the program documentation. */ + +void scf_a_solve(SCF *scf, double x[/*1+n*/], + double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/], + double work2[/*1+n*/], double work3[/*1+n*/]) +{ int n = scf->n; + int n0 = scf->n0; + int nn = scf->nn; + int *pp_ind = scf->pp_ind; + int *qq_inv = scf->qq_inv; + int i, ii; + /* (u1, u2) := inv(P) * (b, 0) */ + for (ii = 1; ii <= n0+nn; ii++) + { i = pp_ind[ii]; +#if 1 /* FIXME: currently P = I */ + xassert(i == ii); +#endif + w[ii] = (i <= n ? x[i] : 0.0); + } + /* v1 := inv(R0) * u1 */ + scf_r0_solve(scf, 0, &w[0]); + /* v2 := u2 - R * v1 */ + scf_r_prod(scf, &w[n0], -1.0, &w[0]); + /* w2 := inv(C) * v2 */ + ifu_a_solve(&scf->ifu, &w[n0], work1); + /* w1 := inv(S0) * (v1 - S * w2) */ + scf_s_prod(scf, &w[0], -1.0, &w[n0]); + scf_s0_solve(scf, 0, &w[0], work1, work2, work3); + /* (x, x~) := inv(Q) * (w1, w2); x~ is not needed */ + for (i = 1; i <= n; i++) + x[i] = w[qq_inv[i]]; + return; +} + +/*********************************************************************** +* scf_at_solve - solve system A'* x = b +* +* This routine solves the system A'* x = b, where A' is a matrix +* transposed to the current matrix A. +* +* On entry the array x should contain elements of the right-hand side +* vector b in locations x[1], ..., x[n], where n is the order of the +* matrix A. On exit the array x will contain elements of the solution +* vector in the same locations. +* +* For details see the program documentation. */ + +void scf_at_solve(SCF *scf, double x[/*1+n*/], + double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/], + double work2[/*1+n*/], double work3[/*1+n*/]) +{ int n = scf->n; + int n0 = scf->n0; + int nn = scf->nn; + int *pp_inv = scf->pp_inv; + int *qq_ind = scf->qq_ind; + int i, ii; + /* (u1, u2) := Q * (b, 0) */ + for (ii = 1; ii <= n0+nn; ii++) + { i = qq_ind[ii]; + w[ii] = (i <= n ? x[i] : 0.0); + } + /* v1 := inv(S0') * u1 */ + scf_s0_solve(scf, 1, &w[0], work1, work2, work3); + /* v2 := inv(C') * (u2 - S'* v1) */ + scf_st_prod(scf, &w[n0], -1.0, &w[0]); + ifu_at_solve(&scf->ifu, &w[n0], work1); + /* w2 := v2 */ + /* nop */ + /* w1 := inv(R0') * (v1 - R'* w2) */ + scf_rt_prod(scf, &w[0], -1.0, &w[n0]); + scf_r0_solve(scf, 1, &w[0]); + /* compute (x, x~) := P * (w1, w2); x~ is not needed */ + for (i = 1; i <= n; i++) + { +#if 1 /* FIXME: currently P = I */ + xassert(pp_inv[i] == i); +#endif + x[i] = w[pp_inv[i]]; + } + return; +} + +/*********************************************************************** +* scf_add_r_row - add new row to matrix R +* +* This routine adds new (nn+1)-th row to matrix R, whose elements are +* specified in locations w[1,...,n0]. */ + +void scf_add_r_row(SCF *scf, const double w[/*1+n0*/]) +{ int n0 = scf->n0; + int nn = scf->nn; + SVA *sva = scf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int rr_ref = scf->rr_ref; + int *rr_ptr = &sva->ptr[rr_ref-1]; + int *rr_len = &sva->len[rr_ref-1]; + int j, len, ptr; + xassert(0 <= nn && nn < scf->nn_max); + /* determine length of new row */ + len = 0; + for (j = 1; j <= n0; j++) + { if (w[j] != 0.0) + len++; + } + /* reserve locations for new row in static part of SVA */ + if (len > 0) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_reserve_cap(sva, rr_ref + nn, len); + } + /* store new row in sparse format */ + ptr = rr_ptr[nn+1]; + for (j = 1; j <= n0; j++) + { if (w[j] != 0.0) + { sv_ind[ptr] = j; + sv_val[ptr] = w[j]; + ptr++; + } + } + xassert(ptr - rr_ptr[nn+1] == len); + rr_len[nn+1] = len; +#ifdef GLP_DEBUG + sva_check_area(sva); +#endif + return; +} + +/*********************************************************************** +* scf_add_s_col - add new column to matrix S +* +* This routine adds new (nn+1)-th column to matrix S, whose elements +* are specified in locations v[1,...,n0]. */ + +void scf_add_s_col(SCF *scf, const double v[/*1+n0*/]) +{ int n0 = scf->n0; + int nn = scf->nn; + SVA *sva = scf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int ss_ref = scf->ss_ref; + int *ss_ptr = &sva->ptr[ss_ref-1]; + int *ss_len = &sva->len[ss_ref-1]; + int i, len, ptr; + xassert(0 <= nn && nn < scf->nn_max); + /* determine length of new column */ + len = 0; + for (i = 1; i <= n0; i++) + { if (v[i] != 0.0) + len++; + } + /* reserve locations for new column in static part of SVA */ + if (len > 0) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_reserve_cap(sva, ss_ref + nn, len); + } + /* store new column in sparse format */ + ptr = ss_ptr[nn+1]; + for (i = 1; i <= n0; i++) + { if (v[i] != 0.0) + { sv_ind[ptr] = i; + sv_val[ptr] = v[i]; + ptr++; + } + } + xassert(ptr - ss_ptr[nn+1] == len); + ss_len[nn+1] = len; +#ifdef GLP_DEBUG + sva_check_area(sva); +#endif + return; +} + +/*********************************************************************** +* scf_update_aug - update factorization of augmented matrix +* +* Given factorization of the current augmented matrix: +* +* ( A0 A1 ) ( R0 ) ( S0 S ) +* ( ) = ( ) ( ), +* ( A2 A3 ) ( R I ) ( C ) +* +* this routine computes factorization of the new augmented matrix: +* +* ( A0 | A1 b ) +* ( ---+------ ) ( A0 A1^ ) ( R0 ) ( S0 S^ ) +* ( A2 | A3 f ) = ( ) = ( ) ( ), +* ( | ) ( A2^ A3^ ) ( R^ I ) ( C^ ) +* ( d' | g' h ) +* +* where b and d are specified n0-vectors, f and g are specified +* nn-vectors, and h is a specified scalar. (Note that corresponding +* arrays are clobbered on exit.) +* +* The parameter upd specifies how to update factorization of the Schur +* complement C: +* +* 1 Bartels-Golub updating. +* +* 2 Givens rotations updating. +* +* The working arrays w1, w2, and w3 are used in the same way as in the +* routine scf_s0_solve. +* +* RETURNS +* +* 0 Factorization has been successfully updated. +* +* 1 Updating limit has been reached. +* +* 2 Updating IFU-factorization of matrix C failed. +* +* For details see the program documentation. */ + +int scf_update_aug(SCF *scf, double b[/*1+n0*/], double d[/*1+n0*/], + double f[/*1+nn*/], double g[/*1+nn*/], double h, int upd, + double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]) +{ int n0 = scf->n0; + int k, ret; + double *v, *w, *x, *y, z; + if (scf->nn == scf->nn_max) + { /* updating limit has been reached */ + return 1; + } + /* v := inv(R0) * b */ + scf_r0_solve(scf, 0, (v = b)); + /* w := inv(S0') * d */ + scf_s0_solve(scf, 1, (w = d), w1, w2, w3); + /* x := f - R * v */ + scf_r_prod(scf, (x = f), -1.0, v); + /* y := g - S'* w */ + scf_st_prod(scf, (y = g), -1.0, w); + /* z := h - v'* w */ + z = h; + for (k = 1; k <= n0; k++) + z -= v[k] * w[k]; + /* new R := R with row w added */ + scf_add_r_row(scf, w); + /* new S := S with column v added */ + scf_add_s_col(scf, v); + /* update IFU-factorization of C */ + switch (upd) + { case 1: + ret = ifu_bg_update(&scf->ifu, x, y, z); + break; + case 2: + ret = ifu_gr_update(&scf->ifu, x, y, z); + break; + default: + xassert(upd != upd); + } + if (ret != 0) + { /* updating IFU-factorization failed */ + return 2; + } + /* increase number of additional rows and columns */ + scf->nn++; + /* expand P and Q */ + k = n0 + scf->nn; + scf->pp_ind[k] = scf->pp_inv[k] = k; + scf->qq_ind[k] = scf->qq_inv[k] = k; + /* factorization has been successfully updated */ + return 0; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/scf.h b/glpk-5.0/src/bflib/scf.h new file mode 100644 index 0000000..08f5971 --- /dev/null +++ b/glpk-5.0/src/bflib/scf.h @@ -0,0 +1,209 @@ +/* scf.h (sparse updatable Schur-complement-based factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef SCF_H +#define SCF_H + +#include "btf.h" +#include "ifu.h" +#include "luf.h" + +/*********************************************************************** +* The structure SCF describes sparse updatable factorization based on +* Schur complement. +* +* The SCF-factorization has the following format: +* +* ( A A1~ ) ( A0 A1 ) ( R0 ) ( S0 S ) +* ( ) = P ( ) Q = P ( ) ( ) Q, (1) +* ( A2~ A3~ ) ( A2 A3 ) ( R I ) ( C ) +* +* where: +* +* A is current (unsymmetric) square matrix (not stored); +* +* A1~, A2~, A3~ are some additional matrices (not stored); +* +* A0 is initial (unsymmetric) square matrix (not stored); +* +* A1, A2, A3 are some additional matrices (not stored); +* +* R0 and S0 are matrices that define factorization of the initial +* matrix A0 = R0 * S0 (stored in an invertable form); +* +* R is a matrix defined from R * S0 = A2, so R = A2 * inv(S0) (stored +* in row-wise sparse format); +* +* S is a matrix defined from R0 * S = A1, so S = inv(R0) * A1 (stored +* in column-wise sparse format); +* +* C is Schur complement (to matrix A0) defined from R * S + C = A3, +* so C = A3 - R * S = A3 - A2 * inv(A0) * A1 (stored in an invertable +* form). +* +* P, Q are permutation matrices (stored in both row- and column-like +* formats). */ + +typedef struct SCF SCF; + +struct SCF +{ /* Schur-complement-based factorization */ + int n; + /* order of current matrix A */ + /*--------------------------------------------------------------*/ + /* initial matrix A0 = R0 * S0 of order n0 in invertable form */ + int n0; + /* order of matrix A0 */ + int type; + /* type of factorization used: + * 1 - LU-factorization (R0 = F0, S0 = V0) + * 2 - BT-factorization (R0 = I, S0 = A0) */ + union + { LUF *luf; /* type = 1 */ + BTF *btf; /* type = 2 */ + } a0; + /* factorization of matrix A0 */ + /*--------------------------------------------------------------*/ + /* augmented matrix (A0, A1; A2, A3) of order n0+nn */ + int nn_max; + /* maximal number of additional rows and columns in the augmented + * matrix (this limits the number of updates) */ + int nn; + /* current number of additional rows and columns in the augmented + * matrix, 0 <= nn <= nn_max */ + SVA *sva; + /* associated sparse vector area (SVA) used to store rows of + * matrix R and columns of matrix S */ + /*--------------------------------------------------------------*/ + /* nn*n0-matrix R in row-wise format */ + int rr_ref; + /* reference number of sparse vector in SVA, which is the first + * row of matrix R */ +#if 0 + 0 + int *rr_ptr = &sva->ptr[rr_ref-1]; + /* rr_ptr[0] is not used; + * rr_ptr[i], 1 <= i <= nn, is pointer to i-th row in SVA; + * rr_ptr[nn+1,...,nn_max] are reserved locations */ + int *rr_len = &sva->len[rr_ref-1]; + /* rr_len[0] is not used; + * rr_len[i], 1 <= i <= nn, is length of i-th row; + * rr_len[nn+1,...,nn_max] are reserved locations */ +#endif + /*--------------------------------------------------------------*/ + /* n0*nn-matrix S in column-wise format */ + int ss_ref; + /* reference number of sparse vector in SVA, which is the first + * column of matrix S */ +#if 0 + 0 + int *ss_ptr = &sva->ptr[ss_ref-1]; + /* ss_ptr[0] is not used; + * ss_ptr[j], 1 <= j <= nn, is pointer to j-th column in SVA; + * ss_ptr[nn+1,...,nn_max] are reserved locations */ + int *ss_len = &sva->len[ss_ref-1]; + /* ss_len[0] is not used; + * ss_len[j], 1 <= j <= nn, is length of j-th column; + * ss_len[nn+1,...,nn_max] are reserved locations */ +#endif + /*--------------------------------------------------------------*/ + /* Schur complement C of order nn in invertable form */ + IFU ifu; + /* IFU-factorization of matrix C */ + /*--------------------------------------------------------------*/ + /* permutation matrix P of order n0+nn */ + int *pp_ind; /* int pp_ind[1+n0+nn_max]; */ + /* pp_ind[i] = j means that P[i,j] = 1 */ + int *pp_inv; /* int pp_inv[1+n0+nn_max]; */ + /* pp_inv[j] = i means that P[i,j] = 1 */ + /*--------------------------------------------------------------*/ + /* permutation matrix Q of order n0+nn */ + int *qq_ind; /* int qq_ind[1+n0+nn_max]; */ + /* qq_ind[i] = j means that Q[i,j] = 1 */ + int *qq_inv; /* int qq_inv[1+n0+nn_max]; */ + /* qq_inv[j] = i means that Q[i,j] = 1 */ +}; + +#define scf_swap_q_cols(j1, j2) \ + do \ + { int i1, i2; \ + i1 = qq_inv[j1], i2 = qq_inv[j2]; \ + qq_ind[i1] = j2, qq_inv[j2] = i1; \ + qq_ind[i2] = j1, qq_inv[j1] = i2; \ + } while (0) +/* swap columns j1 and j2 of permutation matrix Q */ + +#define scf_r0_solve _glp_scf_r0_solve +void scf_r0_solve(SCF *scf, int tr, double x[/*1+n0*/]); +/* solve system R0 * x = b or R0'* x = b */ + +#define scf_s0_solve _glp_scf_s0_solve +void scf_s0_solve(SCF *scf, int tr, double x[/*1+n0*/], + double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]); +/* solve system S0 * x = b or S0'* x = b */ + +#define scf_r_prod _glp_scf_r_prod +void scf_r_prod(SCF *scf, double y[/*1+nn*/], double a, const double + x[/*1+n0*/]); +/* compute product y := y + alpha * R * x */ + +#define scf_rt_prod _glp_scf_rt_prod +void scf_rt_prod(SCF *scf, double y[/*1+n0*/], double a, const double + x[/*1+nn*/]); +/* compute product y := y + alpha * R'* x */ + +#define scf_s_prod _glp_scf_s_prod +void scf_s_prod(SCF *scf, double y[/*1+n0*/], double a, const double + x[/*1+nn*/]); +/* compute product y := y + alpha * S * x */ + +#define scf_st_prod _glp_scf_st_prod +void scf_st_prod(SCF *scf, double y[/*1+nn*/], double a, const double + x[/*1+n0*/]); +/* compute product y := y + alpha * S'* x */ + +#define scf_a_solve _glp_scf_a_solve +void scf_a_solve(SCF *scf, double x[/*1+n*/], + double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/], + double work2[/*1+n*/], double work3[/*1+n*/]); +/* solve system A * x = b */ + +#define scf_at_solve _glp_scf_at_solve +void scf_at_solve(SCF *scf, double x[/*1+n*/], + double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/], + double work2[/*1+n*/], double work3[/*1+n*/]); +/* solve system A'* x = b */ + +#define scf_add_r_row _glp_scf_add_r_row +void scf_add_r_row(SCF *scf, const double w[/*1+n0*/]); +/* add new row to matrix R */ + +#define scf_add_s_col _glp_scf_add_s_col +void scf_add_s_col(SCF *scf, const double v[/*1+n0*/]); +/* add new column to matrix S */ + +#define scf_update_aug _glp_scf_update_aug +int scf_update_aug(SCF *scf, double b[/*1+n0*/], double d[/*1+n0*/], + double f[/*1+nn*/], double g[/*1+nn*/], double h, int upd, + double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]); +/* update factorization of augmented matrix */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/scfint.c b/glpk-5.0/src/bflib/scfint.c new file mode 100644 index 0000000..8fb39c7 --- /dev/null +++ b/glpk-5.0/src/bflib/scfint.c @@ -0,0 +1,253 @@ +/* scfint.c (interface to Schur-complement-based factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "scfint.h" + +SCFINT *scfint_create(int type) +{ /* create interface to SC-factorization */ + SCFINT *fi; + fi = talloc(1, SCFINT); + memset(fi, 0, sizeof(SCFINT)); + switch ((fi->scf.type = type)) + { case 1: + fi->u.lufi = lufint_create(); + break; + case 2: + fi->u.btfi = btfint_create(); + break; + default: + xassert(type != type); + } + return fi; +} + +int scfint_factorize(SCFINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info) +{ /* compute SC-factorization of specified matrix A */ + int nn_max, old_n0_max, n0_max, k, ret; + xassert(n > 0); + fi->valid = 0; + /* get required value of nn_max */ + nn_max = fi->nn_max; + if (nn_max == 0) + nn_max = 100; + xassert(nn_max > 0); + /* compute factorization of specified matrix A */ + switch (fi->scf.type) + { case 1: + old_n0_max = fi->u.lufi->n_max; + fi->u.lufi->sva_n_max = 4 * n + 2 * nn_max; + ret = lufint_factorize(fi->u.lufi, n, col, info); + n0_max = fi->u.lufi->n_max; + fi->scf.sva = fi->u.lufi->sva; + fi->scf.a0.luf = fi->u.lufi->luf; + break; + case 2: + old_n0_max = fi->u.btfi->n_max; + fi->u.btfi->sva_n_max = 6 * n + 2 * nn_max; + ret = btfint_factorize(fi->u.btfi, n, col, info); + n0_max = fi->u.btfi->n_max; + fi->scf.sva = fi->u.btfi->sva; + fi->scf.a0.btf = fi->u.btfi->btf; + break; + default: + xassert(fi != fi); + } + /* allocate/reallocate arrays, if necessary */ + if (old_n0_max < n0_max) + { if (fi->w1 != NULL) + tfree(fi->w1); + if (fi->w2 != NULL) + tfree(fi->w2); + if (fi->w3 != NULL) + tfree(fi->w3); + fi->w1 = talloc(1+n0_max, double); + fi->w2 = talloc(1+n0_max, double); + fi->w3 = talloc(1+n0_max, double); + } + if (fi->scf.nn_max != nn_max) + { if (fi->scf.ifu.f != NULL) + tfree(fi->scf.ifu.f); + if (fi->scf.ifu.u != NULL) + tfree(fi->scf.ifu.u); + fi->scf.ifu.f = talloc(nn_max * nn_max, double); + fi->scf.ifu.u = talloc(nn_max * nn_max, double); + } + if (old_n0_max < n0_max || fi->scf.nn_max != nn_max) + { if (fi->scf.pp_ind != NULL) + tfree(fi->scf.pp_ind); + if (fi->scf.pp_inv != NULL) + tfree(fi->scf.pp_inv); + if (fi->scf.qq_ind != NULL) + tfree(fi->scf.qq_ind); + if (fi->scf.qq_inv != NULL) + tfree(fi->scf.qq_inv); + if (fi->w4 != NULL) + tfree(fi->w4); + if (fi->w5 != NULL) + tfree(fi->w5); + fi->scf.pp_ind = talloc(1+n0_max+nn_max, int); + fi->scf.pp_inv = talloc(1+n0_max+nn_max, int); + fi->scf.qq_ind = talloc(1+n0_max+nn_max, int); + fi->scf.qq_inv = talloc(1+n0_max+nn_max, int); + fi->w4 = talloc(1+n0_max+nn_max, double); + fi->w5 = talloc(1+n0_max+nn_max, double); + } + /* initialize SC-factorization */ + fi->scf.n = n; + fi->scf.n0 = n; + fi->scf.nn_max = nn_max; + fi->scf.nn = 0; + fi->scf.rr_ref = sva_alloc_vecs(fi->scf.sva, nn_max); + fi->scf.ss_ref = sva_alloc_vecs(fi->scf.sva, nn_max); + fi->scf.ifu.n_max = nn_max; + fi->scf.ifu.n = 0; + for (k = 1; k <= n; k++) + { fi->scf.pp_ind[k] = k; + fi->scf.pp_inv[k] = k; + fi->scf.qq_ind[k] = k; + fi->scf.qq_inv[k] = k; + } + /* set validation flag */ + if (ret == 0) + fi->valid = 1; + return ret; +} + +int scfint_update(SCFINT *fi, int upd, int j, int len, const int ind[], + const double val[]) +{ /* update SC-factorization after replacing j-th column of A */ + int n = fi->scf.n; + int n0 = fi->scf.n0; + int nn = fi->scf.nn; + int *pp_ind = fi->scf.pp_ind; + int *qq_ind = fi->scf.qq_ind; + int *qq_inv = fi->scf.qq_inv; + double *bf = fi->w4; + double *dg = fi->w5; + int k, t, ret; + xassert(fi->valid); + xassert(0 <= n && n <= n0+nn); + /* (b, f) := inv(P) * (beta, 0) */ + for (k = 1; k <= n0+nn; k++) + bf[k] = 0.0; + for (t = 1; t <= len; t++) + { k = ind[t]; + xassert(1 <= k && k <= n); +#if 1 /* FIXME: currently P = I */ + xassert(pp_ind[k] == k); +#endif + xassert(bf[k] == 0.0); + xassert(val[t] != 0.0); + bf[k] = val[t]; + } + /* (d, g) := Q * (cj, 0) */ + for (k = 1; k <= n0+nn; k++) + dg[k] = 0.0; + xassert(1 <= j && j <= n); + dg[fi->scf.qq_inv[j]] = 1; + /* update factorization of augmented matrix */ + ret = scf_update_aug(&fi->scf, &bf[0], &dg[0], &bf[n0], &dg[n0], + 0.0, upd, fi->w1, fi->w2, fi->w3); + if (ret == 0) + { /* swap j-th and last columns of new matrix Q */ + scf_swap_q_cols(j, n0+nn+1); + } + else + { /* updating failed */ + fi->valid = 0; + } + return ret; +} + +void scfint_ftran(SCFINT *fi, double x[]) +{ /* solve system A * x = b */ + xassert(fi->valid); + scf_a_solve(&fi->scf, x, fi->w4, fi->w5, fi->w1, fi->w2); + return; +} + +void scfint_btran(SCFINT *fi, double x[]) +{ /* solve system A'* x = b */ + xassert(fi->valid); + scf_at_solve(&fi->scf, x, fi->w4, fi->w5, fi->w1, fi->w2); + return; +} + +double scfint_estimate(SCFINT *fi) +{ /* estimate 1-norm of inv(A) */ + double norm; + xassert(fi->valid); + xassert(fi->scf.n == fi->scf.n0); + switch (fi->scf.type) + { case 1: + norm = luf_estimate_norm(fi->scf.a0.luf, fi->w1, fi->w2); + break; + case 2: + norm = btf_estimate_norm(fi->scf.a0.btf, fi->w1, fi->w2, + fi->w3, fi->w4); + break; + default: + xassert(fi != fi); + } + return norm; +} + +void scfint_delete(SCFINT *fi) +{ /* delete interface to SC-factorization */ + switch (fi->scf.type) + { case 1: + lufint_delete(fi->u.lufi); + break; + case 2: + btfint_delete(fi->u.btfi); + break; + default: + xassert(fi != fi); + } + if (fi->scf.ifu.f != NULL) + tfree(fi->scf.ifu.f); + if (fi->scf.ifu.u != NULL) + tfree(fi->scf.ifu.u); + if (fi->scf.pp_ind != NULL) + tfree(fi->scf.pp_ind); + if (fi->scf.pp_inv != NULL) + tfree(fi->scf.pp_inv); + if (fi->scf.qq_ind != NULL) + tfree(fi->scf.qq_ind); + if (fi->scf.qq_inv != NULL) + tfree(fi->scf.qq_inv); + if (fi->w1 != NULL) + tfree(fi->w1); + if (fi->w2 != NULL) + tfree(fi->w2); + if (fi->w3 != NULL) + tfree(fi->w3); + if (fi->w4 != NULL) + tfree(fi->w4); + if (fi->w5 != NULL) + tfree(fi->w5); + tfree(fi); + return; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/scfint.h b/glpk-5.0/src/bflib/scfint.h new file mode 100644 index 0000000..dc25d96 --- /dev/null +++ b/glpk-5.0/src/bflib/scfint.h @@ -0,0 +1,87 @@ +/* scfint.h (interface to Schur-complement-based factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2013-2014 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef SCFINT_H +#define SCFINT_H + +#include "scf.h" +#include "lufint.h" +#include "btfint.h" + +typedef struct SCFINT SCFINT; + +struct SCFINT +{ /* interface to SC-factorization */ + int valid; + /* factorization is valid only if this flag is set */ + SCF scf; + /* Schur-complement based factorization */ + union + { LUFINT *lufi; /* scf.type = 1 */ + BTFINT *btfi; /* scf.type = 2 */ + } u; + /* interface to factorize initial matrix A0 */ + /*--------------------------------------------------------------*/ + /* working arrays */ + double *w1; /* double w1[1+n0_max]; */ + double *w2; /* double w2[1+n0_max]; */ + double *w3; /* double w3[1+n0_max]; */ + double *w4; /* double w4[1+n0_max+nn_max]; */ + double *w5; /* double w5[1+n0_max+nn_max]; */ + /*--------------------------------------------------------------*/ + /* control parameters */ + int nn_max; + /* required maximal number of updates */ +}; + +#define scfint_create _glp_scfint_create +SCFINT *scfint_create(int type); +/* create interface to SC-factorization */ + +#define scfint_factorize _glp_scfint_factorize +int scfint_factorize(SCFINT *fi, int n, int (*col)(void *info, int j, + int ind[], double val[]), void *info); +/* compute SC-factorization of specified matrix A */ + +#define scfint_update _glp_scfint_update +int scfint_update(SCFINT *fi, int upd, int j, int len, const int ind[], + const double val[]); +/* update SC-factorization after replacing j-th column of A */ + +#define scfint_ftran _glp_scfint_ftran +void scfint_ftran(SCFINT *fi, double x[]); +/* solve system A * x = b */ + +#define scfint_btran _glp_scfint_btran +void scfint_btran(SCFINT *fi, double x[]); +/* solve system A'* x = b */ + +#define scfint_estimate _glp_scfint_estimate +double scfint_estimate(SCFINT *fi); +/* estimate 1-norm of inv(A) */ + +#define scfint_delete _glp_scfint_delete +void scfint_delete(SCFINT *fi); +/* delete interface to SC-factorization */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/sgf.c b/glpk-5.0/src/bflib/sgf.c new file mode 100644 index 0000000..6e9d801 --- /dev/null +++ b/glpk-5.0/src/bflib/sgf.c @@ -0,0 +1,1441 @@ +/* sgf.c (sparse Gaussian factorizer) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2015 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "sgf.h" + +/*********************************************************************** +* sgf_reduce_nuc - initial reordering to minimize nucleus size +* +* On entry to this routine it is assumed that V = A and F = P = Q = I, +* where A is the original matrix to be factorized. It is also assumed +* that matrix V = A is stored in both row- and column-wise formats. +* +* This routine performs (implicit) non-symmetric permutations of rows +* and columns of matrix U = P'* V * Q' to reduce it to the form: +* +* 1 k1 k2 n +* 1 x x x x x x x x x x +* . x x x x x x x x x +* . . x x x x x x x x +* k1 . . . * * * * x x x +* . . . * * * * x x x +* . . . * * * * x x x +* k2 . . . * * * * x x x +* . . . . . . . x x x +* . . . . . . . . x x +* n . . . . . . . . . x +* +* where non-zeros in rows and columns k1, k1+1, ..., k2 constitute so +* called nucleus ('*'), whose size is minimized by the routine. +* +* The numbers k1 and k2 are returned by the routine on exit. Usually, +* if the nucleus exists, 1 <= k1 < k2 <= n. However, if the resultant +* matrix U is upper triangular (has no nucleus), k1 = n+1 and k2 = n. +* +* Note that the routines sgf_choose_pivot and sgf_eliminate perform +* exactly the same transformations (by processing row and columns +* singletons), so preliminary minimization of the nucleus may not be +* used. However, processing row and column singletons by the routines +* sgf_minimize_nuc and sgf_singl_phase is more efficient. */ + +#if 1 /* 21/II-2016 */ +/* Normally this routine returns zero. If the matrix is structurally +* singular, the routine returns non-zero. */ +#endif + +int sgf_reduce_nuc(LUF *luf, int *k1_, int *k2_, int cnt[/*1+n*/], + int list[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *pp_ind = luf->pp_ind; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int *qq_inv = luf->qq_inv; + int i, ii, j, jj, k1, k2, ns, ptr, end; + /* initial nucleus is U = V = A */ + k1 = 1, k2 = n; + /*--------------------------------------------------------------*/ + /* process column singletons */ + /*--------------------------------------------------------------*/ + /* determine initial counts of columns of V and initialize list + * of active column singletons */ + ns = 0; /* number of active column singletons */ + for (j = 1; j <= n; j++) + { if ((cnt[j] = vc_len[j]) == 1) + list[++ns] = j; + } + /* process active column singletons */ + while (ns > 0) + { /* column singleton is in j-th column of V */ + j = list[ns--]; +#if 1 /* 21/II-2016 */ + if (cnt[j] == 0) + { /* j-th column in the current nucleus is actually empty */ + /* this happened because on a previous step in the nucleus + * there were two or more identical column singletons (that + * means structural singularity), so removing one of them + * from the nucleus made other columns empty */ + return 1; + } +#endif + /* find i-th row of V containing column singleton */ + ptr = vc_ptr[j]; + end = ptr + vc_len[j]; + for (; pp_ind[i = sv_ind[ptr]] < k1; ptr++) + /* nop */; + xassert(ptr < end); + /* permute rows and columns of U to move column singleton to + * position u[k1,k1] */ + ii = pp_ind[i]; + luf_swap_u_rows(k1, ii); + jj = qq_inv[j]; + luf_swap_u_cols(k1, jj); + /* nucleus size decreased */ + k1++; + /* walk thru i-th row of V and decrease column counts; this + * may cause new column singletons to appear */ + ptr = vr_ptr[i]; + end = ptr + vr_len[i]; + for (; ptr < end; ptr++) + { if (--(cnt[j = sv_ind[ptr]]) == 1) + list[++ns] = j; + } + } + /* nucleus begins at k1-th row/column of U */ + if (k1 > n) + { /* U is upper triangular; no nucleus exist */ + goto done; + } + /*--------------------------------------------------------------*/ + /* process row singletons */ + /*--------------------------------------------------------------*/ + /* determine initial counts of rows of V and initialize list of + * active row singletons */ + ns = 0; /* number of active row singletons */ + for (i = 1; i <= n; i++) + { if (pp_ind[i] < k1) + { /* corresponding row of U is above its k1-th row; set its + * count to zero to prevent including it in active list */ + cnt[i] = 0; + } + else if ((cnt[i] = vr_len[i]) == 1) + list[++ns] = i; + } + /* process active row singletons */ + while (ns > 0) + { /* row singleton is in i-th row of V */ + i = list[ns--]; +#if 1 /* 21/II-2016 */ + if (cnt[i] == 0) + { /* i-th row in the current nucleus is actually empty */ + /* (see comments above for similar case of empty column) */ + return 2; + } +#endif + /* find j-th column of V containing row singleton */ + ptr = vr_ptr[i]; + end = ptr + vr_len[i]; + for (; qq_inv[j = sv_ind[ptr]] > k2; ptr++) + /* nop */; + xassert(ptr < end); + /* permute rows and columns of U to move row singleton to + * position u[k2,k2] */ + ii = pp_ind[i]; + luf_swap_u_rows(k2, ii); + jj = qq_inv[j]; + luf_swap_u_cols(k2, jj); + /* nucleus size decreased */ + k2--; + /* walk thru j-th column of V and decrease row counts; this + * may cause new row singletons to appear */ + ptr = vc_ptr[j]; + end = ptr + vc_len[j]; + for (; ptr < end; ptr++) + { if (--(cnt[i = sv_ind[ptr]]) == 1) + list[++ns] = i; + } + } + /* nucleus ends at k2-th row/column of U */ + xassert(k1 < k2); +done: *k1_ = k1, *k2_ = k2; + return 0; +} + +/*********************************************************************** +* sgf_singl_phase - compute LU-factorization (singleton phase) +* +* It is assumed that on entry to the routine L = P'* F * P = F = I +* and matrix U = P'* V * Q' has the following structure (provided by +* the routine sgf_reduce_nuc): +* +* 1 k1 k2 n +* 1 a a a b b b b c c c +* . a a b b b b c c c +* . . a b b b b c c c +* k1 . . . * * * * d d d +* . . . * * * * d d d +* . . . * * * * d d d +* k2 . . . * * * * d d d +* . . . . . . . e e e +* . . . . . . . . e e +* n . . . . . . . . . e +* +* First, the routine performs (implicit) symmetric permutations of +* rows and columns of matrix U to place them in the following order: +* +* 1, 2, ..., k1-1; n, n-1, ..., k2+1; k1, k1+1, ..., k2 +* +* This changes the structure of matrix U as follows: +* +* 1 k1 k2' n +* 1 a a a c c c b b b b +* . a a c c c b b b b +* . . a c c c b b b b +* k1 . . . e . . . . . . +* . . . e e . . . . . +* . . . e e e . . . . +* k2'. . . d d d * * * * +* . . . d d d * * * * +* . . . d d d * * * * +* n . . . d d d * * * * +* +* where k2' = n - k2 + k1. +* +* Then the routine performs elementary gaussian transformations to +* eliminate subdiagonal elements in columns k1, ..., k2'-1 of U. The +* effect is the same as if the routine sgf_eliminate would be called +* for k = 1, ..., k2'-1 using diagonal elements u[k,k] as pivots. +* +* After elimination matrices L and U becomes the following: +* +* 1 k1 k2' n 1 k1 k2' n +* 1 1 . . . . . . . . . 1 a a a c c c b b b b +* . 1 . . . . . . . . . a a c c c b b b b +* . . 1 . . . . . . . . . a c c c b b b b +* k1 . . . 1 . . . . . . k1 . . . e . . . . . . +* . . . e'1 . . . . . . . . . e . . . . . +* . . . e'e'1 . . . . . . . . . e . . . . +* k2'. . . d'd'd'1 . . . k2'. . . . . . * * * * +* . . . d'd'd'. 1 . . . . . . . . * * * * +* . . . d'd'd'. . 1 . . . . . . . * * * * +* n . . . d'd'd'. . . 1 n . . . . . . * * * * +* +* matrix L matrix U +* +* where columns k1, ..., k2'-1 of L consist of subdiagonal elements +* of initial matrix U divided by pivots u[k,k]. +* +* On exit the routine returns k2', the elimination step number, from +* which computing of the factorization should be continued. Note that +* k2' = n+1 means that matrix U is already upper triangular. */ + +int sgf_singl_phase(LUF *luf, int k1, int k2, int updat, + int ind[/*1+n*/], double val[/*1+n*/]) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fc_ref = luf->fc_ref; + int *fc_ptr = &sva->ptr[fc_ref-1]; + int *fc_len = &sva->len[fc_ref-1]; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + double *vr_piv = luf->vr_piv; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *pp_ind = luf->pp_ind; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int *qq_inv = luf->qq_inv; + int i, j, k, ptr, ptr1, end, len; + double piv; + /* (see routine sgf_reduce_nuc) */ + xassert((1 <= k1 && k1 < k2 && k2 <= n) + || (k1 == n+1 && k2 == n)); + /* perform symmetric permutations of rows/columns of U */ + for (k = k1; k <= k2; k++) + pp_ind[pp_inv[k]] = qq_inv[qq_ind[k]] = k - k2 + n; + for (k = k2+1; k <= n; k++) + pp_ind[pp_inv[k]] = qq_inv[qq_ind[k]] = n - k + k1; + for (k = 1; k <= n; k++) + pp_inv[pp_ind[k]] = qq_ind[qq_inv[k]] = k; + /* determine k2' */ + k2 = n - k2 + k1; + /* process rows and columns of V corresponding to rows and + * columns 1, ..., k1-1 of U */ + for (k = 1; k < k1; k++) + { /* k-th row of U = i-th row of V */ + i = pp_inv[k]; + /* find pivot u[k,k] = v[i,j] in i-th row of V */ + ptr = vr_ptr[i]; + end = ptr + vr_len[i]; + for (; qq_inv[sv_ind[ptr]] != k; ptr++) + /* nop */; + xassert(ptr < end); + /* store pivot */ + vr_piv[i] = sv_val[ptr]; + /* and remove it from i-th row of V */ + sv_ind[ptr] = sv_ind[end-1]; + sv_val[ptr] = sv_val[end-1]; + vr_len[i]--; + /* clear column of V corresponding to k-th column of U */ + vc_len[qq_ind[k]] = 0; + } + /* clear rows of V corresponding to rows k1, ..., k2'-1 of U */ + for (k = k1; k < k2; k++) + vr_len[pp_inv[k]] = 0; + /* process rows and columns of V corresponding to rows and + * columns k2', ..., n of U */ + for (k = k2; k <= n; k++) + { /* k-th row of U = i-th row of V */ + i = pp_inv[k]; + /* remove elements from i-th row of V that correspond to + * elements u[k,k1], ..., u[k,k2'-1] */ + ptr = ptr1 = vr_ptr[i]; + end = ptr + vr_len[i]; + for (; ptr < end; ptr++) + { if (qq_inv[sv_ind[ptr]] >= k2) + { sv_ind[ptr1] = sv_ind[ptr]; + sv_val[ptr1] = sv_val[ptr]; + ptr1++; + } + } + vr_len[i] = ptr1 - vr_ptr[i]; + /* k-th column of U = j-th column of V */ + j = qq_ind[k]; + /* remove elements from j-th column of V that correspond to + * elements u[1,k], ..., u[k1-1,k] */ + ptr = ptr1 = vc_ptr[j]; + end = ptr + vc_len[j]; + for (; ptr < end; ptr++) + { if (pp_ind[sv_ind[ptr]] >= k2) + /* element value is not needed in this case */ + sv_ind[ptr1++] = sv_ind[ptr]; + } + vc_len[j] = ptr1 - vc_ptr[j]; + } + /* process columns of V corresponding to columns k1, ..., k2'-1 + * of U, build columns of F */ + for (k = k1; k < k2; k++) + { /* k-th column of U = j-th column of V */ + j = qq_ind[k]; + /* remove elements from j-th column of V that correspond to + * pivot (diagonal) element u[k,k] and subdiagonal elements + * u[k+1,k], ..., u[n,k]; subdiagonal elements are stored for + * further addition to matrix F */ + len = 0; + piv = 0.0; + ptr = vc_ptr[j]; + end = ptr + vc_len[j]; + for (; ptr < end; ptr++) + { i = sv_ind[ptr]; /* v[i,j] */ + if (pp_ind[i] == k) + { /* store pivot v[i,j] = u[k,k] */ + piv = vr_piv[i] = sv_val[ptr]; + } + else if (pp_ind[i] > k) + { /* store subdiagonal element v[i,j] = u[i',k] */ + len++; + ind[len] = i; + val[len] = sv_val[ptr]; + } + } + /* clear j-th column of V = k-th column of U */ + vc_len[j] = 0; + /* build k-th column of L = j-th column of F */ + j = pp_inv[k]; + xassert(piv != 0.0); + if (len > 0) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_reserve_cap(sva, fc_ref-1+j, len); + for (ptr = fc_ptr[j], ptr1 = 1; ptr1 <= len; ptr++, ptr1++) + { sv_ind[ptr] = ind[ptr1]; + sv_val[ptr] = val[ptr1] / piv; + } + fc_len[j] = len; + } + } + /* if it is not planned to update matrix V, relocate all its + * non-active rows corresponding to rows 1, ..., k2'-1 of U to + * the right (static) part of SVA */ + if (!updat) + { for (k = 1; k < k2; k++) + { i = pp_inv[k]; + len = vr_len[i]; + if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_make_static(sva, vr_ref-1+i); + } + } + /* elimination steps 1, ..., k2'-1 have been performed */ + return k2; +} + +/*********************************************************************** +* sgf_choose_pivot - choose pivot element v[p,q] +* +* This routine chooses pivot element v[p,q], k <= p, q <= n, in the +* active submatrix of matrix V = P * U * Q, where k is the number of +* current elimination step, 1 <= k <= n. +* +* It is assumed that on entry to the routine matrix U = P'* V * Q' has +* the following partially triangularized form: +* +* 1 k n +* 1 x x x x x x x x x x +* . x x x x x x x x x +* . . x x x x x x x x +* . . . x x x x x x x +* k . . . . * * * * * * +* . . . . * * * * * * +* . . . . * * * * * * +* . . . . * * * * * * +* . . . . * * * * * * +* n . . . . * * * * * * +* +* where rows and columns k, k+1, ..., n belong to the active submatrix +* (its elements are marked by '*'). +* +* Since the matrix U is not stored, the routine works with the matrix +* V = P * U * Q. It is assumed that the row-wise representation +* corresponds to the matrix V, but the column-wise representation +* corresponds to the active submatrix of the matrix V, i.e. elements, +* which are not in the active submatrix, are not included in column +* vectors. It is also assumed that each active row of the matrix V is +* in the set R[len], where len is the number of non-zeros in the row, +* and each active column of the matrix V is in the set C[len], where +* len is the number of non-zeros in the column (in the latter case +* only elements of the active submatrix are counted; such elements are +* marked by '*' on the figure above). +* +* For the reason of numerical stability the routine applies so called +* threshold pivoting proposed by J.Reid. It is assumed that an element +* v[i,j] can be selected as a pivot candidate if it is not very small +* (in magnitude) among other elements in the same row, i.e. if it +* satisfies to the stability condition |v[i,j]| >= tol * max|v[i,*]|, +* where 0 < tol < 1 is a given tolerance. +* +* In order to keep sparsity of the matrix V the routine uses Markowitz +* strategy, trying to choose such element v[p,q], which satisfies to +* the stability condition (see above) and has smallest Markowitz cost +* (nr[p]-1) * (nc[q]-1), where nr[p] and nc[q] are, resp., numbers of +* non-zeros in p-th row and q-th column of the active submatrix. +* +* In order to reduce the search, i.e. not to walk through all elements +* of the active submatrix, the routine uses a technique proposed by +* I.Duff. This technique is based on using the sets R[len] and C[len] +* of active rows and columns. +* +* If the pivot element v[p,q] has been chosen, the routine stores its +* indices to locations *p and *q and returns zero. Otherwise, non-zero +* is returned. */ + +int sgf_choose_pivot(SGF *sgf, int *p_, int *q_) +{ LUF *luf = sgf->luf; + int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *rs_head = sgf->rs_head; + int *rs_next = sgf->rs_next; + int *cs_head = sgf->cs_head; + int *cs_prev = sgf->cs_prev; + int *cs_next = sgf->cs_next; + double *vr_max = sgf->vr_max; + double piv_tol = sgf->piv_tol; + int piv_lim = sgf->piv_lim; + int suhl = sgf->suhl; + int i, i_ptr, i_end, j, j_ptr, j_end, len, min_i, min_j, min_len, + ncand, next_j, p, q; + double best, big, cost, temp; + /* no pivot candidate has been chosen so far */ + p = q = 0, best = DBL_MAX, ncand = 0; + /* if the active submatrix contains a column having the only + * non-zero element (column singleton), choose it as the pivot */ + j = cs_head[1]; + if (j != 0) + { xassert(vc_len[j] == 1); + p = sv_ind[vc_ptr[j]], q = j; + goto done; + } + /* if the active submatrix contains a row having the only + * non-zero element (row singleton), choose it as the pivot */ + i = rs_head[1]; + if (i != 0) + { xassert(vr_len[i] == 1); + p = i, q = sv_ind[vr_ptr[i]]; + goto done; + } + /* the active submatrix contains no singletons; walk thru its + * other non-empty rows and columns */ + for (len = 2; len <= n; len++) + { /* consider active columns containing len non-zeros */ + for (j = cs_head[len]; j != 0; j = next_j) + { /* save the number of next column of the same length */ + next_j = cs_next[j]; + /* find an element in j-th column, which is placed in the + * row with minimal number of non-zeros and satisfies to + * the stability condition (such element may not exist) */ + min_i = min_j = 0, min_len = INT_MAX; + for (j_end = (j_ptr = vc_ptr[j]) + vc_len[j]; + j_ptr < j_end; j_ptr++) + { /* get row index of v[i,j] */ + i = sv_ind[j_ptr]; + /* if i-th row is not shorter, skip v[i,j] */ + if (vr_len[i] >= min_len) + continue; + /* big := max|v[i,*]| */ + if ((big = vr_max[i]) < 0.0) + { /* largest magnitude is unknown; compute it */ + for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + i_ptr < i_end; i_ptr++) + { if ((temp = sv_val[i_ptr]) < 0.0) + temp = -temp; + if (big < temp) + big = temp; + } + xassert(big > 0.0); + vr_max[i] = big; + } + /* find v[i,j] in i-th row */ + for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + sv_ind[i_ptr] != j; i_ptr++) + /* nop */; + xassert(i_ptr < i_end); + /* if |v[i,j]| < piv_tol * max|v[i,*]|, skip v[i,j] */ + if ((temp = sv_val[i_ptr]) < 0.0) + temp = -temp; + if (temp < piv_tol * big) + continue; + /* v[i,j] is a better candidate */ + min_i = i, min_j = j, min_len = vr_len[i]; + /* if Markowitz cost of v[i,j] is not greater than + * (len-1)**2, v[i,j] can be chosen as the pivot right + * now; this heuristic reduces the search and works well + * in many cases */ + if (min_len <= len) + { p = min_i, q = min_j; + goto done; + } + } + /* j-th column has been scanned */ + if (min_i != 0) + { /* element v[min_i,min_j] is a next pivot candidate */ + ncand++; + /* compute its Markowitz cost */ + cost = (double)(min_len - 1) * (double)(len - 1); + /* if this element is better, choose it as the pivot */ + if (cost < best) + p = min_i, q = min_j, best = cost; + /* if piv_lim candidates were considered, terminate + * the search, because it is doubtful that a much better + * candidate will be found */ + if (ncand == piv_lim) + goto done; + } + else if (suhl) + { /* j-th column has no eligible elements that satisfy to + * the stability criterion; Uwe Suhl suggests to exclude + * such column from further considerations until it + * becomes a column singleton; in hard cases this may + * significantly reduce the time needed to choose the + * pivot element */ + sgf_deactivate_col(j); + cs_prev[j] = cs_next[j] = j; + } + } + /* consider active rows containing len non-zeros */ + for (i = rs_head[len]; i != 0; i = rs_next[i]) + { /* big := max|v[i,*]| */ + if ((big = vr_max[i]) < 0.0) + { /* largest magnitude is unknown; compute it */ + for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + i_ptr < i_end; i_ptr++) + { if ((temp = sv_val[i_ptr]) < 0.0) + temp = -temp; + if (big < temp) + big = temp; + } + xassert(big > 0.0); + vr_max[i] = big; + } + /* find an element in i-th row, which is placed in the + * column with minimal number of non-zeros and satisfies to + * the stability condition (such element always exists) */ + min_i = min_j = 0, min_len = INT_MAX; + for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + i_ptr < i_end; i_ptr++) + { /* get column index of v[i,j] */ + j = sv_ind[i_ptr]; + /* if j-th column is not shorter, skip v[i,j] */ + if (vc_len[j] >= min_len) + continue; + /* if |v[i,j]| < piv_tol * max|v[i,*]|, skip v[i,j] */ + if ((temp = sv_val[i_ptr]) < 0.0) + temp = -temp; + if (temp < piv_tol * big) + continue; + /* v[i,j] is a better candidate */ + min_i = i, min_j = j, min_len = vc_len[j]; + /* if Markowitz cost of v[i,j] is not greater than + * (len-1)**2, v[i,j] can be chosen as the pivot right + * now; this heuristic reduces the search and works well + * in many cases */ + if (min_len <= len) + { p = min_i, q = min_j; + goto done; + } + } + /* i-th row has been scanned */ + if (min_i != 0) + { /* element v[min_i,min_j] is a next pivot candidate */ + ncand++; + /* compute its Markowitz cost */ + cost = (double)(len - 1) * (double)(min_len - 1); + /* if this element is better, choose it as the pivot */ + if (cost < best) + p = min_i, q = min_j, best = cost; + /* if piv_lim candidates were considered, terminate + * the search, because it is doubtful that a much better + * candidate will be found */ + if (ncand == piv_lim) + goto done; + } + else + { /* this can never be */ + xassert(min_i != min_i); + } + } + } +done: /* report the pivot to the factorization routine */ + *p_ = p, *q_ = q; + return (p == 0); +} + +/*********************************************************************** +* sgf_eliminate - perform gaussian elimination +* +* This routine performs elementary gaussian transformations in order +* to eliminate subdiagonal elements in k-th column of matrix +* U = P'* V * Q' using pivot element u[k,k], where k is the number of +* current elimination step, 1 <= k <= n. +* +* The parameters p and q specify, resp., row and column indices of the +* pivot element v[p,q] = u[k,k]. +* +* On entry the routine assumes that partially triangularized matrices +* L = P'* F * P and U = P'* V * Q' have the following structure: +* +* 1 k n 1 k n +* 1 1 . . . . . . . . . 1 x x x x x x x x x x +* x 1 . . . . . . . . . x x x x x x x x x +* x x 1 . . . . . . . . . x x x x x x x x +* x x x 1 . . . . . . . . . x x x x x x x +* k x x x x 1 . . . . . k . . . . * * * * * * +* x x x x _ 1 . . . . . . . . # * * * * * +* x x x x _ . 1 . . . . . . . # * * * * * +* x x x x _ . . 1 . . . . . . # * * * * * +* x x x x _ . . . 1 . . . . . # * * * * * +* n x x x x _ . . . . 1 n . . . . # * * * * * +* +* matrix L matrix U +* +* where rows and columns k, k+1, ..., n of matrix U constitute the +* active submatrix. Elements to be eliminated are marked by '#', and +* other elements of the active submatrix are marked by '*'. May note +* that each eliminated non-zero element u[i,k] of matrix U gives +* corresponding non-zero element l[i,k] of matrix L (marked by '_'). +* +* Actually all operations are performed on matrix V. It is assumed +* that the row-wise representation corresponds to matrix V, but the +* column-wise representation corresponds to the active submatrix of +* matrix V (or, more precisely, to its pattern, because only row +* indices for columns of the active submatrix are used on this stage). +* +* Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal +* elements u[i',k] = v[i,q], i'= k+1, k+2, ..., n, the routine applies +* the following elementary gaussian transformations: +* +* (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V), +* +* where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier stored to +* p-th column of matrix F to keep the main equality A = F * V +* (corresponding elements l[i',k] of matrix L are marked by '_' on the +* figure above). +* +* NOTE: On entry to the routine the working arrays flag and work +* should contain zeros. This status is retained by the routine +* on exit. */ + +int sgf_eliminate(SGF *sgf, int p, int q) +{ LUF *luf = sgf->luf; + int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fc_ref = luf->fc_ref; + int *fc_ptr = &sva->ptr[fc_ref-1]; + int *fc_len = &sva->len[fc_ref-1]; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int *vr_cap = &sva->cap[vr_ref-1]; + double *vr_piv = luf->vr_piv; + int vc_ref = luf->vc_ref; + int *vc_ptr = &sva->ptr[vc_ref-1]; + int *vc_len = &sva->len[vc_ref-1]; + int *vc_cap = &sva->cap[vc_ref-1]; + int *rs_head = sgf->rs_head; + int *rs_prev = sgf->rs_prev; + int *rs_next = sgf->rs_next; + int *cs_head = sgf->cs_head; + int *cs_prev = sgf->cs_prev; + int *cs_next = sgf->cs_next; + double *vr_max = sgf->vr_max; + char *flag = sgf->flag; + double *work = sgf->work; + double eps_tol = sgf->eps_tol; + int nnz_diff = 0; + int fill, i, i_ptr, i_end, j, j_ptr, j_end, ptr, len, loc, loc1; + double vpq, fip, vij; + xassert(1 <= p && p <= n); + xassert(1 <= q && q <= n); + /* remove p-th row from the active set; this row will never + * return there */ + sgf_deactivate_row(p); + /* process p-th (pivot) row */ + ptr = 0; + for (i_end = (i_ptr = vr_ptr[p]) + vr_len[p]; + i_ptr < i_end; i_ptr++) + { /* get column index of v[p,j] */ + j = sv_ind[i_ptr]; + if (j == q) + { /* save pointer to pivot v[p,q] */ + ptr = i_ptr; + } + else + { /* store v[p,j], j != q, to working array */ + flag[j] = 1; + work[j] = sv_val[i_ptr]; + } + /* remove j-th column from the active set; q-th column will + * never return there while other columns will return to the + * active set with new length */ + if (cs_next[j] == j) + { /* j-th column was marked by the pivoting routine according + * to Uwe Suhl's suggestion and is already inactive */ + xassert(cs_prev[j] == j); + } + else + sgf_deactivate_col(j); + nnz_diff -= vc_len[j]; + /* find and remove v[p,j] from j-th column */ + for (j_end = (j_ptr = vc_ptr[j]) + vc_len[j]; + sv_ind[j_ptr] != p; j_ptr++) + /* nop */; + xassert(j_ptr < j_end); + sv_ind[j_ptr] = sv_ind[j_end-1]; + vc_len[j]--; + } + /* save pivot v[p,q] and remove it from p-th row */ + xassert(ptr > 0); + vpq = vr_piv[p] = sv_val[ptr]; + sv_ind[ptr] = sv_ind[i_end-1]; + sv_val[ptr] = sv_val[i_end-1]; + vr_len[p]--; + /* if it is not planned to update matrix V, relocate p-th row to + * the right (static) part of SVA */ + if (!sgf->updat) + { len = vr_len[p]; + if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_make_static(sva, vr_ref-1+p); + } + /* copy the pattern (row indices) of q-th column of the active + * submatrix (from which v[p,q] has been just removed) to p-th + * column of matrix F (without unity diagonal element) */ + len = vc_len[q]; + if (len > 0) + { if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_reserve_cap(sva, fc_ref-1+p, len); + memcpy(&sv_ind[fc_ptr[p]], &sv_ind[vc_ptr[q]], + len * sizeof(int)); + fc_len[p] = len; + } + /* make q-th column of the active submatrix empty */ + vc_len[q] = 0; + /* transform non-pivot rows of the active submatrix */ + for (loc = fc_len[p]-1; loc >= 0; loc--) + { /* get row index of v[i,q] = row index of f[i,p] */ + i = sv_ind[fc_ptr[p] + loc]; + xassert(i != p); /* v[p,q] was removed */ + /* remove i-th row from the active set; this row will return + * there with new length */ + sgf_deactivate_row(i); + /* find v[i,q] in i-th row */ + for (i_end = (i_ptr = vr_ptr[i]) + vr_len[i]; + sv_ind[i_ptr] != q; i_ptr++) + /* nop */; + xassert(i_ptr < i_end); + /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] */ + fip = sv_val[fc_ptr[p] + loc] = sv_val[i_ptr] / vpq; + /* remove v[i,q] from i-th row */ + sv_ind[i_ptr] = sv_ind[i_end-1]; + sv_val[i_ptr] = sv_val[i_end-1]; + vr_len[i]--; + /* perform elementary gaussian transformation: + * (i-th row) := (i-th row) - f[i,p] * (p-th row) + * note that p-th row of V, which is in the working array, + * doesn't contain pivot v[p,q], and i-th row of V doesn't + * contain v[i,q] to be eliminated */ + /* walk thru i-th row and transform existing elements */ + fill = vr_len[p]; + for (i_end = (i_ptr = ptr = vr_ptr[i]) + vr_len[i]; + i_ptr < i_end; i_ptr++) + { /* get column index and value of v[i,j] */ + j = sv_ind[i_ptr]; + vij = sv_val[i_ptr]; + if (flag[j]) + { /* v[p,j] != 0 */ + flag[j] = 0, fill--; + /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */ + vij -= fip * work[j]; + if (-eps_tol < vij && vij < +eps_tol) + { /* new v[i,j] is close to zero; remove it from the + * active submatrix, i.e. replace it by exact zero */ + /* find and remove v[i,j] from j-th column */ + for (j_end = (j_ptr = vc_ptr[j]) + vc_len[j]; + sv_ind[j_ptr] != i; j_ptr++) + /* nop */; + xassert(j_ptr < j_end); + sv_ind[j_ptr] = sv_ind[j_end-1]; + vc_len[j]--; + continue; + } + } + /* keep new v[i,j] in i-th row */ + sv_ind[ptr] = j; + sv_val[ptr] = vij; + ptr++; + } + /* (new length of i-th row may decrease because of numerical + * cancellation) */ + vr_len[i] = len = ptr - vr_ptr[i]; + /* now flag[*] is the pattern of the set v[p,*] \ v[i,*], and + * fill is the number of non-zeros in this set */ + if (fill == 0) + { /* no fill-in occurs */ + /* walk thru p-th row and restore the column flags */ + for (i_end = (i_ptr = vr_ptr[p]) + vr_len[p]; + i_ptr < i_end; i_ptr++) + flag[sv_ind[i_ptr]] = 1; /* v[p,j] != 0 */ + goto skip; + } + /* up to fill new non-zero elements may appear in i-th row due + * to fill-in; reserve locations for these elements (note that + * actual length of i-th row is currently stored in len) */ + if (vr_cap[i] < len + fill) + { if (sva->r_ptr - sva->m_ptr < len + fill) + { sva_more_space(sva, len + fill); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, vr_ref-1+i, len + fill, 0); + } + vr_len[i] += fill; + /* walk thru p-th row and add new elements to i-th row */ + for (loc1 = vr_len[p]-1; loc1 >= 0; loc1--) + { /* get column index of v[p,j] */ + j = sv_ind[vr_ptr[p] + loc1]; + if (!flag[j]) + { /* restore j-th column flag */ + flag[j] = 1; + /* v[i,j] was computed earlier on transforming existing + * elements of i-th row */ + continue; + } + /* v[i,j] := 0 - f[i,p] * v[p,j] */ + vij = - fip * work[j]; + if (-eps_tol < vij && vij < +eps_tol) + { /* new v[i,j] is close to zero; do not add it to the + * active submatrix, i.e. replace it by exact zero */ + continue; + } + /* add new v[i,j] to i-th row */ + sv_ind[ptr = vr_ptr[i] + (len++)] = j; + sv_val[ptr] = vij; + /* add new v[i,j] to j-th column */ + if (vc_cap[j] == vc_len[j]) + { /* we reserve extra locations in j-th column to reduce + * further relocations of that column */ +#if 1 /* FIXME */ + /* use control parameter to specify the number of extra + * locations reserved */ + int need = vc_len[j] + 10; +#endif + if (sva->r_ptr - sva->m_ptr < need) + { sva_more_space(sva, need); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_enlarge_cap(sva, vc_ref-1+j, need, 1); + } + sv_ind[vc_ptr[j] + (vc_len[j]++)] = i; + } + /* set final length of i-th row just transformed */ + xassert(len <= vr_len[i]); + vr_len[i] = len; +skip: /* return i-th row to the active set with new length */ + sgf_activate_row(i); + /* since i-th row has been changed, largest magnitude of its + * elements becomes unknown */ + vr_max[i] = -1.0; + } + /* walk thru p-th (pivot) row */ + for (i_end = (i_ptr = vr_ptr[p]) + vr_len[p]; + i_ptr < i_end; i_ptr++) + { /* get column index of v[p,j] */ + j = sv_ind[i_ptr]; + xassert(j != q); /* v[p,q] was removed */ + /* return j-th column to the active set with new length */ + if (cs_next[j] == j && vc_len[j] != 1) + { /* j-th column was marked by the pivoting routine and it is + * still not a column singleton, so leave it incative */ + xassert(cs_prev[j] == j); + } + else + sgf_activate_col(j); + nnz_diff += vc_len[j]; + /* restore zero content of the working arrays */ + flag[j] = 0; + work[j] = 0.0; + } + /* return the difference between the numbers of non-zeros in the + * active submatrix on entry and on exit, resp. */ + return nnz_diff; +} + +/*********************************************************************** +* sgf_dense_lu - compute dense LU-factorization with full pivoting +* +* This routine performs Gaussian elimination with full pivoting to +* compute dense LU-factorization of the specified matrix A of order n +* in the form: +* +* A = P * L * U * Q, (1) +* +* where L is lower triangular matrix with unit diagonal, U is upper +* triangular matrix, P and Q are permutation matrices. +* +* On entry to the routine elements of matrix A = (a[i,j]) should be +* placed in the array elements a[0], ..., a[n^2-1] in dense row-wise +* format. On exit from the routine matrix A is replaced by factors L +* and U as follows: +* +* u[1,1] u[1,2] ... u[1,n-1] u[1,n] +* l[2,1] u[2,2] ... u[2,n-1] u[2,n] +* . . . . . . . . . . . . . . +* l[n-1,1] l[n-1,2] u[n-1,n-1] u[n-1,n] +* l[n,1] l[n,2] ... l[n,n-1] u[n,n] +* +* The unit diagonal elements of L are not stored. +* +* Information on permutations of rows and columns of active submatrix +* during factorization is accumulated by the routine as follows. Every +* time the routine permutes rows i and i' or columns j and j', it also +* permutes elements r[i-1] and r[i'-1] or c[j-1] and c[j'-1], resp. +* Thus, on entry to the routine elements r[0], r[1], ..., r[n-1] and +* c[0], c[1], ..., c[n-1] should be initialized by some integers that +* identify rows and columns of the original matrix A. +* +* If the factorization has been successfully computed, the routine +* returns zero. Otherwise, if on k-th elimination step, 1 <= k <= n, +* all elements of the active submatrix are close to zero, the routine +* returns k, in which case a partial factorization is stored in the +* array a. */ + +int sgf_dense_lu(int n, double a_[], int r[], int c[], double eps) +{ /* non-optimized version */ + int i, j, k, p, q, ref; + double akk, big, temp; +# define a(i,j) a_[(i)*n+(j)] + /* initially U = A, L = P = Q = I */ + /* main elimination loop */ + for (k = 0; k < n; k++) + { /* choose pivot u[p,q], k <= p, q <= n */ + p = q = -1, big = eps; + for (i = k; i < n; i++) + { for (j = k; j < n; j++) + { /* temp = |u[i,j]| */ + if ((temp = a(i,j)) < 0.0) + temp = -temp; + if (big < temp) + p = i, q = j, big = temp; + } + } + if (p < 0) + { /* k-th elimination step failed */ + return k+1; + } + /* permute rows k and p */ + if (k != p) + { for (j = 0; j < n; j++) + temp = a(k,j), a(k,j) = a(p,j), a(p,j) = temp; + ref = r[k], r[k] = r[p], r[p] = ref; + } + /* permute columns k and q */ + if (k != q) + { for (i = 0; i < n; i++) + temp = a(i,k), a(i,k) = a(i,q), a(i,q) = temp; + ref = c[k], c[k] = c[q], c[q] = ref; + } + /* now pivot is in position u[k,k] */ + akk = a(k,k); + /* eliminate subdiagonal elements u[k+1,k], ..., u[n,k] */ + for (i = k+1; i < n; i++) + { if (a(i,k) != 0.0) + { /* gaussian multiplier l[i,k] := u[i,k] / u[k,k] */ + temp = (a(i,k) /= akk); + /* (i-th row) := (i-th row) - l[i,k] * (k-th row) */ + for (j = k+1; j < n; j++) + a(i,j) -= temp * a(k,j); + } + } + } +# undef a + return 0; +} + +/*********************************************************************** +* sgf_dense_phase - compute LU-factorization (dense phase) +* +* This routine performs dense phase of computing LU-factorization. +* +* The aim is two-fold. First, the main factorization routine switches +* to dense phase when the active submatrix is relatively dense, so +* using dense format allows significantly reduces overheads needed to +* maintain sparse data structures. And second, that is more important, +* on dense phase full pivoting is used (rather than partial pivoting) +* that allows improving numerical stability, since round-off errors +* tend to increase on last steps of the elimination process. +* +* On entry the routine assumes that elimination steps 1, 2, ..., k-1 +* have been performed, so partially transformed matrices L = P'* F * P +* and U = P'* V * Q' have the following structure: +* +* 1 k n 1 k n +* 1 1 . . . . . . . . . 1 x x x x x x x x x x +* x 1 . . . . . . . . . x x x x x x x x x +* x x 1 . . . . . . . . . x x x x x x x x +* x x x 1 . . . . . . . . . x x x x x x x +* k x x x x 1 . . . . . k . . . . * * * * * * +* x x x x . 1 . . . . . . . . * * * * * * +* x x x x . . 1 . . . . . . . * * * * * * +* x x x x . . . 1 . . . . . . * * * * * * +* x x x x . . . . 1 . . . . . * * * * * * +* n x x x x . . . . . 1 n . . . . * * * * * * +* +* matrix L matrix U +* +* where rows and columns k, k+1, ..., n of matrix U constitute the +* active submatrix A~, whose elements are marked by '*'. +* +* The routine copies the active submatrix A~ to a working array in +* dense format, compute dense factorization A~ = P~* L~* U~* Q~ using +* full pivoting, and then copies non-zero elements of factors L~ and +* U~ back to factors L and U (more precisely, to factors F and V). +* +* If the factorization has been successfully computed, the routine +* returns zero. Otherwise, if on k-th elimination step, 1 <= k <= n, +* all elements of the active submatrix are close to zero, the routine +* returns k (information on linearly dependent rows/columns in this +* case is provided by matrices P and Q). */ + +int sgf_dense_phase(LUF *luf, int k, int updat) +{ int n = luf->n; + SVA *sva = luf->sva; + int *sv_ind = sva->ind; + double *sv_val = sva->val; + int fc_ref = luf->fc_ref; + int *fc_ptr = &sva->ptr[fc_ref-1]; + int *fc_len = &sva->len[fc_ref-1]; + int *fc_cap = &sva->cap[fc_ref-1]; + int vr_ref = luf->vr_ref; + int *vr_ptr = &sva->ptr[vr_ref-1]; + int *vr_len = &sva->len[vr_ref-1]; + int *vr_cap = &sva->cap[vr_ref-1]; + double *vr_piv = luf->vr_piv; + int vc_ref = luf->vc_ref; + int *vc_len = &sva->len[vc_ref-1]; + int *pp_inv = luf->pp_inv; + int *pp_ind = luf->pp_ind; + int *qq_ind = luf->qq_ind; + int *qq_inv = luf->qq_inv; + int a_end, a_ptr, end, i, ia, ii, j, ja, jj, ka, len, na, ne, + need, ptr; + double *a_; + xassert(1 <= k && k <= n); + /* active columns of V are not longer needed; make them empty */ + for (jj = k; jj <= n; jj++) + { /* jj is number of active column of U = P'* V * Q' */ + vc_len[qq_ind[jj]] = 0; + } + /* determine order of active submatrix A~ of matrix U */ + na = n - k + 1; + xassert(1 <= na && na <= n); + /* determine number of elements in dense triangular factor (L~ or + * U~), except diagonal elements */ + ne = na * (na - 1) / 2; + /* we allocate active submatrix A~ in free (middle) part of SVA; + * to avoid defragmentation that could destroy A~ we also should + * reserve ne locations to build rows of V from rows of U~ and ne + * locations to build columns of F from columns of L~ */ + need = na * na + ne + ne; + if (sva->r_ptr - sva->m_ptr < need) + { sva_more_space(sva, need); + sv_ind = sva->ind; + sv_val = sva->val; + } + /* free (middle) part of SVA is structured as follows: + * end of left (dynamic) part + * ne free locations for new rows of V + * na free locations for active submatrix A~ + * unused locations, if any + * ne free locations for new columns of F + * beginning of right (static) part */ + a_ptr = sva->m_ptr + ne; + a_end = a_ptr + na * na; + /* copy active submatrix A~ from matrix V to working array in + * dense row-wise format */ + a_ = &sva->val[a_ptr]; +# define a(ia, ja) a_[((ia) - 1) * na + ((ja) - 1)] + for (ia = 1; ia <= na; ia++) + { /* clear ia-th row of A~ */ + for (ja = 1; ja <= na; ja++) + a(ia, ja) = 0.0; + /* ia-th row of A~ = (k-1+ia)-th row of U = i-th row of V */ + i = pp_inv[k-1+ia]; + ptr = vr_ptr[i]; + end = ptr + vr_len[i]; + for (; ptr < end; ptr++) + a(ia, qq_inv[sv_ind[ptr]]-k+1) = sv_val[ptr]; + /* i-th row of V is no longer needed; make it empty */ + vr_len[i] = 0; + } + /* compute dense factorization A~ = P~* L~* U~* Q~ */ +#if 1 /* FIXME: epsilon tolerance */ + ka = sgf_dense_lu(na, &a(1, 1), &pp_inv[k], &qq_ind[k], 1e-20); +#endif + /* rows of U with numbers pp_inv[k, k+1, ..., n] were permuted + * due to row permutations of A~; update matrix P using P~ */ + for (ii = k; ii <= n; ii++) + pp_ind[pp_inv[ii]] = ii; + /* columns of U with numbers qq_ind[k, k+1, ..., n] were permuted + * due to column permutations of A~; update matrix Q using Q~ */ + for (jj = k; jj <= n; jj++) + qq_inv[qq_ind[jj]] = jj; + /* check if dense factorization is complete */ + if (ka != 0) + { /* A~ is singular to working precision */ + /* information on linearly dependent rows/columns is provided + * by matrices P and Q */ + xassert(1 <= ka && ka <= na); + return k - 1 + ka; + } + /* build new rows of V from rows of U~ */ + for (ia = 1; ia <= na; ia++) + { /* ia-th row of U~ = (k-1+ia)-th row of U = i-th row of V */ + i = pp_inv[k-1+ia]; + xassert(vr_len[i] == 0); + /* store diagonal element u~[ia,ia] */ + vr_piv[i] = a(ia, ia); + /* determine number of non-zero non-diagonal elements in ia-th + * row of U~ */ + len = 0; + for (ja = ia+1; ja <= na; ja++) + { if (a(ia, ja) != 0.0) + len++; + } + /* reserve len locations for i-th row of matrix V in left + * (dynamic) part of SVA */ + if (vr_cap[i] < len) + { /* there should be enough room in free part of SVA */ + xassert(sva->r_ptr - sva->m_ptr >= len); + sva_enlarge_cap(sva, vr_ref-1+i, len, 0); + /* left part of SVA should not overlap matrix A~ */ + xassert(sva->m_ptr <= a_ptr); + } + /* copy non-zero non-diaginal elements of ia-th row of U~ to + * i-th row of V */ + ptr = vr_ptr[i]; + for (ja = ia+1; ja <= na; ja++) + { if (a(ia, ja) != 0.0) + { sv_ind[ptr] = qq_ind[k-1+ja]; + sv_val[ptr] = a(ia, ja); + ptr++; + } + } + xassert(ptr - vr_ptr[i] == len); + vr_len[i] = len; + } + /* build new columns of F from columns of L~ */ + for (ja = 1; ja <= na; ja++) + { /* ja-th column of L~ = (k-1+ja)-th column of L = j-th column + * of F */ + j = pp_inv[k-1+ja]; + xassert(fc_len[j] == 0); + xassert(fc_cap[j] == 0); + /* determine number of non-zero non-diagonal elements in ja-th + * column of L~ */ + len = 0; + for (ia = ja+1; ia <= na; ia++) + { if (a(ia, ja) != 0.0) + len++; + } + /* reserve len locations for j-th column of matrix F in right + * (static) part of SVA */ + /* there should be enough room in free part of SVA */ + xassert(sva->r_ptr - sva->m_ptr >= len); + if (len > 0) + sva_reserve_cap(sva, fc_ref-1+j, len); + /* right part of SVA should not overlap matrix A~ */ + xassert(a_end <= sva->r_ptr); + /* copy non-zero non-diagonal elements of ja-th column of L~ + * to j-th column of F */ + ptr = fc_ptr[j]; + for (ia = ja+1; ia <= na; ia++) + { if (a(ia, ja) != 0.0) + { sv_ind[ptr] = pp_inv[k-1+ia]; + sv_val[ptr] = a(ia, ja); + ptr++; + } + } + xassert(ptr - fc_ptr[j] == len); + fc_len[j] = len; + } + /* factors L~ and U~ are no longer needed */ +# undef a + /* if it is not planned to update matrix V, relocate all its new + * rows to the right (static) part of SVA */ + if (!updat) + { for (ia = 1; ia <= na; ia++) + { i = pp_inv[k-1+ia]; + len = vr_len[i]; + if (sva->r_ptr - sva->m_ptr < len) + { sva_more_space(sva, len); + sv_ind = sva->ind; + sv_val = sva->val; + } + sva_make_static(sva, vr_ref-1+i); + } + } + return 0; +} + +/*********************************************************************** +* sgf_factorize - compute LU-factorization (main routine) +* +* This routine computes sparse LU-factorization of specified matrix A +* using Gaussian elimination. +* +* On entry to the routine matrix V = A should be stored in column-wise +* format. +* +* If the factorization has been successfully computed, the routine +* returns zero. Otherwise, if on k-th elimination step, 1 <= k <= n, +* all elements of the active submatrix are close to zero, the routine +* returns k (information on linearly dependent rows/columns in this +* case is provided by matrices P and Q). */ + +#if 1 /* 21/II-2016 */ +/* If the matrix A is structurally singular, the routine returns -1. +* NOTE: This case can be detected only if the singl flag is set. */ +#endif + +int sgf_factorize(SGF *sgf, int singl) +{ LUF *luf = sgf->luf; + int n = luf->n; + SVA *sva = luf->sva; + int vr_ref = luf->vr_ref; + int *vr_len = &sva->len[vr_ref-1]; + double *vr_piv = luf->vr_piv; + int vc_ref = luf->vc_ref; + int *vc_len = &sva->len[vc_ref-1]; + int *pp_ind = luf->pp_ind; + int *pp_inv = luf->pp_inv; + int *qq_ind = luf->qq_ind; + int *qq_inv = luf->qq_inv; + int *rs_head = sgf->rs_head; + int *rs_prev = sgf->rs_prev; + int *rs_next = sgf->rs_next; + int *cs_head = sgf->cs_head; + int *cs_prev = sgf->cs_prev; + int *cs_next = sgf->cs_next; + double *vr_max = sgf->vr_max; + char *flag = sgf->flag; + double *work = sgf->work; + int i, j, k, k1, k2, p, q, nnz; + /* build matrix V = A in row-wise format */ + luf_build_v_rows(luf, rs_prev); + /* P := Q := I, so V = U = A, F = L = I */ + for (k = 1; k <= n; k++) + { vr_piv[k] = 0.0; + pp_ind[k] = pp_inv[k] = qq_ind[k] = qq_inv[k] = k; + } +#ifdef GLP_DEBUG + sva_check_area(sva); + luf_check_all(luf, 1); +#endif + /* perform singleton phase, if required */ + if (!singl) + { /* assume that nucleus is entire matrix U */ + k2 = 1; + } + else + { /* minimize nucleus size */ +#if 0 /* 21/II-2016 */ + sgf_reduce_nuc(luf, &k1, &k2, rs_prev, rs_next); +#else + if (sgf_reduce_nuc(luf, &k1, &k2, rs_prev, rs_next)) + return -1; +#endif +#ifdef GLP_DEBUG + xprintf("n = %d; k1 = %d; k2 = %d\n", n, k1, k2); +#endif + /* perform singleton phase */ + k2 = sgf_singl_phase(luf, k1, k2, sgf->updat, rs_prev, work); + } +#ifdef GLP_DEBUG + sva_check_area(sva); + luf_check_all(luf, k2); +#endif + /* initialize working arrays */ + rs_head[0] = cs_head[0] = 0; + for (k = 1; k <= n; k++) + { rs_head[k] = cs_head[k] = 0; + vr_max[k] = -1.0; + flag[k] = 0; + work[k] = 0.0; + } + /* build lists of active rows and columns of matrix V; determine + * number of non-zeros in initial active submatrix */ + nnz = 0; + for (k = k2; k <= n; k++) + { i = pp_inv[k]; + sgf_activate_row(i); + nnz += vr_len[i]; + j = qq_ind[k]; + sgf_activate_col(j); + } + /* main factorization loop */ + for (k = k2; k <= n; k++) + { int na; + double den; + /* calculate density of active submatrix */ + na = n - k + 1; /* order of active submatrix */ +#if 0 /* 21/VIII-2014 */ + den = (double)nnz / (double)(na * na); +#else + den = (double)nnz / ((double)(na) * (double)(na)); +#endif + /* if active submatrix is relatively dense, switch to dense + * phase */ +#if 1 /* FIXME */ + if (na >= 5 && den >= 0.71) + { +#ifdef GLP_DEBUG + xprintf("na = %d; nnz = %d; den = %g\n", na, nnz, den); +#endif + break; + } +#endif + /* choose pivot v[p,q] */ + if (sgf_choose_pivot(sgf, &p, &q) != 0) + return k; /* failure */ + /* u[i,j] = v[p,q], k <= i, j <= n */ + i = pp_ind[p]; + xassert(k <= i && i <= n); + j = qq_inv[q]; + xassert(k <= j && j <= n); + /* move u[i,j] to position u[k,k] by implicit permutations of + * rows and columns of matrix U */ + luf_swap_u_rows(k, i); + luf_swap_u_cols(k, j); + /* perform gaussian elimination */ + nnz += sgf_eliminate(sgf, p, q); + } +#if 1 /* FIXME */ + if (k <= n) + { /* continue computing factorization in dense mode */ +#ifdef GLP_DEBUG + sva_check_area(sva); + luf_check_all(luf, k); +#endif + k = sgf_dense_phase(luf, k, sgf->updat); + if (k != 0) + return k; /* failure */ + } +#endif +#ifdef GLP_DEBUG + sva_check_area(sva); + luf_check_all(luf, n+1); +#endif + /* defragment SVA; currently all columns of V are empty, so they + * will have zero capacity as required by luf_build_v_cols */ + sva_defrag_area(sva); + /* build matrix F in row-wise format */ + luf_build_f_rows(luf, rs_head); + /* build matrix V in column-wise format */ + luf_build_v_cols(luf, sgf->updat, rs_head); + return 0; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/sgf.h b/glpk-5.0/src/bflib/sgf.h new file mode 100644 index 0000000..4bcac21 --- /dev/null +++ b/glpk-5.0/src/bflib/sgf.h @@ -0,0 +1,201 @@ +/* sgf.h (sparse Gaussian factorizer) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef SGF_H +#define SGF_H + +#include "luf.h" + +typedef struct SGF SGF; + +struct SGF +{ /* sparse Gaussian factorizer workspace */ + LUF *luf; + /* LU-factorization being computed */ + /*--------------------------------------------------------------*/ + /* to efficiently choose pivot elements according to Markowitz + * strategy, the search technique proposed by Iain Duff is used; + * it is based on using two families of sets {R[0], ..., R[n]} + * and {C[0], ..., C[n]}, where R[k] and C[k], 0 <= k <= n, are, + * respectively, sets of rows and columns of the active submatrix + * of matrix V having k non-zeros (i.e. whose length is k); each + * set R[k] and C[k] is implemented as a doubly linked list */ + int *rs_head; /* int rs_head[1+n]; */ + /* rs_head[k], 0 <= k <= n, is the number of first row, which + * has k non-zeros in the active submatrix */ + int *rs_prev; /* int rs_prev[1+n]; */ + /* rs_prev[0] is not used; + * rs_prev[i], 1 <= i <= n, is the number of previous row, which + * has the same number of non-zeros as i-th row; + * rs_prev[i] < 0 means that i-th row is inactive */ + int *rs_next; /* int rs_next[1+n]; */ + /* rs_next[0] is not used; + * rs_next[i], 1 <= i <= n, is the number of next row, which has + * the same number of non-zeros as i-th row; + * rs_next[i] < 0 means that i-th row is inactive */ + int *cs_head; /* int cs_head[1+n]; */ + /* cs_head[k], 0 <= k <= n, is the number of first column, which + * has k non-zeros in the active submatrix */ + int *cs_prev; /* int cs_prev[1+n]; */ + /* cs_prev[0] is not used; + * cs_prev[j], 1 <= j <= n, is the number of previous column, + * which has the same number of non-zeros as j-th column; + * cs_prev[j] < 0 means that j-th column is inactive */ + int *cs_next; /* int cs_next[1+n]; */ + /* cs_next[0] is not used; + * cs_next[j], 1 <= j <= n, is the number of next column, which + * has the same number of non-zeros as j-th column; + * cs_next[j] < 0 means that j-th column is inactive */ + /* NOTE: cs_prev[j] = cs_next[j] = j means that j-th column was + * temporarily removed from corresponding set C[k] by the + * pivoting routine according to Uwe Suhl's heuristic */ + /*--------------------------------------------------------------*/ + /* working arrays */ + double *vr_max; /* int vr_max[1+n]; */ + /* vr_max[0] is not used; + * vr_max[i], 1 <= i <= n, is used only if i-th row of matrix V + * is active (i.e. belongs to the active submatrix), and is the + * largest magnitude of elements in that row; if vr_max[i] < 0, + * the largest magnitude is unknown yet */ + char *flag; /* char flag[1+n]; */ + /* boolean working array */ + double *work; /* double work[1+n]; */ + /* floating-point working array */ + /*--------------------------------------------------------------*/ + /* control parameters */ + int updat; + /* if this flag is set, the matrix V is assumed to be updatable; + * in this case factorized (non-active) part of V is stored in + * the left part of SVA rather than in its right part */ + double piv_tol; + /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j] + * of the active submatrix fits to be pivot if it satisfies to + * the stability criterion |v[i,j]| >= piv_tol * max |v[i,*]|, + * i.e. if it is not very small in the magnitude among other + * elements in the same row; decreasing this parameter gives + * better sparsity at the expense of numerical accuracy and vice + * versa */ + int piv_lim; + /* maximal allowable number of pivot candidates to be considered; + * if piv_lim pivot candidates have been considered, the pivoting + * routine terminates the search with the best candidate found */ + int suhl; + /* if this flag is set, the pivoting routine applies a heuristic + * proposed by Uwe Suhl: if a column of the active submatrix has + * no eligible pivot candidates (i.e. all its elements do not + * satisfy to the stability criterion), the routine excludes it + * from futher consideration until it becomes column singleton; + * in many cases this allows reducing the time needed to choose + * the pivot */ + double eps_tol; + /* epsilon tolerance; each element of the active submatrix, whose + * magnitude is less than eps_tol, is replaced by exact zero */ +#if 0 /* FIXME */ + double den_lim; + /* density limit; if the density of the active submatrix reaches + * this limit, the factorization routine switches from sparse to + * dense mode */ +#endif +}; + +#define sgf_activate_row(i) \ + do \ + { int len = vr_len[i]; \ + rs_prev[i] = 0; \ + rs_next[i] = rs_head[len]; \ + if (rs_next[i] != 0) \ + rs_prev[rs_next[i]] = i; \ + rs_head[len] = i; \ + } while (0) +/* include i-th row of matrix V in active set R[len] */ + +#define sgf_deactivate_row(i) \ + do \ + { if (rs_prev[i] == 0) \ + rs_head[vr_len[i]] = rs_next[i]; \ + else \ + rs_next[rs_prev[i]] = rs_next[i]; \ + if (rs_next[i] == 0) \ + ; \ + else \ + rs_prev[rs_next[i]] = rs_prev[i]; \ + rs_prev[i] = rs_next[i] = -1; \ + } while (0) +/* remove i-th row of matrix V from active set R[len] */ + +#define sgf_activate_col(j) \ + do \ + { int len = vc_len[j]; \ + cs_prev[j] = 0; \ + cs_next[j] = cs_head[len]; \ + if (cs_next[j] != 0) \ + cs_prev[cs_next[j]] = j; \ + cs_head[len] = j; \ + } while (0) +/* include j-th column of matrix V in active set C[len] */ + +#define sgf_deactivate_col(j) \ + do \ + { if (cs_prev[j] == 0) \ + cs_head[vc_len[j]] = cs_next[j]; \ + else \ + cs_next[cs_prev[j]] = cs_next[j]; \ + if (cs_next[j] == 0) \ + ; \ + else \ + cs_prev[cs_next[j]] = cs_prev[j]; \ + cs_prev[j] = cs_next[j] = -1; \ + } while (0) +/* remove j-th column of matrix V from active set C[len] */ + +#define sgf_reduce_nuc _glp_sgf_reduce_nuc +int sgf_reduce_nuc(LUF *luf, int *k1, int *k2, int cnt[/*1+n*/], + int list[/*1+n*/]); +/* initial reordering to minimize nucleus size */ + +#define sgf_singl_phase _glp_sgf_singl_phase +int sgf_singl_phase(LUF *luf, int k1, int k2, int updat, + int ind[/*1+n*/], double val[/*1+n*/]); +/* compute LU-factorization (singleton phase) */ + +#define sgf_choose_pivot _glp_sgf_choose_pivot +int sgf_choose_pivot(SGF *sgf, int *p, int *q); +/* choose pivot element v[p,q] */ + +#define sgf_eliminate _glp_sgf_eliminate +int sgf_eliminate(SGF *sgf, int p, int q); +/* perform gaussian elimination */ + +#define sgf_dense_lu _glp_sgf_dense_lu +int sgf_dense_lu(int n, double a[], int r[], int c[], double eps); +/* compute dense LU-factorization with full pivoting */ + +#define sgf_dense_phase _glp_sgf_dense_phase +int sgf_dense_phase(LUF *luf, int k, int updat); +/* compute LU-factorization (dense phase) */ + +#define sgf_factorize _glp_sgf_factorize +int sgf_factorize(SGF *sgf, int singl); +/* compute LU-factorization (main routine) */ + +#endif + +/* eof */ diff --git a/glpk-5.0/src/bflib/sva.c b/glpk-5.0/src/bflib/sva.c new file mode 100644 index 0000000..0065b78 --- /dev/null +++ b/glpk-5.0/src/bflib/sva.c @@ -0,0 +1,570 @@ +/* sva.c (sparse vector area) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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 "sva.h" + +/*********************************************************************** +* sva_create_area - create sparse vector area (SVA) +* +* This routine creates the sparse vector area (SVA), which initially +* is empty. +* +* The parameter n_max specifies the initial number of vectors that can +* be allocated in the SVA, n_max > 0. +* +* The parameter size specifies the initial number of free locations in +* the SVA, size > 0. +* +* On exit the routine returns a pointer to the SVA created. */ + +SVA *sva_create_area(int n_max, int size) +{ SVA *sva; + xassert(0 < n_max && n_max < INT_MAX); + xassert(0 < size && size < INT_MAX); + sva = talloc(1, SVA); + sva->n_max = n_max; + sva->n = 0; + sva->ptr = talloc(1+n_max, int); + sva->len = talloc(1+n_max, int); + sva->cap = talloc(1+n_max, int); + sva->size = size; + sva->m_ptr = 1; + sva->r_ptr = size+1; + sva->head = sva->tail = 0; + sva->prev = talloc(1+n_max, int); + sva->next = talloc(1+n_max, int); + sva->ind = talloc(1+size, int); + sva->val = talloc(1+size, double); + sva->talky = 0; + return sva; +} + +/*********************************************************************** +* sva_alloc_vecs - allocate new vectors in SVA +* +* This routine allocates nnn new empty vectors, nnn > 0, in the sparse +* vector area (SVA). +* +* The new vectors are assigned reference numbers k, k+1, ..., k+nnn-1, +* where k is a reference number assigned to the very first new vector, +* which is returned by the routine on exit. */ + +int sva_alloc_vecs(SVA *sva, int nnn) +{ int n = sva->n; + int n_max = sva->n_max; + int *ptr = sva->ptr; + int *len = sva->len; + int *cap = sva->cap; + int *prev = sva->prev; + int *next = sva->next; + int k, new_n; +#if 1 + if (sva->talky) + xprintf("sva_alloc_vecs: nnn = %d\n", nnn); +#endif + xassert(nnn > 0); + /* determine new number of vectors in SVA */ + new_n = n + nnn; + xassert(new_n > n); + if (n_max < new_n) + { /* enlarge the SVA arrays */ + while (n_max < new_n) + { n_max += n_max; + xassert(n_max > 0); + } + sva->n_max = n_max; + sva->ptr = ptr = trealloc(ptr, 1+n_max, int); + sva->len = len = trealloc(len, 1+n_max, int); + sva->cap = cap = trealloc(cap, 1+n_max, int); + sva->prev = prev = trealloc(prev, 1+n_max, int); + sva->next = next = trealloc(next, 1+n_max, int); + } + /* initialize new vectors */ + sva->n = new_n; + for (k = n+1; k <= new_n; k++) + { ptr[k] = len[k] = cap[k] = 0; + prev[k] = next[k] = -1; + } +#if 1 + if (sva->talky) + xprintf("now sva->n_max = %d, sva->n = %d\n", + sva->n_max, sva->n); +#endif + /* return reference number of very first new vector */ + return n+1; +} + +/*********************************************************************** +* sva_resize_area - change size of SVA storage +* +* This routine increases or decrases the size of the SVA storage by +* reallocating it. +* +* The parameter delta specifies the number of location by which the +* current size of the SVA storage should be increased (if delta > 0) +* or decreased (if delta < 0). Note that if delta is negative, it +* should not be less than the current size of the middle part. +* +* As a result of this operation the size of the middle part of SVA is +* increased/decreased by delta locations. +* +* NOTE: This operation changes ptr[k] for all vectors stored in the +* right part of SVA. */ + +void sva_resize_area(SVA *sva, int delta) +{ int n = sva->n; + int *ptr = sva->ptr; + int size = sva->size; + int m_ptr = sva->m_ptr; + int r_ptr = sva->r_ptr; + int k, r_size; +#if 1 + if (sva->talky) + xprintf("sva_resize_area: delta = %d\n", delta); +#endif + xassert(delta != 0); + /* determine size of the right part, in locations */ + r_size = size - r_ptr + 1; + /* relocate the right part in case of negative delta */ + if (delta < 0) + { xassert(delta >= m_ptr - r_ptr); + sva->r_ptr += delta; + memmove(&sva->ind[sva->r_ptr], &sva->ind[r_ptr], + r_size * sizeof(int)); + memmove(&sva->val[sva->r_ptr], &sva->val[r_ptr], + r_size * sizeof(double)); + } + /* reallocate the storage arrays */ + xassert(delta < INT_MAX - sva->size); + sva->size += delta; + sva->ind = trealloc(sva->ind, 1+sva->size, int); + sva->val = trealloc(sva->val, 1+sva->size, double); + /* relocate the right part in case of positive delta */ + if (delta > 0) + { sva->r_ptr += delta; + memmove(&sva->ind[sva->r_ptr], &sva->ind[r_ptr], + r_size * sizeof(int)); + memmove(&sva->val[sva->r_ptr], &sva->val[r_ptr], + r_size * sizeof(double)); + } + /* update pointers to vectors stored in the right part */ + for (k = 1; k <= n; k++) + { if (ptr[k] >= r_ptr) + ptr[k] += delta; + } +#if 1 + if (sva->talky) + xprintf("now sva->size = %d\n", sva->size); +#endif + return; +} + +/*********************************************************************** +* sva_defrag_area - defragment left part of SVA +* +* This routine performs "garbage" collection to defragment the left +* part of SVA. +* +* NOTE: This operation may change ptr[k] and cap[k] for all vectors +* stored in the left part of SVA. */ + +void sva_defrag_area(SVA *sva) +{ int *ptr = sva->ptr; + int *len = sva->len; + int *cap = sva->cap; + int *prev = sva->prev; + int *next = sva->next; + int *ind = sva->ind; + double *val = sva->val; + int k, next_k, ptr_k, len_k, m_ptr, head, tail; +#if 1 + if (sva->talky) + { xprintf("sva_defrag_area:\n"); + xprintf("before defragmenting = %d %d %d\n", sva->m_ptr - 1, + sva->r_ptr - sva->m_ptr, sva->size + 1 - sva->r_ptr); + } +#endif + m_ptr = 1; + head = tail = 0; + /* walk through the linked list of vectors stored in the left + * part of SVA */ + for (k = sva->head; k != 0; k = next_k) + { /* save number of next vector in the list */ + next_k = next[k]; + /* determine length of k-th vector */ + len_k = len[k]; + if (len_k == 0) + { /* k-th vector is empty; remove it from the left part */ + ptr[k] = cap[k] = 0; + prev[k] = next[k] = -1; + } + else + { /* determine pointer to first location of k-th vector */ + ptr_k = ptr[k]; + xassert(m_ptr <= ptr_k); + /* relocate k-th vector to the beginning of the left part, + * if necessary */ + if (m_ptr < ptr_k) + { memmove(&ind[m_ptr], &ind[ptr_k], + len_k * sizeof(int)); + memmove(&val[m_ptr], &val[ptr_k], + len_k * sizeof(double)); + ptr[k] = m_ptr; + } + /* remove unused locations from k-th vector */ + cap[k] = len_k; + /* the left part of SVA has been enlarged */ + m_ptr += len_k; + /* add k-th vector to the end of the new linked list */ + prev[k] = tail; + next[k] = 0; + if (head == 0) + head = k; + else + next[tail] = k; + tail = k; + } + } + /* set new pointer to the middle part of SVA */ + xassert(m_ptr <= sva->r_ptr); + sva->m_ptr = m_ptr; + /* set new head and tail of the linked list */ + sva->head = head; + sva->tail = tail; +#if 1 + if (sva->talky) + xprintf("after defragmenting = %d %d %d\n", sva->m_ptr - 1, + sva->r_ptr - sva->m_ptr, sva->size + 1 - sva->r_ptr); +#endif + return; +} + +/*********************************************************************** +* sva_more_space - increase size of middle (free) part of SVA +* +* This routine increases the size of the middle (free) part of the +* sparse vector area (SVA). +* +* The parameter m_size specifies the minimal size, in locations, of +* the middle part to be provided. This new size should be greater than +* the current size of the middle part. +* +* First, the routine defragments the left part of SVA. Then, if the +* size of the left part has not sufficiently increased, the routine +* increases the total size of the SVA storage by reallocating it. */ + +void sva_more_space(SVA *sva, int m_size) +{ int size, delta; +#if 1 + if (sva->talky) + xprintf("sva_more_space: m_size = %d\n", m_size); +#endif + xassert(m_size > sva->r_ptr - sva->m_ptr); + /* defragment the left part */ + sva_defrag_area(sva); + /* set, heuristically, the minimal size of the middle part to be + * not less than the size of the defragmented left part */ + if (m_size < sva->m_ptr - 1) + m_size = sva->m_ptr - 1; + /* if there is still not enough room, increase the total size of + * the SVA storage */ + if (sva->r_ptr - sva->m_ptr < m_size) + { size = sva->size; /* new sva size */ + for (;;) + { delta = size - sva->size; + if (sva->r_ptr - sva->m_ptr + delta >= m_size) + break; + size += size; + xassert(size > 0); + } + sva_resize_area(sva, delta); + xassert(sva->r_ptr - sva->m_ptr >= m_size); + } + return; +} + +/*********************************************************************** +* sva_enlarge_cap - enlarge capacity of specified vector +* +* This routine enlarges the current capacity of the specified vector +* by relocating its content. +* +* The parameter k specifies the reference number of the vector whose +* capacity should be enlarged, 1 <= k <= n. This vector should either +* have zero capacity or be stored in the left (dynamic) part of SVA. +* +* The parameter new_cap specifies the new capacity of the vector, +* in locations. This new capacity should be greater than the current +* capacity of the vector. +* +* The parameter skip is a flag. If this flag is set, the routine does +* *not* copy numerical values of elements of the vector on relocating +* its content, i.e. only element indices are copied. +* +* NOTE: On entry to the routine the middle part of SVA should have at +* least new_cap free locations. */ + +void sva_enlarge_cap(SVA *sva, int k, int new_cap, int skip) +{ int *ptr = sva->ptr; + int *len = sva->len; + int *cap = sva->cap; + int *prev = sva->prev; + int *next = sva->next; + int *ind = sva->ind; + double *val = sva->val; + xassert(1 <= k && k <= sva->n); + xassert(new_cap > cap[k]); + /* there should be at least new_cap free locations */ + xassert(sva->r_ptr - sva->m_ptr >= new_cap); + /* relocate the vector */ + if (cap[k] == 0) + { /* the vector is empty */ + xassert(ptr[k] == 0); + xassert(len[k] == 0); + } + else + { /* the vector has non-zero capacity */ + xassert(ptr[k] + len[k] <= sva->m_ptr); + /* copy the current vector content to the beginning of the + * middle part */ + if (len[k] > 0) + { memcpy(&ind[sva->m_ptr], &ind[ptr[k]], + len[k] * sizeof(int)); + if (!skip) + memcpy(&val[sva->m_ptr], &val[ptr[k]], + len[k] * sizeof(double)); + } + /* remove the vector from the linked list */ + if (prev[k] == 0) + sva->head = next[k]; + else + { /* preceding vector exists; increase its capacity */ + cap[prev[k]] += cap[k]; + next[prev[k]] = next[k]; + } + if (next[k] == 0) + sva->tail = prev[k]; + else + prev[next[k]] = prev[k]; + } + /* set new pointer and capacity of the vector */ + ptr[k] = sva->m_ptr; + cap[k] = new_cap; + /* add the vector to the end of the linked list */ + prev[k] = sva->tail; + next[k] = 0; + if (sva->head == 0) + sva->head = k; + else + next[sva->tail] = k; + sva->tail = k; + /* new_cap free locations have been consumed */ + sva->m_ptr += new_cap; + xassert(sva->m_ptr <= sva->r_ptr); + return; +} + +/*********************************************************************** +* sva_reserve_cap - reserve locations for specified vector +* +* This routine reserves locations for the specified vector in the +* right (static) part of SVA. +* +* The parameter k specifies the reference number of the vector (this +* vector should have zero capacity), 1 <= k <= n. +* +* The parameter new_cap specifies a non-zero capacity of the vector, +* in locations. +* +* NOTE: On entry to the routine the middle part of SVA should have at +* least new_cap free locations. */ + +void sva_reserve_cap(SVA *sva, int k, int new_cap) +{ int *ptr = sva->ptr; + int *len = sva->len; + int *cap = sva->cap; + xassert(1 <= k && k <= sva->n); + xassert(new_cap > 0); + xassert(ptr[k] == 0 && len[k] == 0 && cap[k] == 0); + /* there should be at least new_cap free locations */ + xassert(sva->r_ptr - sva->m_ptr >= new_cap); + /* set the pointer and capacity of the vector */ + ptr[k] = sva->r_ptr - new_cap; + cap[k] = new_cap; + /* new_cap free locations have been consumed */ + sva->r_ptr -= new_cap; + return; +} + +/*********************************************************************** +* sva_make_static - relocate specified vector to right part of SVA +* +* Assuming that the specified vector is stored in the left (dynamic) +* part of SVA, this routine makes the vector static by relocating its +* content to the right (static) part of SVA. However, if the specified +* vector has zero capacity, the routine does nothing. +* +* The parameter k specifies the reference number of the vector to be +* relocated, 1 <= k <= n. +* +* NOTE: On entry to the routine the middle part of SVA should have at +* least len[k] free locations, where len[k] is the length of the +* vector to be relocated. */ + +void sva_make_static(SVA *sva, int k) +{ int *ptr = sva->ptr; + int *len = sva->len; + int *cap = sva->cap; + int *prev = sva->prev; + int *next = sva->next; + int *ind = sva->ind; + double *val = sva->val; + int ptr_k, len_k; + xassert(1 <= k && k <= sva->n); + /* if the vector has zero capacity, do nothing */ + if (cap[k] == 0) + { xassert(ptr[k] == 0); + xassert(len[k] == 0); + goto done; + } + /* there should be at least len[k] free locations */ + len_k = len[k]; + xassert(sva->r_ptr - sva->m_ptr >= len_k); + /* remove the vector from the linked list */ + if (prev[k] == 0) + sva->head = next[k]; + else + { /* preceding vector exists; increase its capacity */ + cap[prev[k]] += cap[k]; + next[prev[k]] = next[k]; + } + if (next[k] == 0) + sva->tail = prev[k]; + else + prev[next[k]] = prev[k]; + /* if the vector has zero length, make it empty */ + if (len_k == 0) + { ptr[k] = cap[k] = 0; + goto done; + } + /* copy the vector content to the beginning of the right part */ + ptr_k = sva->r_ptr - len_k; + memcpy(&ind[ptr_k], &ind[ptr[k]], len_k * sizeof(int)); + memcpy(&val[ptr_k], &val[ptr[k]], len_k * sizeof(double)); + /* set new pointer and capacity of the vector */ + ptr[k] = ptr_k; + cap[k] = len_k; + /* len[k] free locations have been consumed */ + sva->r_ptr -= len_k; +done: return; +} + +/*********************************************************************** +* sva_check_area - check sparse vector area (SVA) +* +* This routine checks the SVA data structures for correctness. +* +* NOTE: For testing/debugging only. */ + +void sva_check_area(SVA *sva) +{ int n_max = sva->n_max; + int n = sva->n; + int *ptr = sva->ptr; + int *len = sva->len; + int *cap = sva->cap; + int size = sva->size; + int m_ptr = sva->m_ptr; + int r_ptr = sva->r_ptr; + int head = sva->head; + int tail = sva->tail; + int *prev = sva->prev; + int *next = sva->next; + int k; +#if 0 /* 16/II-2004; SVA may be empty */ + xassert(1 <= n && n <= n_max); +#else + xassert(0 <= n && n <= n_max); +#endif + xassert(1 <= m_ptr && m_ptr <= r_ptr && r_ptr <= size+1); + /* all vectors included the linked list should have non-zero + * capacity and be stored in the left part */ + for (k = head; k != 0; k = next[k]) + { xassert(1 <= k && k <= n); + xassert(cap[k] > 0); + xassert(0 <= len[k] && len[k] <= cap[k]); + if (prev[k] == 0) + xassert(k == head); + else + { xassert(1 <= prev[k] && prev[k] <= n); + xassert(next[prev[k]] == k); + } + if (next[k] == 0) + { xassert(k == tail); + xassert(ptr[k] + cap[k] <= m_ptr); + } + else + { xassert(1 <= next[k] && next[k] <= n); + xassert(prev[next[k]] == k); + xassert(ptr[k] + cap[k] <= ptr[next[k]]); + } + cap[k] = -cap[k]; + } + /* all other vectors should either have zero capacity or be + * stored in the right part */ + for (k = 1; k <= n; k++) + { if (cap[k] < 0) + { /* k-th vector is stored in the left part */ + cap[k] = -cap[k]; + } + else if (cap[k] == 0) + { /* k-th vector has zero capacity */ + xassert(ptr[k] == 0); + xassert(len[k] == 0); + } + else /* cap[k] > 0 */ + { /* k-th vector is stored in the right part */ + xassert(0 <= len[k] && len[k] <= cap[k]); + xassert(r_ptr <= ptr[k] && ptr[k] + cap[k] <= size+1); + } + } + return; +} + +/*********************************************************************** +* sva_delete_area - delete sparse vector area (SVA) +* +* This routine deletes the sparse vector area (SVA) freeing all the +* memory allocated to it. */ + +void sva_delete_area(SVA *sva) +{ tfree(sva->ptr); + tfree(sva->len); + tfree(sva->cap); + tfree(sva->prev); + tfree(sva->next); + tfree(sva->ind); + tfree(sva->val); + tfree(sva); + return; +} + +/* eof */ diff --git a/glpk-5.0/src/bflib/sva.h b/glpk-5.0/src/bflib/sva.h new file mode 100644 index 0000000..2da8275 --- /dev/null +++ b/glpk-5.0/src/bflib/sva.h @@ -0,0 +1,159 @@ +/* sva.h (sparse vector area) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* Copyright (C) 2012-2013 Free Software Foundation, Inc. +* Written by Andrew Makhorin <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/>. +***********************************************************************/ + +#ifndef SVA_H +#define SVA_H + +/*********************************************************************** +* Sparse Vector Area (SVA) is a container for sparse vectors. This +* program object is used mainly on computing factorization, where the +* sparse vectors are rows and columns of sparse matrices. +* +* The SVA storage is a set of locations numbered 1, 2, ..., size, +* where size is the size of SVA, which is the total number of +* locations currently allocated. Each location is identified by its +* pointer p, 1 <= p <= size, and is the pair (ind[p], val[p]), where +* ind[p] and val[p] are, respectively, the index and value fields used +* to store the index and numeric value of a particular vector element. +* +* Each sparse vector is identified by its reference number k, +* 1 <= k <= n, where n is the total number of vectors currently stored +* in SVA, and defined by the triplet (ptr[k], len[k], cap[k]), where: +* ptr[k] is a pointer to the first location of the vector; len[k] is +* the vector length, which is the number of its non-zero elements, +* len[k] >= 0; and cap[k] is the capacity of the vector, which is the +* total number of adjacent locations allocated to that vector, +* cap[k] >= len[k]. Thus, non-zero elements of k-th vector are stored +* in locations ptr[k], ptr[k]+1, ..., ptr[k]+len[k]-1, and locations +* ptr[k]+len[k], ptr[k]+len[k]+1, ..., ptr[k]+cap[k]-1 are reserved. +* +* The SVA storage is divided into three parts as follows: +* +* Locations 1, 2, ..., m_ptr-1 constitute the left (dynamic) part of +* SVA. This part is used to store vectors, whose capacity may change. +* Note that all vectors stored in the left part are also included in +* a doubly linked list, where they are ordered by increasing their +* pointers ptr[k] (this list is needed for efficient implementation +* of the garbage collector used to defragment the left part of SVA); +* +* Locations m_ptr, m_ptr+1, ..., r_ptr-1 are free and constitute the +* middle (free) part of SVA. +* +* Locations r_ptr, r_ptr+1, ..., size constitute the right (static) +* part of SVA. This part is used to store vectors, whose capacity is +* not changed. */ + +typedef struct SVA SVA; + +struct SVA +{ /* sparse vector area */ + int n_max; + /* maximal value of n (enlarged automatically) */ + int n; + /* number of currently allocated vectors, 0 <= n <= n_max */ + int *ptr; /* int ptr[1+n_max]; */ + /* ptr[0] is not used; + * ptr[k], 1 <= i <= n, is pointer to first location of k-th + * vector in the arrays ind and val */ + int *len; /* int len[1+n_max]; */ + /* len[0] is not used; + * len[k], 1 <= k <= n, is length of k-th vector, len[k] >= 0 */ + int *cap; /* int cap[1+n_max]; */ + /* cap[0] is not used; + * cap[k], 1 <= k <= n, is capacity of k-th vector (the number + * of adjacent locations allocated to it), cap[k] >= len[k] */ + /* NOTE: if cap[k] = 0, then ptr[k] = 0 and len[k] = 0 */ + int size; + /* total number of locations in SVA */ + int m_ptr, r_ptr; + /* partitioning pointers that define the left, middle, and right + * parts of SVA (see above); 1 <= m_ptr <= r_ptr <= size+1 */ + int head; + /* number of first (leftmost) vector in the linked list */ + int tail; + /* number of last (rightmost) vector in the linked list */ + int *prev; /* int prev[1+n_max]; */ + /* prev[0] is not used; + * prev[k] is number of vector which precedes k-th vector in the + * linked list; + * prev[k] < 0 means that k-th vector is not in the list */ + int *next; /* int next[1+n_max]; */ + /* next[0] is not used; + * next[k] is number of vector which succedes k-th vector in the + * linked list; + * next[k] < 0 means that k-th vector is not in the list */ + /* NOTE: only vectors having non-zero capacity and stored in the + * left part of SVA are included in this linked list */ + int *ind; /* int ind[1+size]; */ + /* ind[0] is not used; + * ind[p], 1 <= p <= size, is index field of location p */ + double *val; /* double val[1+size]; */ + /* val[0] is not used; + * val[p], 1 <= p <= size, is value field of location p */ +#if 1 + int talky; + /* option to enable talky mode */ +#endif +}; + +#define sva_create_area _glp_sva_create_area +SVA *sva_create_area(int n_max, int size); +/* create sparse vector area (SVA) */ + +#define sva_alloc_vecs _glp_sva_alloc_vecs +int sva_alloc_vecs(SVA *sva, int nnn); +/* allocate new vectors in SVA */ + +#define sva_resize_area _glp_sva_resize_area +void sva_resize_area(SVA *sva, int delta); +/* change size of SVA storage */ + +#define sva_defrag_area _glp_sva_defrag_area +void sva_defrag_area(SVA *sva); +/* defragment left part of SVA */ + +#define sva_more_space _glp_sva_more_space +void sva_more_space(SVA *sva, int m_size); +/* increase size of middle (free) part of SVA */ + +#define sva_enlarge_cap _glp_sva_enlarge_cap +void sva_enlarge_cap(SVA *sva, int k, int new_cap, int skip); +/* enlarge capacity of specified vector */ + +#define sva_reserve_cap _glp_sva_reserve_cap +void sva_reserve_cap(SVA *sva, int k, int new_cap); +/* reserve locations for specified vector */ + +#define sva_make_static _glp_sva_make_static +void sva_make_static(SVA *sva, int k); +/* relocate specified vector to right part of SVA */ + +#define sva_check_area _glp_sva_check_area +void sva_check_area(SVA *sva); +/* check sparse vector area (SVA) */ + +#define sva_delete_area _glp_sva_delete_area +void sva_delete_area(SVA *sva); +/* delete sparse vector area (SVA) */ + +#endif + +/* eof */ |