lemon-project-template-glpk
view deps/glpk/src/glpfhv.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line source
1 /* glpfhv.c (LP basis factorization, FHV eta file version) */
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, 2011 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 ***********************************************************************/
25 #include "glpfhv.h"
26 #include "glpenv.h"
27 #define xfault xerror
29 /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
31 #define M_MAX 100000000 /* = 100*10^6 */
32 /* maximal order of the basis matrix */
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
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. */
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 }
774 /* eof */