1 | /* glpfhv.c (LP basis factorization, FHV eta file version) */ |
---|
2 | |
---|
3 | /*********************************************************************** |
---|
4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
5 | * |
---|
6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
7 | * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
---|
8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
9 | * E-mail: <mao@gnu.org>. |
---|
10 | * |
---|
11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
12 | * under the terms of the GNU General Public License as published by |
---|
13 | * the Free Software Foundation, either version 3 of the License, or |
---|
14 | * (at your option) any later version. |
---|
15 | * |
---|
16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
19 | * License for more details. |
---|
20 | * |
---|
21 | * You should have received a copy of the GNU General Public License |
---|
22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
23 | ***********************************************************************/ |
---|
24 | |
---|
25 | #include "glpfhv.h" |
---|
26 | #include "glpenv.h" |
---|
27 | #define xfault xerror |
---|
28 | |
---|
29 | /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ |
---|
30 | |
---|
31 | #define M_MAX 100000000 /* = 100*10^6 */ |
---|
32 | /* maximal order of the basis matrix */ |
---|
33 | |
---|
34 | /*********************************************************************** |
---|
35 | * NAME |
---|
36 | * |
---|
37 | * fhv_create_it - create LP basis factorization |
---|
38 | * |
---|
39 | * SYNOPSIS |
---|
40 | * |
---|
41 | * #include "glpfhv.h" |
---|
42 | * FHV *fhv_create_it(void); |
---|
43 | * |
---|
44 | * DESCRIPTION |
---|
45 | * |
---|
46 | * The routine fhv_create_it creates a program object, which represents |
---|
47 | * a factorization of LP basis. |
---|
48 | * |
---|
49 | * RETURNS |
---|
50 | * |
---|
51 | * The routine fhv_create_it returns a pointer to the object created. */ |
---|
52 | |
---|
53 | FHV *fhv_create_it(void) |
---|
54 | { FHV *fhv; |
---|
55 | fhv = xmalloc(sizeof(FHV)); |
---|
56 | fhv->m_max = fhv->m = 0; |
---|
57 | fhv->valid = 0; |
---|
58 | fhv->luf = luf_create_it(); |
---|
59 | fhv->hh_max = 50; |
---|
60 | fhv->hh_nfs = 0; |
---|
61 | fhv->hh_ind = fhv->hh_ptr = fhv->hh_len = NULL; |
---|
62 | fhv->p0_row = fhv->p0_col = NULL; |
---|
63 | fhv->cc_ind = NULL; |
---|
64 | fhv->cc_val = NULL; |
---|
65 | fhv->upd_tol = 1e-6; |
---|
66 | fhv->nnz_h = 0; |
---|
67 | return fhv; |
---|
68 | } |
---|
69 | |
---|
70 | /*********************************************************************** |
---|
71 | * NAME |
---|
72 | * |
---|
73 | * fhv_factorize - compute LP basis factorization |
---|
74 | * |
---|
75 | * SYNOPSIS |
---|
76 | * |
---|
77 | * #include "glpfhv.h" |
---|
78 | * int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j, |
---|
79 | * int ind[], double val[]), void *info); |
---|
80 | * |
---|
81 | * DESCRIPTION |
---|
82 | * |
---|
83 | * The routine fhv_factorize computes the factorization of the basis |
---|
84 | * matrix B specified by the routine col. |
---|
85 | * |
---|
86 | * The parameter fhv specified the basis factorization data structure |
---|
87 | * created by the routine fhv_create_it. |
---|
88 | * |
---|
89 | * The parameter m specifies the order of B, m > 0. |
---|
90 | * |
---|
91 | * The formal routine col specifies the matrix B to be factorized. To |
---|
92 | * obtain j-th column of A the routine fhv_factorize calls the routine |
---|
93 | * col with the parameter j (1 <= j <= n). In response the routine col |
---|
94 | * should store row indices and numerical values of non-zero elements |
---|
95 | * of j-th column of B to locations ind[1,...,len] and val[1,...,len], |
---|
96 | * respectively, where len is the number of non-zeros in j-th column |
---|
97 | * returned on exit. Neither zero nor duplicate elements are allowed. |
---|
98 | * |
---|
99 | * The parameter info is a transit pointer passed to the routine col. |
---|
100 | * |
---|
101 | * RETURNS |
---|
102 | * |
---|
103 | * 0 The factorization has been successfully computed. |
---|
104 | * |
---|
105 | * FHV_ESING |
---|
106 | * The specified matrix is singular within the working precision. |
---|
107 | * |
---|
108 | * FHV_ECOND |
---|
109 | * The specified matrix is ill-conditioned. |
---|
110 | * |
---|
111 | * For more details see comments to the routine luf_factorize. |
---|
112 | * |
---|
113 | * ALGORITHM |
---|
114 | * |
---|
115 | * The routine fhv_factorize calls the routine luf_factorize (see the |
---|
116 | * module GLPLUF), which actually computes LU-factorization of the basis |
---|
117 | * matrix B in the form |
---|
118 | * |
---|
119 | * [B] = (F, V, P, Q), |
---|
120 | * |
---|
121 | * where F and V are such matrices that |
---|
122 | * |
---|
123 | * B = F * V, |
---|
124 | * |
---|
125 | * and P and Q are such permutation matrices that the matrix |
---|
126 | * |
---|
127 | * L = P * F * inv(P) |
---|
128 | * |
---|
129 | * is lower triangular with unity diagonal, and the matrix |
---|
130 | * |
---|
131 | * U = P * V * Q |
---|
132 | * |
---|
133 | * is upper triangular. |
---|
134 | * |
---|
135 | * In order to build the complete representation of the factorization |
---|
136 | * (see formula (1) in the file glpfhv.h) the routine fhv_factorize just |
---|
137 | * additionally sets H = I and P0 = P. */ |
---|
138 | |
---|
139 | int fhv_factorize(FHV *fhv, int m, int (*col)(void *info, int j, |
---|
140 | int ind[], double val[]), void *info) |
---|
141 | { int ret; |
---|
142 | if (m < 1) |
---|
143 | xfault("fhv_factorize: m = %d; invalid parameter\n", m); |
---|
144 | if (m > M_MAX) |
---|
145 | xfault("fhv_factorize: m = %d; matrix too big\n", m); |
---|
146 | fhv->m = m; |
---|
147 | /* invalidate the factorization */ |
---|
148 | fhv->valid = 0; |
---|
149 | /* allocate/reallocate arrays, if necessary */ |
---|
150 | if (fhv->hh_ind == NULL) |
---|
151 | fhv->hh_ind = xcalloc(1+fhv->hh_max, sizeof(int)); |
---|
152 | if (fhv->hh_ptr == NULL) |
---|
153 | fhv->hh_ptr = xcalloc(1+fhv->hh_max, sizeof(int)); |
---|
154 | if (fhv->hh_len == NULL) |
---|
155 | fhv->hh_len = xcalloc(1+fhv->hh_max, sizeof(int)); |
---|
156 | if (fhv->m_max < m) |
---|
157 | { if (fhv->p0_row != NULL) xfree(fhv->p0_row); |
---|
158 | if (fhv->p0_col != NULL) xfree(fhv->p0_col); |
---|
159 | if (fhv->cc_ind != NULL) xfree(fhv->cc_ind); |
---|
160 | if (fhv->cc_val != NULL) xfree(fhv->cc_val); |
---|
161 | fhv->m_max = m + 100; |
---|
162 | fhv->p0_row = xcalloc(1+fhv->m_max, sizeof(int)); |
---|
163 | fhv->p0_col = xcalloc(1+fhv->m_max, sizeof(int)); |
---|
164 | fhv->cc_ind = xcalloc(1+fhv->m_max, sizeof(int)); |
---|
165 | fhv->cc_val = xcalloc(1+fhv->m_max, sizeof(double)); |
---|
166 | } |
---|
167 | /* try to factorize the basis matrix */ |
---|
168 | switch (luf_factorize(fhv->luf, m, col, info)) |
---|
169 | { case 0: |
---|
170 | break; |
---|
171 | case LUF_ESING: |
---|
172 | ret = FHV_ESING; |
---|
173 | goto done; |
---|
174 | case LUF_ECOND: |
---|
175 | ret = FHV_ECOND; |
---|
176 | goto done; |
---|
177 | default: |
---|
178 | xassert(fhv != fhv); |
---|
179 | } |
---|
180 | /* the basis matrix has been successfully factorized */ |
---|
181 | fhv->valid = 1; |
---|
182 | /* H := I */ |
---|
183 | fhv->hh_nfs = 0; |
---|
184 | /* P0 := P */ |
---|
185 | memcpy(&fhv->p0_row[1], &fhv->luf->pp_row[1], sizeof(int) * m); |
---|
186 | memcpy(&fhv->p0_col[1], &fhv->luf->pp_col[1], sizeof(int) * m); |
---|
187 | /* currently H has no factors */ |
---|
188 | fhv->nnz_h = 0; |
---|
189 | ret = 0; |
---|
190 | done: /* return to the calling program */ |
---|
191 | return ret; |
---|
192 | } |
---|
193 | |
---|
194 | /*********************************************************************** |
---|
195 | * NAME |
---|
196 | * |
---|
197 | * fhv_h_solve - solve system H*x = b or H'*x = b |
---|
198 | * |
---|
199 | * SYNOPSIS |
---|
200 | * |
---|
201 | * #include "glpfhv.h" |
---|
202 | * void fhv_h_solve(FHV *fhv, int tr, double x[]); |
---|
203 | * |
---|
204 | * DESCRIPTION |
---|
205 | * |
---|
206 | * The routine fhv_h_solve solves either the system H*x = b (if the |
---|
207 | * flag tr is zero) or the system H'*x = b (if the flag tr is non-zero), |
---|
208 | * where the matrix H is a component of the factorization specified by |
---|
209 | * the parameter fhv, H' is a matrix transposed to H. |
---|
210 | * |
---|
211 | * On entry the array x should contain elements of the right-hand side |
---|
212 | * vector b in locations x[1], ..., x[m], where m is the order of the |
---|
213 | * matrix H. On exit this array will contain elements of the solution |
---|
214 | * vector x in the same locations. */ |
---|
215 | |
---|
216 | void fhv_h_solve(FHV *fhv, int tr, double x[]) |
---|
217 | { int nfs = fhv->hh_nfs; |
---|
218 | int *hh_ind = fhv->hh_ind; |
---|
219 | int *hh_ptr = fhv->hh_ptr; |
---|
220 | int *hh_len = fhv->hh_len; |
---|
221 | int *sv_ind = fhv->luf->sv_ind; |
---|
222 | double *sv_val = fhv->luf->sv_val; |
---|
223 | int i, k, beg, end, ptr; |
---|
224 | double temp; |
---|
225 | if (!fhv->valid) |
---|
226 | xfault("fhv_h_solve: the factorization is not valid\n"); |
---|
227 | if (!tr) |
---|
228 | { /* solve the system H*x = b */ |
---|
229 | for (k = 1; k <= nfs; k++) |
---|
230 | { i = hh_ind[k]; |
---|
231 | temp = x[i]; |
---|
232 | beg = hh_ptr[k]; |
---|
233 | end = beg + hh_len[k] - 1; |
---|
234 | for (ptr = beg; ptr <= end; ptr++) |
---|
235 | temp -= sv_val[ptr] * x[sv_ind[ptr]]; |
---|
236 | x[i] = temp; |
---|
237 | } |
---|
238 | } |
---|
239 | else |
---|
240 | { /* solve the system H'*x = b */ |
---|
241 | for (k = nfs; k >= 1; k--) |
---|
242 | { i = hh_ind[k]; |
---|
243 | temp = x[i]; |
---|
244 | if (temp == 0.0) continue; |
---|
245 | beg = hh_ptr[k]; |
---|
246 | end = beg + hh_len[k] - 1; |
---|
247 | for (ptr = beg; ptr <= end; ptr++) |
---|
248 | x[sv_ind[ptr]] -= sv_val[ptr] * temp; |
---|
249 | } |
---|
250 | } |
---|
251 | return; |
---|
252 | } |
---|
253 | |
---|
254 | /*********************************************************************** |
---|
255 | * NAME |
---|
256 | * |
---|
257 | * fhv_ftran - perform forward transformation (solve system B*x = b) |
---|
258 | * |
---|
259 | * SYNOPSIS |
---|
260 | * |
---|
261 | * #include "glpfhv.h" |
---|
262 | * void fhv_ftran(FHV *fhv, double x[]); |
---|
263 | * |
---|
264 | * DESCRIPTION |
---|
265 | * |
---|
266 | * The routine fhv_ftran performs forward transformation, i.e. solves |
---|
267 | * the system B*x = b, where B is the basis matrix, x is the vector of |
---|
268 | * unknowns to be computed, b is the vector of right-hand sides. |
---|
269 | * |
---|
270 | * On entry elements of the vector b should be stored in dense format |
---|
271 | * in locations x[1], ..., x[m], where m is the number of rows. On exit |
---|
272 | * the routine stores elements of the vector x in the same locations. */ |
---|
273 | |
---|
274 | void fhv_ftran(FHV *fhv, double x[]) |
---|
275 | { int *pp_row = fhv->luf->pp_row; |
---|
276 | int *pp_col = fhv->luf->pp_col; |
---|
277 | int *p0_row = fhv->p0_row; |
---|
278 | int *p0_col = fhv->p0_col; |
---|
279 | if (!fhv->valid) |
---|
280 | xfault("fhv_ftran: the factorization is not valid\n"); |
---|
281 | /* B = F*H*V, therefore inv(B) = inv(V)*inv(H)*inv(F) */ |
---|
282 | fhv->luf->pp_row = p0_row; |
---|
283 | fhv->luf->pp_col = p0_col; |
---|
284 | luf_f_solve(fhv->luf, 0, x); |
---|
285 | fhv->luf->pp_row = pp_row; |
---|
286 | fhv->luf->pp_col = pp_col; |
---|
287 | fhv_h_solve(fhv, 0, x); |
---|
288 | luf_v_solve(fhv->luf, 0, x); |
---|
289 | return; |
---|
290 | } |
---|
291 | |
---|
292 | /*********************************************************************** |
---|
293 | * NAME |
---|
294 | * |
---|
295 | * fhv_btran - perform backward transformation (solve system B'*x = b) |
---|
296 | * |
---|
297 | * SYNOPSIS |
---|
298 | * |
---|
299 | * #include "glpfhv.h" |
---|
300 | * void fhv_btran(FHV *fhv, double x[]); |
---|
301 | * |
---|
302 | * DESCRIPTION |
---|
303 | * |
---|
304 | * The routine fhv_btran performs backward transformation, i.e. solves |
---|
305 | * the system B'*x = b, where B' is a matrix transposed to the basis |
---|
306 | * matrix B, x is the vector of unknowns to be computed, b is the vector |
---|
307 | * of right-hand sides. |
---|
308 | * |
---|
309 | * On entry elements of the vector b should be stored in dense format |
---|
310 | * in locations x[1], ..., x[m], where m is the number of rows. On exit |
---|
311 | * the routine stores elements of the vector x in the same locations. */ |
---|
312 | |
---|
313 | void fhv_btran(FHV *fhv, double x[]) |
---|
314 | { int *pp_row = fhv->luf->pp_row; |
---|
315 | int *pp_col = fhv->luf->pp_col; |
---|
316 | int *p0_row = fhv->p0_row; |
---|
317 | int *p0_col = fhv->p0_col; |
---|
318 | if (!fhv->valid) |
---|
319 | xfault("fhv_btran: the factorization is not valid\n"); |
---|
320 | /* B = F*H*V, therefore inv(B') = inv(F')*inv(H')*inv(V') */ |
---|
321 | luf_v_solve(fhv->luf, 1, x); |
---|
322 | fhv_h_solve(fhv, 1, x); |
---|
323 | fhv->luf->pp_row = p0_row; |
---|
324 | fhv->luf->pp_col = p0_col; |
---|
325 | luf_f_solve(fhv->luf, 1, x); |
---|
326 | fhv->luf->pp_row = pp_row; |
---|
327 | fhv->luf->pp_col = pp_col; |
---|
328 | return; |
---|
329 | } |
---|
330 | |
---|
331 | /*********************************************************************** |
---|
332 | * NAME |
---|
333 | * |
---|
334 | * fhv_update_it - update LP basis factorization |
---|
335 | * |
---|
336 | * SYNOPSIS |
---|
337 | * |
---|
338 | * #include "glpfhv.h" |
---|
339 | * int fhv_update_it(FHV *fhv, int j, int len, const int ind[], |
---|
340 | * const double val[]); |
---|
341 | * |
---|
342 | * DESCRIPTION |
---|
343 | * |
---|
344 | * The routine fhv_update_it updates the factorization of the basis |
---|
345 | * matrix B after replacing its j-th column by a new vector. |
---|
346 | * |
---|
347 | * The parameter j specifies the number of column of B, which has been |
---|
348 | * replaced, 1 <= j <= m, where m is the order of B. |
---|
349 | * |
---|
350 | * Row indices and numerical values of non-zero elements of the new |
---|
351 | * column of B should be placed in locations ind[1], ..., ind[len] and |
---|
352 | * val[1], ..., val[len], resp., where len is the number of non-zeros |
---|
353 | * in the column. Neither zero nor duplicate elements are allowed. |
---|
354 | * |
---|
355 | * RETURNS |
---|
356 | * |
---|
357 | * 0 The factorization has been successfully updated. |
---|
358 | * |
---|
359 | * FHV_ESING |
---|
360 | * The adjacent basis matrix is structurally singular, since after |
---|
361 | * changing j-th column of matrix V by the new column (see algorithm |
---|
362 | * below) the case k1 > k2 occured. |
---|
363 | * |
---|
364 | * FHV_ECHECK |
---|
365 | * The factorization is inaccurate, since after transforming k2-th |
---|
366 | * row of matrix U = P*V*Q, its diagonal element u[k2,k2] is zero or |
---|
367 | * close to zero, |
---|
368 | * |
---|
369 | * FHV_ELIMIT |
---|
370 | * Maximal number of H factors has been reached. |
---|
371 | * |
---|
372 | * FHV_EROOM |
---|
373 | * Overflow of the sparse vector area. |
---|
374 | * |
---|
375 | * In case of non-zero return code the factorization becomes invalid. |
---|
376 | * It should not be used until it has been recomputed with the routine |
---|
377 | * fhv_factorize. |
---|
378 | * |
---|
379 | * ALGORITHM |
---|
380 | * |
---|
381 | * The routine fhv_update_it is based on the transformation proposed by |
---|
382 | * Forrest and Tomlin. |
---|
383 | * |
---|
384 | * Let j-th column of the basis matrix B have been replaced by new |
---|
385 | * column B[j]. In order to keep the equality B = F*H*V j-th column of |
---|
386 | * matrix V should be replaced by the column inv(F*H)*B[j]. |
---|
387 | * |
---|
388 | * From the standpoint of matrix U = P*V*Q, replacement of j-th column |
---|
389 | * of matrix V is equivalent to replacement of k1-th column of matrix U, |
---|
390 | * where k1 is determined by permutation matrix Q. Thus, matrix U loses |
---|
391 | * its upper triangular form and becomes the following: |
---|
392 | * |
---|
393 | * 1 k1 k2 m |
---|
394 | * 1 x x * x x x x x x x |
---|
395 | * . x * x x x x x x x |
---|
396 | * k1 . . * x x x x x x x |
---|
397 | * . . * x x x x x x x |
---|
398 | * . . * . x x x x x x |
---|
399 | * . . * . . x x x x x |
---|
400 | * . . * . . . x x x x |
---|
401 | * k2 . . * . . . . x x x |
---|
402 | * . . . . . . . . x x |
---|
403 | * m . . . . . . . . . x |
---|
404 | * |
---|
405 | * where row index k2 corresponds to the lowest non-zero element of |
---|
406 | * k1-th column. |
---|
407 | * |
---|
408 | * The routine moves rows and columns k1+1, k1+2, ..., k2 of matrix U |
---|
409 | * by one position to the left and upwards and moves k1-th row and k1-th |
---|
410 | * column to position k2. As the result of such symmetric permutations |
---|
411 | * matrix U becomes the following: |
---|
412 | * |
---|
413 | * 1 k1 k2 m |
---|
414 | * 1 x x x x x x x * x x |
---|
415 | * . x x x x x x * x x |
---|
416 | * k1 . . x x x x x * x x |
---|
417 | * . . . x x x x * x x |
---|
418 | * . . . . x x x * x x |
---|
419 | * . . . . . x x * x x |
---|
420 | * . . . . . . x * x x |
---|
421 | * k2 . . x x x x x * x x |
---|
422 | * . . . . . . . . x x |
---|
423 | * m . . . . . . . . . x |
---|
424 | * |
---|
425 | * Then the routine performs gaussian elimination to eliminate elements |
---|
426 | * u[k2,k1], u[k2,k1+1], ..., u[k2,k2-1] using diagonal elements |
---|
427 | * u[k1,k1], u[k1+1,k1+1], ..., u[k2-1,k2-1] as pivots in the same way |
---|
428 | * as described in comments to the routine luf_factorize (see the module |
---|
429 | * GLPLUF). Note that actually all operations are performed on matrix V, |
---|
430 | * not on matrix U. During the elimination process the routine permutes |
---|
431 | * neither rows nor columns, so only k2-th row of matrix U is changed. |
---|
432 | * |
---|
433 | * To keep the main equality B = F*H*V, each time when the routine |
---|
434 | * applies elementary gaussian transformation to the transformed row of |
---|
435 | * matrix V (which corresponds to k2-th row of matrix U), it also adds |
---|
436 | * a new element (gaussian multiplier) to the current row-like factor |
---|
437 | * of matrix H, which corresponds to the transformed row of matrix V. */ |
---|
438 | |
---|
439 | int fhv_update_it(FHV *fhv, int j, int len, const int ind[], |
---|
440 | const double val[]) |
---|
441 | { int m = fhv->m; |
---|
442 | LUF *luf = fhv->luf; |
---|
443 | int *vr_ptr = luf->vr_ptr; |
---|
444 | int *vr_len = luf->vr_len; |
---|
445 | int *vr_cap = luf->vr_cap; |
---|
446 | double *vr_piv = luf->vr_piv; |
---|
447 | int *vc_ptr = luf->vc_ptr; |
---|
448 | int *vc_len = luf->vc_len; |
---|
449 | int *vc_cap = luf->vc_cap; |
---|
450 | int *pp_row = luf->pp_row; |
---|
451 | int *pp_col = luf->pp_col; |
---|
452 | int *qq_row = luf->qq_row; |
---|
453 | int *qq_col = luf->qq_col; |
---|
454 | int *sv_ind = luf->sv_ind; |
---|
455 | double *sv_val = luf->sv_val; |
---|
456 | double *work = luf->work; |
---|
457 | double eps_tol = luf->eps_tol; |
---|
458 | int *hh_ind = fhv->hh_ind; |
---|
459 | int *hh_ptr = fhv->hh_ptr; |
---|
460 | int *hh_len = fhv->hh_len; |
---|
461 | int *p0_row = fhv->p0_row; |
---|
462 | int *p0_col = fhv->p0_col; |
---|
463 | int *cc_ind = fhv->cc_ind; |
---|
464 | double *cc_val = fhv->cc_val; |
---|
465 | double upd_tol = fhv->upd_tol; |
---|
466 | int i, i_beg, i_end, i_ptr, j_beg, j_end, j_ptr, k, k1, k2, p, q, |
---|
467 | p_beg, p_end, p_ptr, ptr, ret; |
---|
468 | double f, temp; |
---|
469 | if (!fhv->valid) |
---|
470 | xfault("fhv_update_it: the factorization is not valid\n"); |
---|
471 | if (!(1 <= j && j <= m)) |
---|
472 | xfault("fhv_update_it: j = %d; column number out of range\n", |
---|
473 | j); |
---|
474 | /* check if the new factor of matrix H can be created */ |
---|
475 | if (fhv->hh_nfs == fhv->hh_max) |
---|
476 | { /* maximal number of updates has been reached */ |
---|
477 | fhv->valid = 0; |
---|
478 | ret = FHV_ELIMIT; |
---|
479 | goto done; |
---|
480 | } |
---|
481 | /* convert new j-th column of B to dense format */ |
---|
482 | for (i = 1; i <= m; i++) |
---|
483 | cc_val[i] = 0.0; |
---|
484 | for (k = 1; k <= len; k++) |
---|
485 | { i = ind[k]; |
---|
486 | if (!(1 <= i && i <= m)) |
---|
487 | xfault("fhv_update_it: ind[%d] = %d; row number out of rang" |
---|
488 | "e\n", k, i); |
---|
489 | if (cc_val[i] != 0.0) |
---|
490 | xfault("fhv_update_it: ind[%d] = %d; duplicate row index no" |
---|
491 | "t allowed\n", k, i); |
---|
492 | if (val[k] == 0.0) |
---|
493 | xfault("fhv_update_it: val[%d] = %g; zero element not allow" |
---|
494 | "ed\n", k, val[k]); |
---|
495 | cc_val[i] = val[k]; |
---|
496 | } |
---|
497 | /* new j-th column of V := inv(F * H) * (new B[j]) */ |
---|
498 | fhv->luf->pp_row = p0_row; |
---|
499 | fhv->luf->pp_col = p0_col; |
---|
500 | luf_f_solve(fhv->luf, 0, cc_val); |
---|
501 | fhv->luf->pp_row = pp_row; |
---|
502 | fhv->luf->pp_col = pp_col; |
---|
503 | fhv_h_solve(fhv, 0, cc_val); |
---|
504 | /* convert new j-th column of V to sparse format */ |
---|
505 | len = 0; |
---|
506 | for (i = 1; i <= m; i++) |
---|
507 | { temp = cc_val[i]; |
---|
508 | if (temp == 0.0 || fabs(temp) < eps_tol) continue; |
---|
509 | len++, cc_ind[len] = i, cc_val[len] = temp; |
---|
510 | } |
---|
511 | /* clear old content of j-th column of matrix V */ |
---|
512 | j_beg = vc_ptr[j]; |
---|
513 | j_end = j_beg + vc_len[j] - 1; |
---|
514 | for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) |
---|
515 | { /* get row index of v[i,j] */ |
---|
516 | i = sv_ind[j_ptr]; |
---|
517 | /* find v[i,j] in the i-th row */ |
---|
518 | i_beg = vr_ptr[i]; |
---|
519 | i_end = i_beg + vr_len[i] - 1; |
---|
520 | for (i_ptr = i_beg; sv_ind[i_ptr] != j; i_ptr++) /* nop */; |
---|
521 | xassert(i_ptr <= i_end); |
---|
522 | /* remove v[i,j] from the i-th row */ |
---|
523 | sv_ind[i_ptr] = sv_ind[i_end]; |
---|
524 | sv_val[i_ptr] = sv_val[i_end]; |
---|
525 | vr_len[i]--; |
---|
526 | } |
---|
527 | /* now j-th column of matrix V is empty */ |
---|
528 | luf->nnz_v -= vc_len[j]; |
---|
529 | vc_len[j] = 0; |
---|
530 | /* add new elements of j-th column of matrix V to corresponding |
---|
531 | row lists; determine indices k1 and k2 */ |
---|
532 | k1 = qq_row[j], k2 = 0; |
---|
533 | for (ptr = 1; ptr <= len; ptr++) |
---|
534 | { /* get row index of v[i,j] */ |
---|
535 | i = cc_ind[ptr]; |
---|
536 | /* at least one unused location is needed in i-th row */ |
---|
537 | if (vr_len[i] + 1 > vr_cap[i]) |
---|
538 | { if (luf_enlarge_row(luf, i, vr_len[i] + 10)) |
---|
539 | { /* overflow of the sparse vector area */ |
---|
540 | fhv->valid = 0; |
---|
541 | luf->new_sva = luf->sv_size + luf->sv_size; |
---|
542 | xassert(luf->new_sva > luf->sv_size); |
---|
543 | ret = FHV_EROOM; |
---|
544 | goto done; |
---|
545 | } |
---|
546 | } |
---|
547 | /* add v[i,j] to i-th row */ |
---|
548 | i_ptr = vr_ptr[i] + vr_len[i]; |
---|
549 | sv_ind[i_ptr] = j; |
---|
550 | sv_val[i_ptr] = cc_val[ptr]; |
---|
551 | vr_len[i]++; |
---|
552 | /* adjust index k2 */ |
---|
553 | if (k2 < pp_col[i]) k2 = pp_col[i]; |
---|
554 | } |
---|
555 | /* capacity of j-th column (which is currently empty) should be |
---|
556 | not less than len locations */ |
---|
557 | if (vc_cap[j] < len) |
---|
558 | { if (luf_enlarge_col(luf, j, len)) |
---|
559 | { /* overflow of the sparse vector area */ |
---|
560 | fhv->valid = 0; |
---|
561 | luf->new_sva = luf->sv_size + luf->sv_size; |
---|
562 | xassert(luf->new_sva > luf->sv_size); |
---|
563 | ret = FHV_EROOM; |
---|
564 | goto done; |
---|
565 | } |
---|
566 | } |
---|
567 | /* add new elements of matrix V to j-th column list */ |
---|
568 | j_ptr = vc_ptr[j]; |
---|
569 | memmove(&sv_ind[j_ptr], &cc_ind[1], len * sizeof(int)); |
---|
570 | memmove(&sv_val[j_ptr], &cc_val[1], len * sizeof(double)); |
---|
571 | vc_len[j] = len; |
---|
572 | luf->nnz_v += len; |
---|
573 | /* if k1 > k2, diagonal element u[k2,k2] of matrix U is zero and |
---|
574 | therefore the adjacent basis matrix is structurally singular */ |
---|
575 | if (k1 > k2) |
---|
576 | { fhv->valid = 0; |
---|
577 | ret = FHV_ESING; |
---|
578 | goto done; |
---|
579 | } |
---|
580 | /* perform implicit symmetric permutations of rows and columns of |
---|
581 | matrix U */ |
---|
582 | i = pp_row[k1], j = qq_col[k1]; |
---|
583 | for (k = k1; k < k2; k++) |
---|
584 | { pp_row[k] = pp_row[k+1], pp_col[pp_row[k]] = k; |
---|
585 | qq_col[k] = qq_col[k+1], qq_row[qq_col[k]] = k; |
---|
586 | } |
---|
587 | pp_row[k2] = i, pp_col[i] = k2; |
---|
588 | qq_col[k2] = j, qq_row[j] = k2; |
---|
589 | /* now i-th row of the matrix V is k2-th row of matrix U; since |
---|
590 | no pivoting is used, only this row will be transformed */ |
---|
591 | /* copy elements of i-th row of matrix V to the working array and |
---|
592 | remove these elements from matrix V */ |
---|
593 | for (j = 1; j <= m; j++) work[j] = 0.0; |
---|
594 | i_beg = vr_ptr[i]; |
---|
595 | i_end = i_beg + vr_len[i] - 1; |
---|
596 | for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) |
---|
597 | { /* get column index of v[i,j] */ |
---|
598 | j = sv_ind[i_ptr]; |
---|
599 | /* store v[i,j] to the working array */ |
---|
600 | work[j] = sv_val[i_ptr]; |
---|
601 | /* find v[i,j] in the j-th column */ |
---|
602 | j_beg = vc_ptr[j]; |
---|
603 | j_end = j_beg + vc_len[j] - 1; |
---|
604 | for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++) /* nop */; |
---|
605 | xassert(j_ptr <= j_end); |
---|
606 | /* remove v[i,j] from the j-th column */ |
---|
607 | sv_ind[j_ptr] = sv_ind[j_end]; |
---|
608 | sv_val[j_ptr] = sv_val[j_end]; |
---|
609 | vc_len[j]--; |
---|
610 | } |
---|
611 | /* now i-th row of matrix V is empty */ |
---|
612 | luf->nnz_v -= vr_len[i]; |
---|
613 | vr_len[i] = 0; |
---|
614 | /* create the next row-like factor of the matrix H; this factor |
---|
615 | corresponds to i-th (transformed) row */ |
---|
616 | fhv->hh_nfs++; |
---|
617 | hh_ind[fhv->hh_nfs] = i; |
---|
618 | /* hh_ptr[] will be set later */ |
---|
619 | hh_len[fhv->hh_nfs] = 0; |
---|
620 | /* up to (k2 - k1) free locations are needed to add new elements |
---|
621 | to the non-trivial row of the row-like factor */ |
---|
622 | if (luf->sv_end - luf->sv_beg < k2 - k1) |
---|
623 | { luf_defrag_sva(luf); |
---|
624 | if (luf->sv_end - luf->sv_beg < k2 - k1) |
---|
625 | { /* overflow of the sparse vector area */ |
---|
626 | fhv->valid = luf->valid = 0; |
---|
627 | luf->new_sva = luf->sv_size + luf->sv_size; |
---|
628 | xassert(luf->new_sva > luf->sv_size); |
---|
629 | ret = FHV_EROOM; |
---|
630 | goto done; |
---|
631 | } |
---|
632 | } |
---|
633 | /* eliminate subdiagonal elements of matrix U */ |
---|
634 | for (k = k1; k < k2; k++) |
---|
635 | { /* v[p,q] = u[k,k] */ |
---|
636 | p = pp_row[k], q = qq_col[k]; |
---|
637 | /* this is the crucial point, where even tiny non-zeros should |
---|
638 | not be dropped */ |
---|
639 | if (work[q] == 0.0) continue; |
---|
640 | /* compute gaussian multiplier f = v[i,q] / v[p,q] */ |
---|
641 | f = work[q] / vr_piv[p]; |
---|
642 | /* perform gaussian transformation: |
---|
643 | (i-th row) := (i-th row) - f * (p-th row) |
---|
644 | in order to eliminate v[i,q] = u[k2,k] */ |
---|
645 | p_beg = vr_ptr[p]; |
---|
646 | p_end = p_beg + vr_len[p] - 1; |
---|
647 | for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++) |
---|
648 | work[sv_ind[p_ptr]] -= f * sv_val[p_ptr]; |
---|
649 | /* store new element (gaussian multiplier that corresponds to |
---|
650 | p-th row) in the current row-like factor */ |
---|
651 | luf->sv_end--; |
---|
652 | sv_ind[luf->sv_end] = p; |
---|
653 | sv_val[luf->sv_end] = f; |
---|
654 | hh_len[fhv->hh_nfs]++; |
---|
655 | } |
---|
656 | /* set pointer to the current row-like factor of the matrix H |
---|
657 | (if no elements were added to this factor, it is unity matrix |
---|
658 | and therefore can be discarded) */ |
---|
659 | if (hh_len[fhv->hh_nfs] == 0) |
---|
660 | fhv->hh_nfs--; |
---|
661 | else |
---|
662 | { hh_ptr[fhv->hh_nfs] = luf->sv_end; |
---|
663 | fhv->nnz_h += hh_len[fhv->hh_nfs]; |
---|
664 | } |
---|
665 | /* store new pivot which corresponds to u[k2,k2] */ |
---|
666 | vr_piv[i] = work[qq_col[k2]]; |
---|
667 | /* new elements of i-th row of matrix V (which are non-diagonal |
---|
668 | elements u[k2,k2+1], ..., u[k2,m] of matrix U = P*V*Q) now are |
---|
669 | contained in the working array; add them to matrix V */ |
---|
670 | len = 0; |
---|
671 | for (k = k2+1; k <= m; k++) |
---|
672 | { /* get column index and value of v[i,j] = u[k2,k] */ |
---|
673 | j = qq_col[k]; |
---|
674 | temp = work[j]; |
---|
675 | /* if v[i,j] is close to zero, skip it */ |
---|
676 | if (fabs(temp) < eps_tol) continue; |
---|
677 | /* at least one unused location is needed in j-th column */ |
---|
678 | if (vc_len[j] + 1 > vc_cap[j]) |
---|
679 | { if (luf_enlarge_col(luf, j, vc_len[j] + 10)) |
---|
680 | { /* overflow of the sparse vector area */ |
---|
681 | fhv->valid = 0; |
---|
682 | luf->new_sva = luf->sv_size + luf->sv_size; |
---|
683 | xassert(luf->new_sva > luf->sv_size); |
---|
684 | ret = FHV_EROOM; |
---|
685 | goto done; |
---|
686 | } |
---|
687 | } |
---|
688 | /* add v[i,j] to j-th column */ |
---|
689 | j_ptr = vc_ptr[j] + vc_len[j]; |
---|
690 | sv_ind[j_ptr] = i; |
---|
691 | sv_val[j_ptr] = temp; |
---|
692 | vc_len[j]++; |
---|
693 | /* also store v[i,j] to the auxiliary array */ |
---|
694 | len++, cc_ind[len] = j, cc_val[len] = temp; |
---|
695 | } |
---|
696 | /* capacity of i-th row (which is currently empty) should be not |
---|
697 | less than len locations */ |
---|
698 | if (vr_cap[i] < len) |
---|
699 | { if (luf_enlarge_row(luf, i, len)) |
---|
700 | { /* overflow of the sparse vector area */ |
---|
701 | fhv->valid = 0; |
---|
702 | luf->new_sva = luf->sv_size + luf->sv_size; |
---|
703 | xassert(luf->new_sva > luf->sv_size); |
---|
704 | ret = FHV_EROOM; |
---|
705 | goto done; |
---|
706 | } |
---|
707 | } |
---|
708 | /* add new elements to i-th row list */ |
---|
709 | i_ptr = vr_ptr[i]; |
---|
710 | memmove(&sv_ind[i_ptr], &cc_ind[1], len * sizeof(int)); |
---|
711 | memmove(&sv_val[i_ptr], &cc_val[1], len * sizeof(double)); |
---|
712 | vr_len[i] = len; |
---|
713 | luf->nnz_v += len; |
---|
714 | /* updating is finished; check that diagonal element u[k2,k2] is |
---|
715 | not very small in absolute value among other elements in k2-th |
---|
716 | row and k2-th column of matrix U = P*V*Q */ |
---|
717 | /* temp = max(|u[k2,*]|, |u[*,k2]|) */ |
---|
718 | temp = 0.0; |
---|
719 | /* walk through k2-th row of U which is i-th row of V */ |
---|
720 | i = pp_row[k2]; |
---|
721 | i_beg = vr_ptr[i]; |
---|
722 | i_end = i_beg + vr_len[i] - 1; |
---|
723 | for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++) |
---|
724 | if (temp < fabs(sv_val[i_ptr])) temp = fabs(sv_val[i_ptr]); |
---|
725 | /* walk through k2-th column of U which is j-th column of V */ |
---|
726 | j = qq_col[k2]; |
---|
727 | j_beg = vc_ptr[j]; |
---|
728 | j_end = j_beg + vc_len[j] - 1; |
---|
729 | for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++) |
---|
730 | if (temp < fabs(sv_val[j_ptr])) temp = fabs(sv_val[j_ptr]); |
---|
731 | /* check that u[k2,k2] is not very small */ |
---|
732 | if (fabs(vr_piv[i]) < upd_tol * temp) |
---|
733 | { /* the factorization seems to be inaccurate and therefore must |
---|
734 | be recomputed */ |
---|
735 | fhv->valid = 0; |
---|
736 | ret = FHV_ECHECK; |
---|
737 | goto done; |
---|
738 | } |
---|
739 | /* the factorization has been successfully updated */ |
---|
740 | ret = 0; |
---|
741 | done: /* return to the calling program */ |
---|
742 | return ret; |
---|
743 | } |
---|
744 | |
---|
745 | /*********************************************************************** |
---|
746 | * NAME |
---|
747 | * |
---|
748 | * fhv_delete_it - delete LP basis factorization |
---|
749 | * |
---|
750 | * SYNOPSIS |
---|
751 | * |
---|
752 | * #include "glpfhv.h" |
---|
753 | * void fhv_delete_it(FHV *fhv); |
---|
754 | * |
---|
755 | * DESCRIPTION |
---|
756 | * |
---|
757 | * The routine fhv_delete_it deletes LP basis factorization specified |
---|
758 | * by the parameter fhv and frees all memory allocated to this program |
---|
759 | * object. */ |
---|
760 | |
---|
761 | void fhv_delete_it(FHV *fhv) |
---|
762 | { luf_delete_it(fhv->luf); |
---|
763 | if (fhv->hh_ind != NULL) xfree(fhv->hh_ind); |
---|
764 | if (fhv->hh_ptr != NULL) xfree(fhv->hh_ptr); |
---|
765 | if (fhv->hh_len != NULL) xfree(fhv->hh_len); |
---|
766 | if (fhv->p0_row != NULL) xfree(fhv->p0_row); |
---|
767 | if (fhv->p0_col != NULL) xfree(fhv->p0_col); |
---|
768 | if (fhv->cc_ind != NULL) xfree(fhv->cc_ind); |
---|
769 | if (fhv->cc_val != NULL) xfree(fhv->cc_val); |
---|
770 | xfree(fhv); |
---|
771 | return; |
---|
772 | } |
---|
773 | |
---|
774 | /* eof */ |
---|