lemon-project-template-glpk
view deps/glpk/src/glpluf.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 /* glpluf.c (LU-factorization) */
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 "glpenv.h"
26 #include "glpluf.h"
27 #define xfault xerror
29 /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */
31 #define N_MAX 100000000 /* = 100*10^6 */
32 /* maximal order of the original matrix */
34 /***********************************************************************
35 * NAME
36 *
37 * luf_create_it - create LU-factorization
38 *
39 * SYNOPSIS
40 *
41 * #include "glpluf.h"
42 * LUF *luf_create_it(void);
43 *
44 * DESCRIPTION
45 *
46 * The routine luf_create_it creates a program object, which represents
47 * LU-factorization of a square matrix.
48 *
49 * RETURNS
50 *
51 * The routine luf_create_it returns a pointer to the object created. */
53 LUF *luf_create_it(void)
54 { LUF *luf;
55 luf = xmalloc(sizeof(LUF));
56 luf->n_max = luf->n = 0;
57 luf->valid = 0;
58 luf->fr_ptr = luf->fr_len = NULL;
59 luf->fc_ptr = luf->fc_len = NULL;
60 luf->vr_ptr = luf->vr_len = luf->vr_cap = NULL;
61 luf->vr_piv = NULL;
62 luf->vc_ptr = luf->vc_len = luf->vc_cap = NULL;
63 luf->pp_row = luf->pp_col = NULL;
64 luf->qq_row = luf->qq_col = NULL;
65 luf->sv_size = 0;
66 luf->sv_beg = luf->sv_end = 0;
67 luf->sv_ind = NULL;
68 luf->sv_val = NULL;
69 luf->sv_head = luf->sv_tail = 0;
70 luf->sv_prev = luf->sv_next = NULL;
71 luf->vr_max = NULL;
72 luf->rs_head = luf->rs_prev = luf->rs_next = NULL;
73 luf->cs_head = luf->cs_prev = luf->cs_next = NULL;
74 luf->flag = NULL;
75 luf->work = NULL;
76 luf->new_sva = 0;
77 luf->piv_tol = 0.10;
78 luf->piv_lim = 4;
79 luf->suhl = 1;
80 luf->eps_tol = 1e-15;
81 luf->max_gro = 1e+10;
82 luf->nnz_a = luf->nnz_f = luf->nnz_v = 0;
83 luf->max_a = luf->big_v = 0.0;
84 luf->rank = 0;
85 return luf;
86 }
88 /***********************************************************************
89 * NAME
90 *
91 * luf_defrag_sva - defragment the sparse vector area
92 *
93 * SYNOPSIS
94 *
95 * #include "glpluf.h"
96 * void luf_defrag_sva(LUF *luf);
97 *
98 * DESCRIPTION
99 *
100 * The routine luf_defrag_sva defragments the sparse vector area (SVA)
101 * gathering all unused locations in one continuous extent. In order to
102 * do that the routine moves all unused locations from the left part of
103 * SVA (which contains rows and columns of the matrix V) to the middle
104 * part (which contains free locations). This is attained by relocating
105 * elements of rows and columns of the matrix V toward the beginning of
106 * the left part.
107 *
108 * NOTE that this "garbage collection" involves changing row and column
109 * pointers of the matrix V. */
111 void luf_defrag_sva(LUF *luf)
112 { int n = luf->n;
113 int *vr_ptr = luf->vr_ptr;
114 int *vr_len = luf->vr_len;
115 int *vr_cap = luf->vr_cap;
116 int *vc_ptr = luf->vc_ptr;
117 int *vc_len = luf->vc_len;
118 int *vc_cap = luf->vc_cap;
119 int *sv_ind = luf->sv_ind;
120 double *sv_val = luf->sv_val;
121 int *sv_next = luf->sv_next;
122 int sv_beg = 1;
123 int i, j, k;
124 /* skip rows and columns, which do not need to be relocated */
125 for (k = luf->sv_head; k != 0; k = sv_next[k])
126 { if (k <= n)
127 { /* i-th row of the matrix V */
128 i = k;
129 if (vr_ptr[i] != sv_beg) break;
130 vr_cap[i] = vr_len[i];
131 sv_beg += vr_cap[i];
132 }
133 else
134 { /* j-th column of the matrix V */
135 j = k - n;
136 if (vc_ptr[j] != sv_beg) break;
137 vc_cap[j] = vc_len[j];
138 sv_beg += vc_cap[j];
139 }
140 }
141 /* relocate other rows and columns in order to gather all unused
142 locations in one continuous extent */
143 for (k = k; k != 0; k = sv_next[k])
144 { if (k <= n)
145 { /* i-th row of the matrix V */
146 i = k;
147 memmove(&sv_ind[sv_beg], &sv_ind[vr_ptr[i]],
148 vr_len[i] * sizeof(int));
149 memmove(&sv_val[sv_beg], &sv_val[vr_ptr[i]],
150 vr_len[i] * sizeof(double));
151 vr_ptr[i] = sv_beg;
152 vr_cap[i] = vr_len[i];
153 sv_beg += vr_cap[i];
154 }
155 else
156 { /* j-th column of the matrix V */
157 j = k - n;
158 memmove(&sv_ind[sv_beg], &sv_ind[vc_ptr[j]],
159 vc_len[j] * sizeof(int));
160 memmove(&sv_val[sv_beg], &sv_val[vc_ptr[j]],
161 vc_len[j] * sizeof(double));
162 vc_ptr[j] = sv_beg;
163 vc_cap[j] = vc_len[j];
164 sv_beg += vc_cap[j];
165 }
166 }
167 /* set new pointer to the beginning of the free part */
168 luf->sv_beg = sv_beg;
169 return;
170 }
172 /***********************************************************************
173 * NAME
174 *
175 * luf_enlarge_row - enlarge row capacity
176 *
177 * SYNOPSIS
178 *
179 * #include "glpluf.h"
180 * int luf_enlarge_row(LUF *luf, int i, int cap);
181 *
182 * DESCRIPTION
183 *
184 * The routine luf_enlarge_row enlarges capacity of the i-th row of the
185 * matrix V to cap locations (assuming that its current capacity is less
186 * than cap). In order to do that the routine relocates elements of the
187 * i-th row to the end of the left part of SVA (which contains rows and
188 * columns of the matrix V) and then expands the left part by allocating
189 * cap free locations from the free part. If there are less than cap
190 * free locations, the routine defragments the sparse vector area.
191 *
192 * Due to "garbage collection" this operation may change row and column
193 * pointers of the matrix V.
194 *
195 * RETURNS
196 *
197 * If no error occured, the routine returns zero. Otherwise, in case of
198 * overflow of the sparse vector area, the routine returns non-zero. */
200 int luf_enlarge_row(LUF *luf, int i, int cap)
201 { int n = luf->n;
202 int *vr_ptr = luf->vr_ptr;
203 int *vr_len = luf->vr_len;
204 int *vr_cap = luf->vr_cap;
205 int *vc_cap = luf->vc_cap;
206 int *sv_ind = luf->sv_ind;
207 double *sv_val = luf->sv_val;
208 int *sv_prev = luf->sv_prev;
209 int *sv_next = luf->sv_next;
210 int ret = 0;
211 int cur, k, kk;
212 xassert(1 <= i && i <= n);
213 xassert(vr_cap[i] < cap);
214 /* if there are less than cap free locations, defragment SVA */
215 if (luf->sv_end - luf->sv_beg < cap)
216 { luf_defrag_sva(luf);
217 if (luf->sv_end - luf->sv_beg < cap)
218 { ret = 1;
219 goto done;
220 }
221 }
222 /* save current capacity of the i-th row */
223 cur = vr_cap[i];
224 /* copy existing elements to the beginning of the free part */
225 memmove(&sv_ind[luf->sv_beg], &sv_ind[vr_ptr[i]],
226 vr_len[i] * sizeof(int));
227 memmove(&sv_val[luf->sv_beg], &sv_val[vr_ptr[i]],
228 vr_len[i] * sizeof(double));
229 /* set new pointer and new capacity of the i-th row */
230 vr_ptr[i] = luf->sv_beg;
231 vr_cap[i] = cap;
232 /* set new pointer to the beginning of the free part */
233 luf->sv_beg += cap;
234 /* now the i-th row starts in the rightmost location among other
235 rows and columns of the matrix V, so its node should be moved
236 to the end of the row/column linked list */
237 k = i;
238 /* remove the i-th row node from the linked list */
239 if (sv_prev[k] == 0)
240 luf->sv_head = sv_next[k];
241 else
242 { /* capacity of the previous row/column can be increased at the
243 expense of old locations of the i-th row */
244 kk = sv_prev[k];
245 if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
246 sv_next[sv_prev[k]] = sv_next[k];
247 }
248 if (sv_next[k] == 0)
249 luf->sv_tail = sv_prev[k];
250 else
251 sv_prev[sv_next[k]] = sv_prev[k];
252 /* insert the i-th row node to the end of the linked list */
253 sv_prev[k] = luf->sv_tail;
254 sv_next[k] = 0;
255 if (sv_prev[k] == 0)
256 luf->sv_head = k;
257 else
258 sv_next[sv_prev[k]] = k;
259 luf->sv_tail = k;
260 done: return ret;
261 }
263 /***********************************************************************
264 * NAME
265 *
266 * luf_enlarge_col - enlarge column capacity
267 *
268 * SYNOPSIS
269 *
270 * #include "glpluf.h"
271 * int luf_enlarge_col(LUF *luf, int j, int cap);
272 *
273 * DESCRIPTION
274 *
275 * The routine luf_enlarge_col enlarges capacity of the j-th column of
276 * the matrix V to cap locations (assuming that its current capacity is
277 * less than cap). In order to do that the routine relocates elements
278 * of the j-th column to the end of the left part of SVA (which contains
279 * rows and columns of the matrix V) and then expands the left part by
280 * allocating cap free locations from the free part. If there are less
281 * than cap free locations, the routine defragments the sparse vector
282 * area.
283 *
284 * Due to "garbage collection" this operation may change row and column
285 * pointers of the matrix V.
286 *
287 * RETURNS
288 *
289 * If no error occured, the routine returns zero. Otherwise, in case of
290 * overflow of the sparse vector area, the routine returns non-zero. */
292 int luf_enlarge_col(LUF *luf, int j, int cap)
293 { int n = luf->n;
294 int *vr_cap = luf->vr_cap;
295 int *vc_ptr = luf->vc_ptr;
296 int *vc_len = luf->vc_len;
297 int *vc_cap = luf->vc_cap;
298 int *sv_ind = luf->sv_ind;
299 double *sv_val = luf->sv_val;
300 int *sv_prev = luf->sv_prev;
301 int *sv_next = luf->sv_next;
302 int ret = 0;
303 int cur, k, kk;
304 xassert(1 <= j && j <= n);
305 xassert(vc_cap[j] < cap);
306 /* if there are less than cap free locations, defragment SVA */
307 if (luf->sv_end - luf->sv_beg < cap)
308 { luf_defrag_sva(luf);
309 if (luf->sv_end - luf->sv_beg < cap)
310 { ret = 1;
311 goto done;
312 }
313 }
314 /* save current capacity of the j-th column */
315 cur = vc_cap[j];
316 /* copy existing elements to the beginning of the free part */
317 memmove(&sv_ind[luf->sv_beg], &sv_ind[vc_ptr[j]],
318 vc_len[j] * sizeof(int));
319 memmove(&sv_val[luf->sv_beg], &sv_val[vc_ptr[j]],
320 vc_len[j] * sizeof(double));
321 /* set new pointer and new capacity of the j-th column */
322 vc_ptr[j] = luf->sv_beg;
323 vc_cap[j] = cap;
324 /* set new pointer to the beginning of the free part */
325 luf->sv_beg += cap;
326 /* now the j-th column starts in the rightmost location among
327 other rows and columns of the matrix V, so its node should be
328 moved to the end of the row/column linked list */
329 k = n + j;
330 /* remove the j-th column node from the linked list */
331 if (sv_prev[k] == 0)
332 luf->sv_head = sv_next[k];
333 else
334 { /* capacity of the previous row/column can be increased at the
335 expense of old locations of the j-th column */
336 kk = sv_prev[k];
337 if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
338 sv_next[sv_prev[k]] = sv_next[k];
339 }
340 if (sv_next[k] == 0)
341 luf->sv_tail = sv_prev[k];
342 else
343 sv_prev[sv_next[k]] = sv_prev[k];
344 /* insert the j-th column node to the end of the linked list */
345 sv_prev[k] = luf->sv_tail;
346 sv_next[k] = 0;
347 if (sv_prev[k] == 0)
348 luf->sv_head = k;
349 else
350 sv_next[sv_prev[k]] = k;
351 luf->sv_tail = k;
352 done: return ret;
353 }
355 /***********************************************************************
356 * reallocate - reallocate LU-factorization arrays
357 *
358 * This routine reallocates arrays, whose size depends of n, the order
359 * of the matrix A to be factorized. */
361 static void reallocate(LUF *luf, int n)
362 { int n_max = luf->n_max;
363 luf->n = n;
364 if (n <= n_max) goto done;
365 if (luf->fr_ptr != NULL) xfree(luf->fr_ptr);
366 if (luf->fr_len != NULL) xfree(luf->fr_len);
367 if (luf->fc_ptr != NULL) xfree(luf->fc_ptr);
368 if (luf->fc_len != NULL) xfree(luf->fc_len);
369 if (luf->vr_ptr != NULL) xfree(luf->vr_ptr);
370 if (luf->vr_len != NULL) xfree(luf->vr_len);
371 if (luf->vr_cap != NULL) xfree(luf->vr_cap);
372 if (luf->vr_piv != NULL) xfree(luf->vr_piv);
373 if (luf->vc_ptr != NULL) xfree(luf->vc_ptr);
374 if (luf->vc_len != NULL) xfree(luf->vc_len);
375 if (luf->vc_cap != NULL) xfree(luf->vc_cap);
376 if (luf->pp_row != NULL) xfree(luf->pp_row);
377 if (luf->pp_col != NULL) xfree(luf->pp_col);
378 if (luf->qq_row != NULL) xfree(luf->qq_row);
379 if (luf->qq_col != NULL) xfree(luf->qq_col);
380 if (luf->sv_prev != NULL) xfree(luf->sv_prev);
381 if (luf->sv_next != NULL) xfree(luf->sv_next);
382 if (luf->vr_max != NULL) xfree(luf->vr_max);
383 if (luf->rs_head != NULL) xfree(luf->rs_head);
384 if (luf->rs_prev != NULL) xfree(luf->rs_prev);
385 if (luf->rs_next != NULL) xfree(luf->rs_next);
386 if (luf->cs_head != NULL) xfree(luf->cs_head);
387 if (luf->cs_prev != NULL) xfree(luf->cs_prev);
388 if (luf->cs_next != NULL) xfree(luf->cs_next);
389 if (luf->flag != NULL) xfree(luf->flag);
390 if (luf->work != NULL) xfree(luf->work);
391 luf->n_max = n_max = n + 100;
392 luf->fr_ptr = xcalloc(1+n_max, sizeof(int));
393 luf->fr_len = xcalloc(1+n_max, sizeof(int));
394 luf->fc_ptr = xcalloc(1+n_max, sizeof(int));
395 luf->fc_len = xcalloc(1+n_max, sizeof(int));
396 luf->vr_ptr = xcalloc(1+n_max, sizeof(int));
397 luf->vr_len = xcalloc(1+n_max, sizeof(int));
398 luf->vr_cap = xcalloc(1+n_max, sizeof(int));
399 luf->vr_piv = xcalloc(1+n_max, sizeof(double));
400 luf->vc_ptr = xcalloc(1+n_max, sizeof(int));
401 luf->vc_len = xcalloc(1+n_max, sizeof(int));
402 luf->vc_cap = xcalloc(1+n_max, sizeof(int));
403 luf->pp_row = xcalloc(1+n_max, sizeof(int));
404 luf->pp_col = xcalloc(1+n_max, sizeof(int));
405 luf->qq_row = xcalloc(1+n_max, sizeof(int));
406 luf->qq_col = xcalloc(1+n_max, sizeof(int));
407 luf->sv_prev = xcalloc(1+n_max+n_max, sizeof(int));
408 luf->sv_next = xcalloc(1+n_max+n_max, sizeof(int));
409 luf->vr_max = xcalloc(1+n_max, sizeof(double));
410 luf->rs_head = xcalloc(1+n_max, sizeof(int));
411 luf->rs_prev = xcalloc(1+n_max, sizeof(int));
412 luf->rs_next = xcalloc(1+n_max, sizeof(int));
413 luf->cs_head = xcalloc(1+n_max, sizeof(int));
414 luf->cs_prev = xcalloc(1+n_max, sizeof(int));
415 luf->cs_next = xcalloc(1+n_max, sizeof(int));
416 luf->flag = xcalloc(1+n_max, sizeof(int));
417 luf->work = xcalloc(1+n_max, sizeof(double));
418 done: return;
419 }
421 /***********************************************************************
422 * initialize - initialize LU-factorization data structures
423 *
424 * This routine initializes data structures for subsequent computing
425 * the LU-factorization of a given matrix A, which is specified by the
426 * formal routine col. On exit V = A and F = P = Q = I, where I is the
427 * unity matrix. (Row-wise representation of the matrix F is not used
428 * at the factorization stage and therefore is not initialized.)
429 *
430 * If no error occured, the routine returns zero. Otherwise, in case of
431 * overflow of the sparse vector area, the routine returns non-zero. */
433 static int initialize(LUF *luf, int (*col)(void *info, int j, int rn[],
434 double aj[]), void *info)
435 { int n = luf->n;
436 int *fc_ptr = luf->fc_ptr;
437 int *fc_len = luf->fc_len;
438 int *vr_ptr = luf->vr_ptr;
439 int *vr_len = luf->vr_len;
440 int *vr_cap = luf->vr_cap;
441 int *vc_ptr = luf->vc_ptr;
442 int *vc_len = luf->vc_len;
443 int *vc_cap = luf->vc_cap;
444 int *pp_row = luf->pp_row;
445 int *pp_col = luf->pp_col;
446 int *qq_row = luf->qq_row;
447 int *qq_col = luf->qq_col;
448 int *sv_ind = luf->sv_ind;
449 double *sv_val = luf->sv_val;
450 int *sv_prev = luf->sv_prev;
451 int *sv_next = luf->sv_next;
452 double *vr_max = luf->vr_max;
453 int *rs_head = luf->rs_head;
454 int *rs_prev = luf->rs_prev;
455 int *rs_next = luf->rs_next;
456 int *cs_head = luf->cs_head;
457 int *cs_prev = luf->cs_prev;
458 int *cs_next = luf->cs_next;
459 int *flag = luf->flag;
460 double *work = luf->work;
461 int ret = 0;
462 int i, i_ptr, j, j_beg, j_end, k, len, nnz, sv_beg, sv_end, ptr;
463 double big, val;
464 /* free all locations of the sparse vector area */
465 sv_beg = 1;
466 sv_end = luf->sv_size + 1;
467 /* (row-wise representation of the matrix F is not initialized,
468 because it is not used at the factorization stage) */
469 /* build the matrix F in column-wise format (initially F = I) */
470 for (j = 1; j <= n; j++)
471 { fc_ptr[j] = sv_end;
472 fc_len[j] = 0;
473 }
474 /* clear rows of the matrix V; clear the flag array */
475 for (i = 1; i <= n; i++)
476 vr_len[i] = vr_cap[i] = 0, flag[i] = 0;
477 /* build the matrix V in column-wise format (initially V = A);
478 count non-zeros in rows of this matrix; count total number of
479 non-zeros; compute largest of absolute values of elements */
480 nnz = 0;
481 big = 0.0;
482 for (j = 1; j <= n; j++)
483 { int *rn = pp_row;
484 double *aj = work;
485 /* obtain j-th column of the matrix A */
486 len = col(info, j, rn, aj);
487 if (!(0 <= len && len <= n))
488 xfault("luf_factorize: j = %d; len = %d; invalid column len"
489 "gth\n", j, len);
490 /* check for free locations */
491 if (sv_end - sv_beg < len)
492 { /* overflow of the sparse vector area */
493 ret = 1;
494 goto done;
495 }
496 /* set pointer to the j-th column */
497 vc_ptr[j] = sv_beg;
498 /* set length of the j-th column */
499 vc_len[j] = vc_cap[j] = len;
500 /* count total number of non-zeros */
501 nnz += len;
502 /* walk through elements of the j-th column */
503 for (ptr = 1; ptr <= len; ptr++)
504 { /* get row index and numerical value of a[i,j] */
505 i = rn[ptr];
506 val = aj[ptr];
507 if (!(1 <= i && i <= n))
508 xfault("luf_factorize: i = %d; j = %d; invalid row index"
509 "\n", i, j);
510 if (flag[i])
511 xfault("luf_factorize: i = %d; j = %d; duplicate element"
512 " not allowed\n", i, j);
513 if (val == 0.0)
514 xfault("luf_factorize: i = %d; j = %d; zero element not "
515 "allowed\n", i, j);
516 /* add new element v[i,j] = a[i,j] to j-th column */
517 sv_ind[sv_beg] = i;
518 sv_val[sv_beg] = val;
519 sv_beg++;
520 /* big := max(big, |a[i,j]|) */
521 if (val < 0.0) val = - val;
522 if (big < val) big = val;
523 /* mark non-zero in the i-th position of the j-th column */
524 flag[i] = 1;
525 /* increase length of the i-th row */
526 vr_cap[i]++;
527 }
528 /* reset all non-zero marks */
529 for (ptr = 1; ptr <= len; ptr++) flag[rn[ptr]] = 0;
530 }
531 /* allocate rows of the matrix V */
532 for (i = 1; i <= n; i++)
533 { /* get length of the i-th row */
534 len = vr_cap[i];
535 /* check for free locations */
536 if (sv_end - sv_beg < len)
537 { /* overflow of the sparse vector area */
538 ret = 1;
539 goto done;
540 }
541 /* set pointer to the i-th row */
542 vr_ptr[i] = sv_beg;
543 /* reserve locations for the i-th row */
544 sv_beg += len;
545 }
546 /* build the matrix V in row-wise format using representation of
547 this matrix in column-wise format */
548 for (j = 1; j <= n; j++)
549 { /* walk through elements of the j-th column */
550 j_beg = vc_ptr[j];
551 j_end = j_beg + vc_len[j] - 1;
552 for (k = j_beg; k <= j_end; k++)
553 { /* get row index and numerical value of v[i,j] */
554 i = sv_ind[k];
555 val = sv_val[k];
556 /* store element in the i-th row */
557 i_ptr = vr_ptr[i] + vr_len[i];
558 sv_ind[i_ptr] = j;
559 sv_val[i_ptr] = val;
560 /* increase count of the i-th row */
561 vr_len[i]++;
562 }
563 }
564 /* initialize the matrices P and Q (initially P = Q = I) */
565 for (k = 1; k <= n; k++)
566 pp_row[k] = pp_col[k] = qq_row[k] = qq_col[k] = k;
567 /* set sva partitioning pointers */
568 luf->sv_beg = sv_beg;
569 luf->sv_end = sv_end;
570 /* the initial physical order of rows and columns of the matrix V
571 is n+1, ..., n+n, 1, ..., n (firstly columns, then rows) */
572 luf->sv_head = n+1;
573 luf->sv_tail = n;
574 for (i = 1; i <= n; i++)
575 { sv_prev[i] = i-1;
576 sv_next[i] = i+1;
577 }
578 sv_prev[1] = n+n;
579 sv_next[n] = 0;
580 for (j = 1; j <= n; j++)
581 { sv_prev[n+j] = n+j-1;
582 sv_next[n+j] = n+j+1;
583 }
584 sv_prev[n+1] = 0;
585 sv_next[n+n] = 1;
586 /* clear working arrays */
587 for (k = 1; k <= n; k++)
588 { flag[k] = 0;
589 work[k] = 0.0;
590 }
591 /* initialize some statistics */
592 luf->nnz_a = nnz;
593 luf->nnz_f = 0;
594 luf->nnz_v = nnz;
595 luf->max_a = big;
596 luf->big_v = big;
597 luf->rank = -1;
598 /* initially the active submatrix is the entire matrix V */
599 /* largest of absolute values of elements in each active row is
600 unknown yet */
601 for (i = 1; i <= n; i++) vr_max[i] = -1.0;
602 /* build linked lists of active rows */
603 for (len = 0; len <= n; len++) rs_head[len] = 0;
604 for (i = 1; i <= n; i++)
605 { len = vr_len[i];
606 rs_prev[i] = 0;
607 rs_next[i] = rs_head[len];
608 if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
609 rs_head[len] = i;
610 }
611 /* build linked lists of active columns */
612 for (len = 0; len <= n; len++) cs_head[len] = 0;
613 for (j = 1; j <= n; j++)
614 { len = vc_len[j];
615 cs_prev[j] = 0;
616 cs_next[j] = cs_head[len];
617 if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
618 cs_head[len] = j;
619 }
620 done: /* return to the factorizing routine */
621 return ret;
622 }
624 /***********************************************************************
625 * find_pivot - choose a pivot element
626 *
627 * This routine chooses a pivot element in the active submatrix of the
628 * matrix U = P*V*Q.
629 *
630 * It is assumed that on entry the matrix U has the following partially
631 * triangularized form:
632 *
633 * 1 k n
634 * 1 x x x x x x x x x x
635 * . x x x x x x x x x
636 * . . x x x x x x x x
637 * . . . x x x x x x x
638 * k . . . . * * * * * *
639 * . . . . * * * * * *
640 * . . . . * * * * * *
641 * . . . . * * * * * *
642 * . . . . * * * * * *
643 * n . . . . * * * * * *
644 *
645 * where rows and columns k, k+1, ..., n belong to the active submatrix
646 * (elements of the active submatrix are marked by '*').
647 *
648 * Since the matrix U = P*V*Q is not stored, the routine works with the
649 * matrix V. It is assumed that the row-wise representation corresponds
650 * to the matrix V, but the column-wise representation corresponds to
651 * the active submatrix of the matrix V, i.e. elements of the matrix V,
652 * which doesn't belong to the active submatrix, are missing from the
653 * column linked lists. It is also assumed that each active row of the
654 * matrix V is in the set R[len], where len is number of non-zeros in
655 * the row, and each active column of the matrix V is in the set C[len],
656 * where len is number of non-zeros in the column (in the latter case
657 * only elements of the active submatrix are counted; such elements are
658 * marked by '*' on the figure above).
659 *
660 * For the reason of numerical stability the routine applies so called
661 * threshold pivoting proposed by J.Reid. It is assumed that an element
662 * v[i,j] can be selected as a pivot candidate if it is not very small
663 * (in absolute value) among other elements in the same row, i.e. if it
664 * satisfies to the stability condition |v[i,j]| >= tol * max|v[i,*]|,
665 * where 0 < tol < 1 is a given tolerance.
666 *
667 * In order to keep sparsity of the matrix V the routine uses Markowitz
668 * strategy, trying to choose such element v[p,q], which satisfies to
669 * the stability condition (see above) and has smallest Markowitz cost
670 * (nr[p]-1) * (nc[q]-1), where nr[p] and nc[q] are numbers of non-zero
671 * elements, respectively, in the p-th row and in the q-th column of the
672 * active submatrix.
673 *
674 * In order to reduce the search, i.e. not to walk through all elements
675 * of the active submatrix, the routine exploits a technique proposed by
676 * I.Duff. This technique is based on using the sets R[len] and C[len]
677 * of active rows and columns.
678 *
679 * If the pivot element v[p,q] has been chosen, the routine stores its
680 * indices to the locations *p and *q and returns zero. Otherwise, if
681 * the active submatrix is empty and therefore the pivot element can't
682 * be chosen, the routine returns non-zero. */
684 static int find_pivot(LUF *luf, int *_p, int *_q)
685 { int n = luf->n;
686 int *vr_ptr = luf->vr_ptr;
687 int *vr_len = luf->vr_len;
688 int *vc_ptr = luf->vc_ptr;
689 int *vc_len = luf->vc_len;
690 int *sv_ind = luf->sv_ind;
691 double *sv_val = luf->sv_val;
692 double *vr_max = luf->vr_max;
693 int *rs_head = luf->rs_head;
694 int *rs_next = luf->rs_next;
695 int *cs_head = luf->cs_head;
696 int *cs_prev = luf->cs_prev;
697 int *cs_next = luf->cs_next;
698 double piv_tol = luf->piv_tol;
699 int piv_lim = luf->piv_lim;
700 int suhl = luf->suhl;
701 int p, q, len, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr,
702 ncand, next_j, min_p, min_q, min_len;
703 double best, cost, big, temp;
704 /* initially no pivot candidates have been found so far */
705 p = q = 0, best = DBL_MAX, ncand = 0;
706 /* if in the active submatrix there is a column that has the only
707 non-zero (column singleton), choose it as pivot */
708 j = cs_head[1];
709 if (j != 0)
710 { xassert(vc_len[j] == 1);
711 p = sv_ind[vc_ptr[j]], q = j;
712 goto done;
713 }
714 /* if in the active submatrix there is a row that has the only
715 non-zero (row singleton), choose it as pivot */
716 i = rs_head[1];
717 if (i != 0)
718 { xassert(vr_len[i] == 1);
719 p = i, q = sv_ind[vr_ptr[i]];
720 goto done;
721 }
722 /* there are no singletons in the active submatrix; walk through
723 other non-empty rows and columns */
724 for (len = 2; len <= n; len++)
725 { /* consider active columns that have len non-zeros */
726 for (j = cs_head[len]; j != 0; j = next_j)
727 { /* the j-th column has len non-zeros */
728 j_beg = vc_ptr[j];
729 j_end = j_beg + vc_len[j] - 1;
730 /* save pointer to the next column with the same length */
731 next_j = cs_next[j];
732 /* find an element in the j-th column, which is placed in a
733 row with minimal number of non-zeros and satisfies to the
734 stability condition (such element may not exist) */
735 min_p = min_q = 0, min_len = INT_MAX;
736 for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
737 { /* get row index of v[i,j] */
738 i = sv_ind[j_ptr];
739 i_beg = vr_ptr[i];
740 i_end = i_beg + vr_len[i] - 1;
741 /* if the i-th row is not shorter than that one, where
742 minimal element is currently placed, skip v[i,j] */
743 if (vr_len[i] >= min_len) continue;
744 /* determine the largest of absolute values of elements
745 in the i-th row */
746 big = vr_max[i];
747 if (big < 0.0)
748 { /* the largest value is unknown yet; compute it */
749 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
750 { temp = sv_val[i_ptr];
751 if (temp < 0.0) temp = - temp;
752 if (big < temp) big = temp;
753 }
754 vr_max[i] = big;
755 }
756 /* find v[i,j] in the i-th row */
757 for (i_ptr = vr_ptr[i]; sv_ind[i_ptr] != j; i_ptr++);
758 xassert(i_ptr <= i_end);
759 /* if v[i,j] doesn't satisfy to the stability condition,
760 skip it */
761 temp = sv_val[i_ptr];
762 if (temp < 0.0) temp = - temp;
763 if (temp < piv_tol * big) continue;
764 /* v[i,j] is better than the current minimal element */
765 min_p = i, min_q = j, min_len = vr_len[i];
766 /* if Markowitz cost of the current minimal element is
767 not greater than (len-1)**2, it can be chosen right
768 now; this heuristic reduces the search and works well
769 in many cases */
770 if (min_len <= len)
771 { p = min_p, q = min_q;
772 goto done;
773 }
774 }
775 /* the j-th column has been scanned */
776 if (min_p != 0)
777 { /* the minimal element is a next pivot candidate */
778 ncand++;
779 /* compute its Markowitz cost */
780 cost = (double)(min_len - 1) * (double)(len - 1);
781 /* choose between the minimal element and the current
782 candidate */
783 if (cost < best) p = min_p, q = min_q, best = cost;
784 /* if piv_lim candidates have been considered, there are
785 doubts that a much better candidate exists; therefore
786 it's time to terminate the search */
787 if (ncand == piv_lim) goto done;
788 }
789 else
790 { /* the j-th column has no elements, which satisfy to the
791 stability condition; Uwe Suhl suggests to exclude such
792 column from the further consideration until it becomes
793 a column singleton; in hard cases this significantly
794 reduces a time needed for pivot searching */
795 if (suhl)
796 { /* remove the j-th column from the active set */
797 if (cs_prev[j] == 0)
798 cs_head[len] = cs_next[j];
799 else
800 cs_next[cs_prev[j]] = cs_next[j];
801 if (cs_next[j] == 0)
802 /* nop */;
803 else
804 cs_prev[cs_next[j]] = cs_prev[j];
805 /* the following assignment is used to avoid an error
806 when the routine eliminate (see below) will try to
807 remove the j-th column from the active set */
808 cs_prev[j] = cs_next[j] = j;
809 }
810 }
811 }
812 /* consider active rows that have len non-zeros */
813 for (i = rs_head[len]; i != 0; i = rs_next[i])
814 { /* the i-th row has len non-zeros */
815 i_beg = vr_ptr[i];
816 i_end = i_beg + vr_len[i] - 1;
817 /* determine the largest of absolute values of elements in
818 the i-th row */
819 big = vr_max[i];
820 if (big < 0.0)
821 { /* the largest value is unknown yet; compute it */
822 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
823 { temp = sv_val[i_ptr];
824 if (temp < 0.0) temp = - temp;
825 if (big < temp) big = temp;
826 }
827 vr_max[i] = big;
828 }
829 /* find an element in the i-th row, which is placed in a
830 column with minimal number of non-zeros and satisfies to
831 the stability condition (such element always exists) */
832 min_p = min_q = 0, min_len = INT_MAX;
833 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
834 { /* get column index of v[i,j] */
835 j = sv_ind[i_ptr];
836 /* if the j-th column is not shorter than that one, where
837 minimal element is currently placed, skip v[i,j] */
838 if (vc_len[j] >= min_len) continue;
839 /* if v[i,j] doesn't satisfy to the stability condition,
840 skip it */
841 temp = sv_val[i_ptr];
842 if (temp < 0.0) temp = - temp;
843 if (temp < piv_tol * big) continue;
844 /* v[i,j] is better than the current minimal element */
845 min_p = i, min_q = j, min_len = vc_len[j];
846 /* if Markowitz cost of the current minimal element is
847 not greater than (len-1)**2, it can be chosen right
848 now; this heuristic reduces the search and works well
849 in many cases */
850 if (min_len <= len)
851 { p = min_p, q = min_q;
852 goto done;
853 }
854 }
855 /* the i-th row has been scanned */
856 if (min_p != 0)
857 { /* the minimal element is a next pivot candidate */
858 ncand++;
859 /* compute its Markowitz cost */
860 cost = (double)(len - 1) * (double)(min_len - 1);
861 /* choose between the minimal element and the current
862 candidate */
863 if (cost < best) p = min_p, q = min_q, best = cost;
864 /* if piv_lim candidates have been considered, there are
865 doubts that a much better candidate exists; therefore
866 it's time to terminate the search */
867 if (ncand == piv_lim) goto done;
868 }
869 else
870 { /* this can't be because this can never be */
871 xassert(min_p != min_p);
872 }
873 }
874 }
875 done: /* bring the pivot to the factorizing routine */
876 *_p = p, *_q = q;
877 return (p == 0);
878 }
880 /***********************************************************************
881 * eliminate - perform gaussian elimination.
882 *
883 * This routine performs elementary gaussian transformations in order
884 * to eliminate subdiagonal elements in the k-th column of the matrix
885 * U = P*V*Q using the pivot element u[k,k], where k is the number of
886 * the current elimination step.
887 *
888 * The parameters p and q are, respectively, row and column indices of
889 * the element v[p,q], which corresponds to the element u[k,k].
890 *
891 * Each time when the routine applies the elementary transformation to
892 * a non-pivot row of the matrix V, it stores the corresponding element
893 * to the matrix F in order to keep the main equality A = F*V.
894 *
895 * The routine assumes that on entry the matrices L = P*F*inv(P) and
896 * U = P*V*Q are the following:
897 *
898 * 1 k 1 k n
899 * 1 1 . . . . . . . . . 1 x x x x x x x x x x
900 * x 1 . . . . . . . . . x x x x x x x x x
901 * x x 1 . . . . . . . . . x x x x x x x x
902 * x x x 1 . . . . . . . . . x x x x x x x
903 * k x x x x 1 . . . . . k . . . . * * * * * *
904 * x x x x _ 1 . . . . . . . . # * * * * *
905 * x x x x _ . 1 . . . . . . . # * * * * *
906 * x x x x _ . . 1 . . . . . . # * * * * *
907 * x x x x _ . . . 1 . . . . . # * * * * *
908 * n x x x x _ . . . . 1 n . . . . # * * * * *
909 *
910 * matrix L matrix U
911 *
912 * where rows and columns of the matrix U with numbers k, k+1, ..., n
913 * form the active submatrix (eliminated elements are marked by '#' and
914 * other elements of the active submatrix are marked by '*'). Note that
915 * each eliminated non-zero element u[i,k] of the matrix U gives the
916 * corresponding element l[i,k] of the matrix L (marked by '_').
917 *
918 * Actually all operations are performed on the matrix V. Should note
919 * that the row-wise representation corresponds to the matrix V, but the
920 * column-wise representation corresponds to the active submatrix of the
921 * matrix V, i.e. elements of the matrix V, which doesn't belong to the
922 * active submatrix, are missing from the column linked lists.
923 *
924 * Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal
925 * elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies
926 * the following elementary gaussian transformations:
927 *
928 * (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V),
929 *
930 * where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier.
931 *
932 * Additionally, in order to keep the main equality A = F*V, each time
933 * when the routine applies the transformation to i-th row of the matrix
934 * V, it also adds f[i,p] as a new element to the matrix F.
935 *
936 * IMPORTANT: On entry the working arrays flag and work should contain
937 * zeros. This status is provided by the routine on exit.
938 *
939 * If no error occured, the routine returns zero. Otherwise, in case of
940 * overflow of the sparse vector area, the routine returns non-zero. */
942 static int eliminate(LUF *luf, int p, int q)
943 { int n = luf->n;
944 int *fc_ptr = luf->fc_ptr;
945 int *fc_len = luf->fc_len;
946 int *vr_ptr = luf->vr_ptr;
947 int *vr_len = luf->vr_len;
948 int *vr_cap = luf->vr_cap;
949 double *vr_piv = luf->vr_piv;
950 int *vc_ptr = luf->vc_ptr;
951 int *vc_len = luf->vc_len;
952 int *vc_cap = luf->vc_cap;
953 int *sv_ind = luf->sv_ind;
954 double *sv_val = luf->sv_val;
955 int *sv_prev = luf->sv_prev;
956 int *sv_next = luf->sv_next;
957 double *vr_max = luf->vr_max;
958 int *rs_head = luf->rs_head;
959 int *rs_prev = luf->rs_prev;
960 int *rs_next = luf->rs_next;
961 int *cs_head = luf->cs_head;
962 int *cs_prev = luf->cs_prev;
963 int *cs_next = luf->cs_next;
964 int *flag = luf->flag;
965 double *work = luf->work;
966 double eps_tol = luf->eps_tol;
967 /* at this stage the row-wise representation of the matrix F is
968 not used, so fr_len can be used as a working array */
969 int *ndx = luf->fr_len;
970 int ret = 0;
971 int len, fill, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, k,
972 p_beg, p_end, p_ptr, q_beg, q_end, q_ptr;
973 double fip, val, vpq, temp;
974 xassert(1 <= p && p <= n);
975 xassert(1 <= q && q <= n);
976 /* remove the p-th (pivot) row from the active set; this row will
977 never return there */
978 if (rs_prev[p] == 0)
979 rs_head[vr_len[p]] = rs_next[p];
980 else
981 rs_next[rs_prev[p]] = rs_next[p];
982 if (rs_next[p] == 0)
983 ;
984 else
985 rs_prev[rs_next[p]] = rs_prev[p];
986 /* remove the q-th (pivot) column from the active set; this column
987 will never return there */
988 if (cs_prev[q] == 0)
989 cs_head[vc_len[q]] = cs_next[q];
990 else
991 cs_next[cs_prev[q]] = cs_next[q];
992 if (cs_next[q] == 0)
993 ;
994 else
995 cs_prev[cs_next[q]] = cs_prev[q];
996 /* find the pivot v[p,q] = u[k,k] in the p-th row */
997 p_beg = vr_ptr[p];
998 p_end = p_beg + vr_len[p] - 1;
999 for (p_ptr = p_beg; sv_ind[p_ptr] != q; p_ptr++) /* nop */;
1000 xassert(p_ptr <= p_end);
1001 /* store value of the pivot */
1002 vpq = (vr_piv[p] = sv_val[p_ptr]);
1003 /* remove the pivot from the p-th row */
1004 sv_ind[p_ptr] = sv_ind[p_end];
1005 sv_val[p_ptr] = sv_val[p_end];
1006 vr_len[p]--;
1007 p_end--;
1008 /* find the pivot v[p,q] = u[k,k] in the q-th column */
1009 q_beg = vc_ptr[q];
1010 q_end = q_beg + vc_len[q] - 1;
1011 for (q_ptr = q_beg; sv_ind[q_ptr] != p; q_ptr++) /* nop */;
1012 xassert(q_ptr <= q_end);
1013 /* remove the pivot from the q-th column */
1014 sv_ind[q_ptr] = sv_ind[q_end];
1015 vc_len[q]--;
1016 q_end--;
1017 /* walk through the p-th (pivot) row, which doesn't contain the
1018 pivot v[p,q] already, and do the following... */
1019 for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
1020 { /* get column index of v[p,j] */
1021 j = sv_ind[p_ptr];
1022 /* store v[p,j] to the working array */
1023 flag[j] = 1;
1024 work[j] = sv_val[p_ptr];
1025 /* remove the j-th column from the active set; this column will
1026 return there later with new length */
1027 if (cs_prev[j] == 0)
1028 cs_head[vc_len[j]] = cs_next[j];
1029 else
1030 cs_next[cs_prev[j]] = cs_next[j];
1031 if (cs_next[j] == 0)
1032 ;
1033 else
1034 cs_prev[cs_next[j]] = cs_prev[j];
1035 /* find v[p,j] in the j-th column */
1036 j_beg = vc_ptr[j];
1037 j_end = j_beg + vc_len[j] - 1;
1038 for (j_ptr = j_beg; sv_ind[j_ptr] != p; j_ptr++) /* nop */;
1039 xassert(j_ptr <= j_end);
1040 /* since v[p,j] leaves the active submatrix, remove it from the
1041 j-th column; however, v[p,j] is kept in the p-th row */
1042 sv_ind[j_ptr] = sv_ind[j_end];
1043 vc_len[j]--;
1044 }
1045 /* walk through the q-th (pivot) column, which doesn't contain the
1046 pivot v[p,q] already, and perform gaussian elimination */
1047 while (q_beg <= q_end)
1048 { /* element v[i,q] should be eliminated */
1049 /* get row index of v[i,q] */
1050 i = sv_ind[q_beg];
1051 /* remove the i-th row from the active set; later this row will
1052 return there with new length */
1053 if (rs_prev[i] == 0)
1054 rs_head[vr_len[i]] = rs_next[i];
1055 else
1056 rs_next[rs_prev[i]] = rs_next[i];
1057 if (rs_next[i] == 0)
1058 ;
1059 else
1060 rs_prev[rs_next[i]] = rs_prev[i];
1061 /* find v[i,q] in the i-th row */
1062 i_beg = vr_ptr[i];
1063 i_end = i_beg + vr_len[i] - 1;
1064 for (i_ptr = i_beg; sv_ind[i_ptr] != q; i_ptr++) /* nop */;
1065 xassert(i_ptr <= i_end);
1066 /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] */
1067 fip = sv_val[i_ptr] / vpq;
1068 /* since v[i,q] should be eliminated, remove it from the i-th
1069 row */
1070 sv_ind[i_ptr] = sv_ind[i_end];
1071 sv_val[i_ptr] = sv_val[i_end];
1072 vr_len[i]--;
1073 i_end--;
1074 /* and from the q-th column */
1075 sv_ind[q_beg] = sv_ind[q_end];
1076 vc_len[q]--;
1077 q_end--;
1078 /* perform gaussian transformation:
1079 (i-th row) := (i-th row) - f[i,p] * (p-th row)
1080 note that now the p-th row, which is in the working array,
1081 doesn't contain the pivot v[p,q], and the i-th row doesn't
1082 contain the eliminated element v[i,q] */
1083 /* walk through the i-th row and transform existing non-zero
1084 elements */
1085 fill = vr_len[p];
1086 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
1087 { /* get column index of v[i,j] */
1088 j = sv_ind[i_ptr];
1089 /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
1090 if (flag[j])
1091 { /* v[p,j] != 0 */
1092 temp = (sv_val[i_ptr] -= fip * work[j]);
1093 if (temp < 0.0) temp = - temp;
1094 flag[j] = 0;
1095 fill--; /* since both v[i,j] and v[p,j] exist */
1096 if (temp == 0.0 || temp < eps_tol)
1097 { /* new v[i,j] is closer to zero; replace it by exact
1098 zero, i.e. remove it from the active submatrix */
1099 /* remove v[i,j] from the i-th row */
1100 sv_ind[i_ptr] = sv_ind[i_end];
1101 sv_val[i_ptr] = sv_val[i_end];
1102 vr_len[i]--;
1103 i_ptr--;
1104 i_end--;
1105 /* find v[i,j] in the j-th column */
1106 j_beg = vc_ptr[j];
1107 j_end = j_beg + vc_len[j] - 1;
1108 for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++);
1109 xassert(j_ptr <= j_end);
1110 /* remove v[i,j] from the j-th column */
1111 sv_ind[j_ptr] = sv_ind[j_end];
1112 vc_len[j]--;
1113 }
1114 else
1115 { /* v_big := max(v_big, |v[i,j]|) */
1116 if (luf->big_v < temp) luf->big_v = temp;
1117 }
1118 }
1119 }
1120 /* now flag is the pattern of the set v[p,*] \ v[i,*], and fill
1121 is number of non-zeros in this set; therefore up to fill new
1122 non-zeros may appear in the i-th row */
1123 if (vr_len[i] + fill > vr_cap[i])
1124 { /* enlarge the i-th row */
1125 if (luf_enlarge_row(luf, i, vr_len[i] + fill))
1126 { /* overflow of the sparse vector area */
1127 ret = 1;
1128 goto done;
1129 }
1130 /* defragmentation may change row and column pointers of the
1131 matrix V */
1132 p_beg = vr_ptr[p];
1133 p_end = p_beg + vr_len[p] - 1;
1134 q_beg = vc_ptr[q];
1135 q_end = q_beg + vc_len[q] - 1;
1136 }
1137 /* walk through the p-th (pivot) row and create new elements
1138 of the i-th row that appear due to fill-in; column indices
1139 of these new elements are accumulated in the array ndx */
1140 len = 0;
1141 for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
1142 { /* get column index of v[p,j], which may cause fill-in */
1143 j = sv_ind[p_ptr];
1144 if (flag[j])
1145 { /* compute new non-zero v[i,j] = 0 - f[i,p] * v[p,j] */
1146 temp = (val = - fip * work[j]);
1147 if (temp < 0.0) temp = - temp;
1148 if (temp == 0.0 || temp < eps_tol)
1149 /* if v[i,j] is closer to zero; just ignore it */;
1150 else
1151 { /* add v[i,j] to the i-th row */
1152 i_ptr = vr_ptr[i] + vr_len[i];
1153 sv_ind[i_ptr] = j;
1154 sv_val[i_ptr] = val;
1155 vr_len[i]++;
1156 /* remember column index of v[i,j] */
1157 ndx[++len] = j;
1158 /* big_v := max(big_v, |v[i,j]|) */
1159 if (luf->big_v < temp) luf->big_v = temp;
1160 }
1161 }
1162 else
1163 { /* there is no fill-in, because v[i,j] already exists in
1164 the i-th row; restore the flag of the element v[p,j],
1165 which was reset before */
1166 flag[j] = 1;
1167 }
1168 }
1169 /* add new non-zeros v[i,j] to the corresponding columns */
1170 for (k = 1; k <= len; k++)
1171 { /* get column index of new non-zero v[i,j] */
1172 j = ndx[k];
1173 /* one free location is needed in the j-th column */
1174 if (vc_len[j] + 1 > vc_cap[j])
1175 { /* enlarge the j-th column */
1176 if (luf_enlarge_col(luf, j, vc_len[j] + 10))
1177 { /* overflow of the sparse vector area */
1178 ret = 1;
1179 goto done;
1180 }
1181 /* defragmentation may change row and column pointers of
1182 the matrix V */
1183 p_beg = vr_ptr[p];
1184 p_end = p_beg + vr_len[p] - 1;
1185 q_beg = vc_ptr[q];
1186 q_end = q_beg + vc_len[q] - 1;
1187 }
1188 /* add new non-zero v[i,j] to the j-th column */
1189 j_ptr = vc_ptr[j] + vc_len[j];
1190 sv_ind[j_ptr] = i;
1191 vc_len[j]++;
1192 }
1193 /* now the i-th row has been completely transformed, therefore
1194 it can return to the active set with new length */
1195 rs_prev[i] = 0;
1196 rs_next[i] = rs_head[vr_len[i]];
1197 if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
1198 rs_head[vr_len[i]] = i;
1199 /* the largest of absolute values of elements in the i-th row
1200 is currently unknown */
1201 vr_max[i] = -1.0;
1202 /* at least one free location is needed to store the gaussian
1203 multiplier */
1204 if (luf->sv_end - luf->sv_beg < 1)
1205 { /* there are no free locations at all; defragment SVA */
1206 luf_defrag_sva(luf);
1207 if (luf->sv_end - luf->sv_beg < 1)
1208 { /* overflow of the sparse vector area */
1209 ret = 1;
1210 goto done;
1211 }
1212 /* defragmentation may change row and column pointers of the
1213 matrix V */
1214 p_beg = vr_ptr[p];
1215 p_end = p_beg + vr_len[p] - 1;
1216 q_beg = vc_ptr[q];
1217 q_end = q_beg + vc_len[q] - 1;
1218 }
1219 /* add the element f[i,p], which is the gaussian multiplier,
1220 to the matrix F */
1221 luf->sv_end--;
1222 sv_ind[luf->sv_end] = i;
1223 sv_val[luf->sv_end] = fip;
1224 fc_len[p]++;
1225 /* end of elimination loop */
1226 }
1227 /* at this point the q-th (pivot) column should be empty */
1228 xassert(vc_len[q] == 0);
1229 /* reset capacity of the q-th column */
1230 vc_cap[q] = 0;
1231 /* remove node of the q-th column from the addressing list */
1232 k = n + q;
1233 if (sv_prev[k] == 0)
1234 luf->sv_head = sv_next[k];
1235 else
1236 sv_next[sv_prev[k]] = sv_next[k];
1237 if (sv_next[k] == 0)
1238 luf->sv_tail = sv_prev[k];
1239 else
1240 sv_prev[sv_next[k]] = sv_prev[k];
1241 /* the p-th column of the matrix F has been completely built; set
1242 its pointer */
1243 fc_ptr[p] = luf->sv_end;
1244 /* walk through the p-th (pivot) row and do the following... */
1245 for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
1246 { /* get column index of v[p,j] */
1247 j = sv_ind[p_ptr];
1248 /* erase v[p,j] from the working array */
1249 flag[j] = 0;
1250 work[j] = 0.0;
1251 /* the j-th column has been completely transformed, therefore
1252 it can return to the active set with new length; however
1253 the special case c_prev[j] = c_next[j] = j means that the
1254 routine find_pivot excluded the j-th column from the active
1255 set due to Uwe Suhl's rule, and therefore in this case the
1256 column can return to the active set only if it is a column
1257 singleton */
1258 if (!(vc_len[j] != 1 && cs_prev[j] == j && cs_next[j] == j))
1259 { cs_prev[j] = 0;
1260 cs_next[j] = cs_head[vc_len[j]];
1261 if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
1262 cs_head[vc_len[j]] = j;
1263 }
1264 }
1265 done: /* return to the factorizing routine */
1266 return ret;
1267 }
1269 /***********************************************************************
1270 * build_v_cols - build the matrix V in column-wise format
1271 *
1272 * This routine builds the column-wise representation of the matrix V
1273 * using its row-wise representation.
1274 *
1275 * If no error occured, the routine returns zero. Otherwise, in case of
1276 * overflow of the sparse vector area, the routine returns non-zero. */
1278 static int build_v_cols(LUF *luf)
1279 { int n = luf->n;
1280 int *vr_ptr = luf->vr_ptr;
1281 int *vr_len = luf->vr_len;
1282 int *vc_ptr = luf->vc_ptr;
1283 int *vc_len = luf->vc_len;
1284 int *vc_cap = luf->vc_cap;
1285 int *sv_ind = luf->sv_ind;
1286 double *sv_val = luf->sv_val;
1287 int *sv_prev = luf->sv_prev;
1288 int *sv_next = luf->sv_next;
1289 int ret = 0;
1290 int i, i_beg, i_end, i_ptr, j, j_ptr, k, nnz;
1291 /* it is assumed that on entry all columns of the matrix V are
1292 empty, i.e. vc_len[j] = vc_cap[j] = 0 for all j = 1, ..., n,
1293 and have been removed from the addressing list */
1294 /* count non-zeros in columns of the matrix V; count total number
1295 of non-zeros in this matrix */
1296 nnz = 0;
1297 for (i = 1; i <= n; i++)
1298 { /* walk through elements of the i-th row and count non-zeros
1299 in the corresponding columns */
1300 i_beg = vr_ptr[i];
1301 i_end = i_beg + vr_len[i] - 1;
1302 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
1303 vc_cap[sv_ind[i_ptr]]++;
1304 /* count total number of non-zeros */
1305 nnz += vr_len[i];
1306 }
1307 /* store total number of non-zeros */
1308 luf->nnz_v = nnz;
1309 /* check for free locations */
1310 if (luf->sv_end - luf->sv_beg < nnz)
1311 { /* overflow of the sparse vector area */
1312 ret = 1;
1313 goto done;
1314 }
1315 /* allocate columns of the matrix V */
1316 for (j = 1; j <= n; j++)
1317 { /* set pointer to the j-th column */
1318 vc_ptr[j] = luf->sv_beg;
1319 /* reserve locations for the j-th column */
1320 luf->sv_beg += vc_cap[j];
1321 }
1322 /* build the matrix V in column-wise format using this matrix in
1323 row-wise format */
1324 for (i = 1; i <= n; i++)
1325 { /* walk through elements of the i-th row */
1326 i_beg = vr_ptr[i];
1327 i_end = i_beg + vr_len[i] - 1;
1328 for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
1329 { /* get column index */
1330 j = sv_ind[i_ptr];
1331 /* store element in the j-th column */
1332 j_ptr = vc_ptr[j] + vc_len[j];
1333 sv_ind[j_ptr] = i;
1334 sv_val[j_ptr] = sv_val[i_ptr];
1335 /* increase length of the j-th column */
1336 vc_len[j]++;
1337 }
1338 }
1339 /* now columns are placed in the sparse vector area behind rows
1340 in the order n+1, n+2, ..., n+n; so insert column nodes in the
1341 addressing list using this order */
1342 for (k = n+1; k <= n+n; k++)
1343 { sv_prev[k] = k-1;
1344 sv_next[k] = k+1;
1345 }
1346 sv_prev[n+1] = luf->sv_tail;
1347 sv_next[luf->sv_tail] = n+1;
1348 sv_next[n+n] = 0;
1349 luf->sv_tail = n+n;
1350 done: /* return to the factorizing routine */
1351 return ret;
1352 }
1354 /***********************************************************************
1355 * build_f_rows - build the matrix F in row-wise format
1356 *
1357 * This routine builds the row-wise representation of the matrix F using
1358 * its column-wise representation.
1359 *
1360 * If no error occured, the routine returns zero. Otherwise, in case of
1361 * overflow of the sparse vector area, the routine returns non-zero. */
1363 static int build_f_rows(LUF *luf)
1364 { int n = luf->n;
1365 int *fr_ptr = luf->fr_ptr;
1366 int *fr_len = luf->fr_len;
1367 int *fc_ptr = luf->fc_ptr;
1368 int *fc_len = luf->fc_len;
1369 int *sv_ind = luf->sv_ind;
1370 double *sv_val = luf->sv_val;
1371 int ret = 0;
1372 int i, j, j_beg, j_end, j_ptr, ptr, nnz;
1373 /* clear rows of the matrix F */
1374 for (i = 1; i <= n; i++) fr_len[i] = 0;
1375 /* count non-zeros in rows of the matrix F; count total number of
1376 non-zeros in this matrix */
1377 nnz = 0;
1378 for (j = 1; j <= n; j++)
1379 { /* walk through elements of the j-th column and count non-zeros
1380 in the corresponding rows */
1381 j_beg = fc_ptr[j];
1382 j_end = j_beg + fc_len[j] - 1;
1383 for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
1384 fr_len[sv_ind[j_ptr]]++;
1385 /* increase total number of non-zeros */
1386 nnz += fc_len[j];
1387 }
1388 /* store total number of non-zeros */
1389 luf->nnz_f = nnz;
1390 /* check for free locations */
1391 if (luf->sv_end - luf->sv_beg < nnz)
1392 { /* overflow of the sparse vector area */
1393 ret = 1;
1394 goto done;
1395 }
1396 /* allocate rows of the matrix F */
1397 for (i = 1; i <= n; i++)
1398 { /* set pointer to the end of the i-th row; later this pointer
1399 will be set to the beginning of the i-th row */
1400 fr_ptr[i] = luf->sv_end;
1401 /* reserve locations for the i-th row */
1402 luf->sv_end -= fr_len[i];
1403 }
1404 /* build the matrix F in row-wise format using this matrix in
1405 column-wise format */
1406 for (j = 1; j <= n; j++)
1407 { /* walk through elements of the j-th column */
1408 j_beg = fc_ptr[j];
1409 j_end = j_beg + fc_len[j] - 1;
1410 for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
1411 { /* get row index */
1412 i = sv_ind[j_ptr];
1413 /* store element in the i-th row */
1414 ptr = --fr_ptr[i];
1415 sv_ind[ptr] = j;
1416 sv_val[ptr] = sv_val[j_ptr];
1417 }
1418 }
1419 done: /* return to the factorizing routine */
1420 return ret;
1421 }
1423 /***********************************************************************
1424 * NAME
1425 *
1426 * luf_factorize - compute LU-factorization
1427 *
1428 * SYNOPSIS
1429 *
1430 * #include "glpluf.h"
1431 * int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
1432 * int ind[], double val[]), void *info);
1433 *
1434 * DESCRIPTION
1435 *
1436 * The routine luf_factorize computes LU-factorization of a specified
1437 * square matrix A.
1438 *
1439 * The parameter luf specifies LU-factorization program object created
1440 * by the routine luf_create_it.
1441 *
1442 * The parameter n specifies the order of A, n > 0.
1443 *
1444 * The formal routine col specifies the matrix A to be factorized. To
1445 * obtain j-th column of A the routine luf_factorize calls the routine
1446 * col with the parameter j (1 <= j <= n). In response the routine col
1447 * should store row indices and numerical values of non-zero elements
1448 * of j-th column of A to locations ind[1,...,len] and val[1,...,len],
1449 * respectively, where len is the number of non-zeros in j-th column
1450 * returned on exit. Neither zero nor duplicate elements are allowed.
1451 *
1452 * The parameter info is a transit pointer passed to the routine col.
1453 *
1454 * RETURNS
1455 *
1456 * 0 LU-factorization has been successfully computed.
1457 *
1458 * LUF_ESING
1459 * The specified matrix is singular within the working precision.
1460 * (On some elimination step the active submatrix is exactly zero,
1461 * so no pivot can be chosen.)
1462 *
1463 * LUF_ECOND
1464 * The specified matrix is ill-conditioned.
1465 * (On some elimination step too intensive growth of elements of the
1466 * active submatix has been detected.)
1467 *
1468 * If matrix A is well scaled, the return code LUF_ECOND may also mean
1469 * that the threshold pivoting tolerance piv_tol should be increased.
1470 *
1471 * In case of non-zero return code the factorization becomes invalid.
1472 * It should not be used in other operations until the cause of failure
1473 * has been eliminated and the factorization has been recomputed again
1474 * with the routine luf_factorize.
1475 *
1476 * REPAIRING SINGULAR MATRIX
1477 *
1478 * If the routine luf_factorize returns non-zero code, it provides all
1479 * necessary information that can be used for "repairing" the matrix A,
1480 * where "repairing" means replacing linearly dependent columns of the
1481 * matrix A by appropriate columns of the unity matrix. This feature is
1482 * needed when this routine is used for factorizing the basis matrix
1483 * within the simplex method procedure.
1484 *
1485 * On exit linearly dependent columns of the (partially transformed)
1486 * matrix U have numbers rank+1, rank+2, ..., n, where rank is estimated
1487 * rank of the matrix A stored by the routine to the member luf->rank.
1488 * The correspondence between columns of A and U is the same as between
1489 * columns of V and U. Thus, linearly dependent columns of the matrix A
1490 * have numbers qq_col[rank+1], qq_col[rank+2], ..., qq_col[n], where
1491 * qq_col is the column-like representation of the permutation matrix Q.
1492 * It is understood that each j-th linearly dependent column of the
1493 * matrix U should be replaced by the unity vector, where all elements
1494 * are zero except the unity diagonal element u[j,j]. On the other hand
1495 * j-th row of the matrix U corresponds to the row of the matrix V (and
1496 * therefore of the matrix A) with the number pp_row[j], where pp_row is
1497 * the row-like representation of the permutation matrix P. Thus, each
1498 * j-th linearly dependent column of the matrix U should be replaced by
1499 * column of the unity matrix with the number pp_row[j].
1500 *
1501 * The code that repairs the matrix A may look like follows:
1502 *
1503 * for (j = rank+1; j <= n; j++)
1504 * { replace the column qq_col[j] of the matrix A by the column
1505 * pp_row[j] of the unity matrix;
1506 * }
1507 *
1508 * where rank, pp_row, and qq_col are members of the structure LUF. */
1510 int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
1511 int ind[], double val[]), void *info)
1512 { int *pp_row, *pp_col, *qq_row, *qq_col;
1513 double max_gro = luf->max_gro;
1514 int i, j, k, p, q, t, ret;
1515 if (n < 1)
1516 xfault("luf_factorize: n = %d; invalid parameter\n", n);
1517 if (n > N_MAX)
1518 xfault("luf_factorize: n = %d; matrix too big\n", n);
1519 /* invalidate the factorization */
1520 luf->valid = 0;
1521 /* reallocate arrays, if necessary */
1522 reallocate(luf, n);
1523 pp_row = luf->pp_row;
1524 pp_col = luf->pp_col;
1525 qq_row = luf->qq_row;
1526 qq_col = luf->qq_col;
1527 /* estimate initial size of the SVA, if not specified */
1528 if (luf->sv_size == 0 && luf->new_sva == 0)
1529 luf->new_sva = 5 * (n + 10);
1530 more: /* reallocate the sparse vector area, if required */
1531 if (luf->new_sva > 0)
1532 { if (luf->sv_ind != NULL) xfree(luf->sv_ind);
1533 if (luf->sv_val != NULL) xfree(luf->sv_val);
1534 luf->sv_size = luf->new_sva;
1535 luf->sv_ind = xcalloc(1+luf->sv_size, sizeof(int));
1536 luf->sv_val = xcalloc(1+luf->sv_size, sizeof(double));
1537 luf->new_sva = 0;
1538 }
1539 /* initialize LU-factorization data structures */
1540 if (initialize(luf, col, info))
1541 { /* overflow of the sparse vector area */
1542 luf->new_sva = luf->sv_size + luf->sv_size;
1543 xassert(luf->new_sva > luf->sv_size);
1544 goto more;
1545 }
1546 /* main elimination loop */
1547 for (k = 1; k <= n; k++)
1548 { /* choose a pivot element v[p,q] */
1549 if (find_pivot(luf, &p, &q))
1550 { /* no pivot can be chosen, because the active submatrix is
1551 exactly zero */
1552 luf->rank = k - 1;
1553 ret = LUF_ESING;
1554 goto done;
1555 }
1556 /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th
1557 rows and k-th and j'-th columns of the matrix U = P*V*Q to
1558 move the element u[i',j'] to the position u[k,k] */
1559 i = pp_col[p], j = qq_row[q];
1560 xassert(k <= i && i <= n && k <= j && j <= n);
1561 /* permute k-th and i-th rows of the matrix U */
1562 t = pp_row[k];
1563 pp_row[i] = t, pp_col[t] = i;
1564 pp_row[k] = p, pp_col[p] = k;
1565 /* permute k-th and j-th columns of the matrix U */
1566 t = qq_col[k];
1567 qq_col[j] = t, qq_row[t] = j;
1568 qq_col[k] = q, qq_row[q] = k;
1569 /* eliminate subdiagonal elements of k-th column of the matrix
1570 U = P*V*Q using the pivot element u[k,k] = v[p,q] */
1571 if (eliminate(luf, p, q))
1572 { /* overflow of the sparse vector area */
1573 luf->new_sva = luf->sv_size + luf->sv_size;
1574 xassert(luf->new_sva > luf->sv_size);
1575 goto more;
1576 }
1577 /* check relative growth of elements of the matrix V */
1578 if (luf->big_v > max_gro * luf->max_a)
1579 { /* the growth is too intensive, therefore most probably the
1580 matrix A is ill-conditioned */
1581 luf->rank = k - 1;
1582 ret = LUF_ECOND;
1583 goto done;
1584 }
1585 }
1586 /* now the matrix U = P*V*Q is upper triangular, the matrix V has
1587 been built in row-wise format, and the matrix F has been built
1588 in column-wise format */
1589 /* defragment the sparse vector area in order to merge all free
1590 locations in one continuous extent */
1591 luf_defrag_sva(luf);
1592 /* build the matrix V in column-wise format */
1593 if (build_v_cols(luf))
1594 { /* overflow of the sparse vector area */
1595 luf->new_sva = luf->sv_size + luf->sv_size;
1596 xassert(luf->new_sva > luf->sv_size);
1597 goto more;
1598 }
1599 /* build the matrix F in row-wise format */
1600 if (build_f_rows(luf))
1601 { /* overflow of the sparse vector area */
1602 luf->new_sva = luf->sv_size + luf->sv_size;
1603 xassert(luf->new_sva > luf->sv_size);
1604 goto more;
1605 }
1606 /* the LU-factorization has been successfully computed */
1607 luf->valid = 1;
1608 luf->rank = n;
1609 ret = 0;
1610 /* if there are few free locations in the sparse vector area, try
1611 increasing its size in the future */
1612 t = 3 * (n + luf->nnz_v) + 2 * luf->nnz_f;
1613 if (luf->sv_size < t)
1614 { luf->new_sva = luf->sv_size;
1615 while (luf->new_sva < t)
1616 { k = luf->new_sva;
1617 luf->new_sva = k + k;
1618 xassert(luf->new_sva > k);
1619 }
1620 }
1621 done: /* return to the calling program */
1622 return ret;
1623 }
1625 /***********************************************************************
1626 * NAME
1627 *
1628 * luf_f_solve - solve system F*x = b or F'*x = b
1629 *
1630 * SYNOPSIS
1631 *
1632 * #include "glpluf.h"
1633 * void luf_f_solve(LUF *luf, int tr, double x[]);
1634 *
1635 * DESCRIPTION
1636 *
1637 * The routine luf_f_solve solves either the system F*x = b (if the
1638 * flag tr is zero) or the system F'*x = b (if the flag tr is non-zero),
1639 * where the matrix F is a component of LU-factorization specified by
1640 * the parameter luf, F' is a matrix transposed to F.
1641 *
1642 * On entry the array x should contain elements of the right-hand side
1643 * vector b in locations x[1], ..., x[n], where n is the order of the
1644 * matrix F. On exit this array will contain elements of the solution
1645 * vector x in the same locations. */
1647 void luf_f_solve(LUF *luf, int tr, double x[])
1648 { int n = luf->n;
1649 int *fr_ptr = luf->fr_ptr;
1650 int *fr_len = luf->fr_len;
1651 int *fc_ptr = luf->fc_ptr;
1652 int *fc_len = luf->fc_len;
1653 int *pp_row = luf->pp_row;
1654 int *sv_ind = luf->sv_ind;
1655 double *sv_val = luf->sv_val;
1656 int i, j, k, beg, end, ptr;
1657 double xk;
1658 if (!luf->valid)
1659 xfault("luf_f_solve: LU-factorization is not valid\n");
1660 if (!tr)
1661 { /* solve the system F*x = b */
1662 for (j = 1; j <= n; j++)
1663 { k = pp_row[j];
1664 xk = x[k];
1665 if (xk != 0.0)
1666 { beg = fc_ptr[k];
1667 end = beg + fc_len[k] - 1;
1668 for (ptr = beg; ptr <= end; ptr++)
1669 x[sv_ind[ptr]] -= sv_val[ptr] * xk;
1670 }
1671 }
1672 }
1673 else
1674 { /* solve the system F'*x = b */
1675 for (i = n; i >= 1; i--)
1676 { k = pp_row[i];
1677 xk = x[k];
1678 if (xk != 0.0)
1679 { beg = fr_ptr[k];
1680 end = beg + fr_len[k] - 1;
1681 for (ptr = beg; ptr <= end; ptr++)
1682 x[sv_ind[ptr]] -= sv_val[ptr] * xk;
1683 }
1684 }
1685 }
1686 return;
1687 }
1689 /***********************************************************************
1690 * NAME
1691 *
1692 * luf_v_solve - solve system V*x = b or V'*x = b
1693 *
1694 * SYNOPSIS
1695 *
1696 * #include "glpluf.h"
1697 * void luf_v_solve(LUF *luf, int tr, double x[]);
1698 *
1699 * DESCRIPTION
1700 *
1701 * The routine luf_v_solve solves either the system V*x = b (if the
1702 * flag tr is zero) or the system V'*x = b (if the flag tr is non-zero),
1703 * where the matrix V is a component of LU-factorization specified by
1704 * the parameter luf, V' is a matrix transposed to V.
1705 *
1706 * On entry the array x should contain elements of the right-hand side
1707 * vector b in locations x[1], ..., x[n], where n is the order of the
1708 * matrix V. On exit this array will contain elements of the solution
1709 * vector x in the same locations. */
1711 void luf_v_solve(LUF *luf, int tr, double x[])
1712 { int n = luf->n;
1713 int *vr_ptr = luf->vr_ptr;
1714 int *vr_len = luf->vr_len;
1715 double *vr_piv = luf->vr_piv;
1716 int *vc_ptr = luf->vc_ptr;
1717 int *vc_len = luf->vc_len;
1718 int *pp_row = luf->pp_row;
1719 int *qq_col = luf->qq_col;
1720 int *sv_ind = luf->sv_ind;
1721 double *sv_val = luf->sv_val;
1722 double *b = luf->work;
1723 int i, j, k, beg, end, ptr;
1724 double temp;
1725 if (!luf->valid)
1726 xfault("luf_v_solve: LU-factorization is not valid\n");
1727 for (k = 1; k <= n; k++) b[k] = x[k], x[k] = 0.0;
1728 if (!tr)
1729 { /* solve the system V*x = b */
1730 for (k = n; k >= 1; k--)
1731 { i = pp_row[k], j = qq_col[k];
1732 temp = b[i];
1733 if (temp != 0.0)
1734 { x[j] = (temp /= vr_piv[i]);
1735 beg = vc_ptr[j];
1736 end = beg + vc_len[j] - 1;
1737 for (ptr = beg; ptr <= end; ptr++)
1738 b[sv_ind[ptr]] -= sv_val[ptr] * temp;
1739 }
1740 }
1741 }
1742 else
1743 { /* solve the system V'*x = b */
1744 for (k = 1; k <= n; k++)
1745 { i = pp_row[k], j = qq_col[k];
1746 temp = b[j];
1747 if (temp != 0.0)
1748 { x[i] = (temp /= vr_piv[i]);
1749 beg = vr_ptr[i];
1750 end = beg + vr_len[i] - 1;
1751 for (ptr = beg; ptr <= end; ptr++)
1752 b[sv_ind[ptr]] -= sv_val[ptr] * temp;
1753 }
1754 }
1755 }
1756 return;
1757 }
1759 /***********************************************************************
1760 * NAME
1761 *
1762 * luf_a_solve - solve system A*x = b or A'*x = b
1763 *
1764 * SYNOPSIS
1765 *
1766 * #include "glpluf.h"
1767 * void luf_a_solve(LUF *luf, int tr, double x[]);
1768 *
1769 * DESCRIPTION
1770 *
1771 * The routine luf_a_solve solves either the system A*x = b (if the
1772 * flag tr is zero) or the system A'*x = b (if the flag tr is non-zero),
1773 * where the parameter luf specifies LU-factorization of the matrix A,
1774 * A' is a matrix transposed to A.
1775 *
1776 * On entry the array x should contain elements of the right-hand side
1777 * vector b in locations x[1], ..., x[n], where n is the order of the
1778 * matrix A. On exit this array will contain elements of the solution
1779 * vector x in the same locations. */
1781 void luf_a_solve(LUF *luf, int tr, double x[])
1782 { if (!luf->valid)
1783 xfault("luf_a_solve: LU-factorization is not valid\n");
1784 if (!tr)
1785 { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
1786 luf_f_solve(luf, 0, x);
1787 luf_v_solve(luf, 0, x);
1788 }
1789 else
1790 { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
1791 luf_v_solve(luf, 1, x);
1792 luf_f_solve(luf, 1, x);
1793 }
1794 return;
1795 }
1797 /***********************************************************************
1798 * NAME
1799 *
1800 * luf_delete_it - delete LU-factorization
1801 *
1802 * SYNOPSIS
1803 *
1804 * #include "glpluf.h"
1805 * void luf_delete_it(LUF *luf);
1806 *
1807 * DESCRIPTION
1808 *
1809 * The routine luf_delete deletes LU-factorization specified by the
1810 * parameter luf and frees all the memory allocated to this program
1811 * object. */
1813 void luf_delete_it(LUF *luf)
1814 { if (luf->fr_ptr != NULL) xfree(luf->fr_ptr);
1815 if (luf->fr_len != NULL) xfree(luf->fr_len);
1816 if (luf->fc_ptr != NULL) xfree(luf->fc_ptr);
1817 if (luf->fc_len != NULL) xfree(luf->fc_len);
1818 if (luf->vr_ptr != NULL) xfree(luf->vr_ptr);
1819 if (luf->vr_len != NULL) xfree(luf->vr_len);
1820 if (luf->vr_cap != NULL) xfree(luf->vr_cap);
1821 if (luf->vr_piv != NULL) xfree(luf->vr_piv);
1822 if (luf->vc_ptr != NULL) xfree(luf->vc_ptr);
1823 if (luf->vc_len != NULL) xfree(luf->vc_len);
1824 if (luf->vc_cap != NULL) xfree(luf->vc_cap);
1825 if (luf->pp_row != NULL) xfree(luf->pp_row);
1826 if (luf->pp_col != NULL) xfree(luf->pp_col);
1827 if (luf->qq_row != NULL) xfree(luf->qq_row);
1828 if (luf->qq_col != NULL) xfree(luf->qq_col);
1829 if (luf->sv_ind != NULL) xfree(luf->sv_ind);
1830 if (luf->sv_val != NULL) xfree(luf->sv_val);
1831 if (luf->sv_prev != NULL) xfree(luf->sv_prev);
1832 if (luf->sv_next != NULL) xfree(luf->sv_next);
1833 if (luf->vr_max != NULL) xfree(luf->vr_max);
1834 if (luf->rs_head != NULL) xfree(luf->rs_head);
1835 if (luf->rs_prev != NULL) xfree(luf->rs_prev);
1836 if (luf->rs_next != NULL) xfree(luf->rs_next);
1837 if (luf->cs_head != NULL) xfree(luf->cs_head);
1838 if (luf->cs_prev != NULL) xfree(luf->cs_prev);
1839 if (luf->cs_next != NULL) xfree(luf->cs_next);
1840 if (luf->flag != NULL) xfree(luf->flag);
1841 if (luf->work != NULL) xfree(luf->work);
1842 xfree(luf);
1843 return;
1844 }
1846 /* eof */