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