diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/bflib/fhv.c')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/bflib/fhv.c | 586 |
1 files changed, 586 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/bflib/fhv.c b/test/monniaux/glpk-4.65/src/bflib/fhv.c new file mode 100644 index 00000000..e4bdf855 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/bflib/fhv.c @@ -0,0 +1,586 @@ +/* fhv.c (sparse updatable FHV-factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2012-2013 Andrew Makhorin, Department for Applied +* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights +* reserved. E-mail: <mao@gnu.org>. +* +* GLPK is free software: you can redistribute it and/or modify it +* under the terms of the GNU General Public License as published by +* the Free Software Foundation, either version 3 of the License, or +* (at your option) any later version. +* +* GLPK is distributed in the hope that it will be useful, but WITHOUT +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public +* License for more details. +* +* You should have received a copy of the GNU General Public License +* along with GLPK. If not, see <http://www.gnu.org/licenses/>. +***********************************************************************/ + +#include "env.h" +#include "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 */ |