1 /* fhvint.c (interface to FHV-factorization) */
2
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 * Copyright (C) 2012-2014 Free Software Foundation, Inc.
6 * Written by Andrew Makhorin <mao@gnu.org>.
7 *
8 * GLPK is free software: you can redistribute it and/or modify it
9 * under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * GLPK is distributed in the hope that it will be useful, but WITHOUT
14 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16 * License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
20 ***********************************************************************/
21
22 #include "env.h"
23 #include "fhvint.h"
24
fhvint_create(void)25 FHVINT *fhvint_create(void)
26 { /* create interface to FHV-factorization */
27 FHVINT *fi;
28 fi = talloc(1, FHVINT);
29 memset(fi, 0, sizeof(FHVINT));
30 fi->lufi = lufint_create();
31 return fi;
32 }
33
fhvint_factorize(FHVINT * fi,int n,int (* col)(void * info,int j,int ind[],double val[]),void * info)34 int fhvint_factorize(FHVINT *fi, int n, int (*col)(void *info, int j,
35 int ind[], double val[]), void *info)
36 { /* compute FHV-factorization of specified matrix A */
37 int nfs_max, old_n_max, n_max, k, ret;
38 xassert(n > 0);
39 fi->valid = 0;
40 /* get required value of nfs_max */
41 nfs_max = fi->nfs_max;
42 if (nfs_max == 0)
43 nfs_max = 100;
44 xassert(nfs_max > 0);
45 /* compute factorization of specified matrix A */
46 old_n_max = fi->lufi->n_max;
47 fi->lufi->sva_n_max = 4 * n + nfs_max;
48 fi->lufi->sgf_updat = 1;
49 ret = lufint_factorize(fi->lufi, n, col, info);
50 n_max = fi->lufi->n_max;
51 /* allocate/reallocate arrays, if necessary */
52 if (fi->fhv.nfs_max != nfs_max)
53 { if (fi->fhv.hh_ind != NULL)
54 tfree(fi->fhv.hh_ind);
55 fi->fhv.hh_ind = talloc(1+nfs_max, int);
56 }
57 if (old_n_max < n_max)
58 { if (fi->fhv.p0_ind != NULL)
59 tfree(fi->fhv.p0_ind);
60 if (fi->fhv.p0_inv != NULL)
61 tfree(fi->fhv.p0_inv);
62 fi->fhv.p0_ind = talloc(1+n_max, int);
63 fi->fhv.p0_inv = talloc(1+n_max, int);
64 }
65 /* initialize FHV-factorization */
66 fi->fhv.luf = fi->lufi->luf;
67 fi->fhv.nfs_max = nfs_max;
68 /* H := I */
69 fi->fhv.nfs = 0;
70 fi->fhv.hh_ref = sva_alloc_vecs(fi->lufi->sva, nfs_max);
71 /* P0 := P */
72 for (k = 1; k <= n; k++)
73 { fi->fhv.p0_ind[k] = fi->fhv.luf->pp_ind[k];
74 fi->fhv.p0_inv[k] = fi->fhv.luf->pp_inv[k];
75 }
76 /* set validation flag */
77 if (ret == 0)
78 fi->valid = 1;
79 return ret;
80 }
81
fhvint_update(FHVINT * fi,int j,int len,const int ind[],const double val[])82 int fhvint_update(FHVINT *fi, int j, int len, const int ind[],
83 const double val[])
84 { /* update FHV-factorization after replacing j-th column of A */
85 SGF *sgf = fi->lufi->sgf;
86 int *ind1 = sgf->rs_next;
87 double *val1 = sgf->vr_max;
88 double *work = sgf->work;
89 int ret;
90 xassert(fi->valid);
91 ret = fhv_ft_update(&fi->fhv, j, len, ind, val, ind1, val1, work);
92 if (ret != 0)
93 fi->valid = 0;
94 return ret;
95 }
96
fhvint_ftran(FHVINT * fi,double x[])97 void fhvint_ftran(FHVINT *fi, double x[])
98 { /* solve system A * x = b */
99 FHV *fhv = &fi->fhv;
100 LUF *luf = fhv->luf;
101 int n = luf->n;
102 int *pp_ind = luf->pp_ind;
103 int *pp_inv = luf->pp_inv;
104 SGF *sgf = fi->lufi->sgf;
105 double *work = sgf->work;
106 xassert(fi->valid);
107 /* A = F * H * V */
108 /* x = inv(A) * b = inv(V) * inv(H) * inv(F) * b */
109 luf->pp_ind = fhv->p0_ind;
110 luf->pp_inv = fhv->p0_inv;
111 luf_f_solve(luf, x);
112 luf->pp_ind = pp_ind;
113 luf->pp_inv = pp_inv;
114 fhv_h_solve(fhv, x);
115 luf_v_solve(luf, x, work);
116 memcpy(&x[1], &work[1], n * sizeof(double));
117 return;
118 }
119
fhvint_btran(FHVINT * fi,double x[])120 void fhvint_btran(FHVINT *fi, double x[])
121 { /* solve system A'* x = b */
122 FHV *fhv = &fi->fhv;
123 LUF *luf = fhv->luf;
124 int n = luf->n;
125 int *pp_ind = luf->pp_ind;
126 int *pp_inv = luf->pp_inv;
127 SGF *sgf = fi->lufi->sgf;
128 double *work = sgf->work;
129 xassert(fi->valid);
130 /* A' = (F * H * V)' = V'* H'* F' */
131 /* x = inv(A') * b = inv(F') * inv(H') * inv(V') * b */
132 luf_vt_solve(luf, x, work);
133 fhv_ht_solve(fhv, work);
134 luf->pp_ind = fhv->p0_ind;
135 luf->pp_inv = fhv->p0_inv;
136 luf_ft_solve(luf, work);
137 luf->pp_ind = pp_ind;
138 luf->pp_inv = pp_inv;
139 memcpy(&x[1], &work[1], n * sizeof(double));
140 return;
141 }
142
fhvint_estimate(FHVINT * fi)143 double fhvint_estimate(FHVINT *fi)
144 { /* estimate 1-norm of inv(A) */
145 double norm;
146 xassert(fi->valid);
147 xassert(fi->fhv.nfs == 0);
148 norm = luf_estimate_norm(fi->fhv.luf, fi->lufi->sgf->vr_max,
149 fi->lufi->sgf->work);
150 return norm;
151 }
152
fhvint_delete(FHVINT * fi)153 void fhvint_delete(FHVINT *fi)
154 { /* delete interface to FHV-factorization */
155 lufint_delete(fi->lufi);
156 if (fi->fhv.hh_ind != NULL)
157 tfree(fi->fhv.hh_ind);
158 if (fi->fhv.p0_ind != NULL)
159 tfree(fi->fhv.p0_ind);
160 if (fi->fhv.p0_inv != NULL)
161 tfree(fi->fhv.p0_inv);
162 tfree(fi);
163 return;
164 }
165
166 /* eof */
167