|
1 /* glpluf.c (LU-factorization) */ |
|
2 |
|
3 /*********************************************************************** |
|
4 * This code is part of GLPK (GNU Linear Programming Kit). |
|
5 * |
|
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
|
7 * 2009, 2010 Andrew Makhorin, Department for Applied Informatics, |
|
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
|
9 * E-mail: <mao@gnu.org>. |
|
10 * |
|
11 * GLPK is free software: you can redistribute it and/or modify it |
|
12 * under the terms of the GNU General Public License as published by |
|
13 * the Free Software Foundation, either version 3 of the License, or |
|
14 * (at your option) any later version. |
|
15 * |
|
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT |
|
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
|
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
|
19 * License for more details. |
|
20 * |
|
21 * You should have received a copy of the GNU General Public License |
|
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
|
23 ***********************************************************************/ |
|
24 |
|
25 #include "glpenv.h" |
|
26 #include "glpluf.h" |
|
27 #define xfault xerror |
|
28 |
|
29 /* CAUTION: DO NOT CHANGE THE LIMIT BELOW */ |
|
30 |
|
31 #define N_MAX 100000000 /* = 100*10^6 */ |
|
32 /* maximal order of the original matrix */ |
|
33 |
|
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. */ |
|
52 |
|
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 } |
|
87 |
|
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. */ |
|
110 |
|
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 } |
|
171 |
|
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. */ |
|
199 |
|
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 } |
|
262 |
|
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. */ |
|
291 |
|
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 } |
|
354 |
|
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. */ |
|
360 |
|
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 } |
|
420 |
|
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. */ |
|
432 |
|
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 } |
|
623 |
|
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. */ |
|
683 |
|
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 } |
|
879 |
|
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. */ |
|
941 |
|
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 } |
|
1268 |
|
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. */ |
|
1277 |
|
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 } |
|
1353 |
|
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. */ |
|
1362 |
|
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 } |
|
1422 |
|
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. */ |
|
1509 |
|
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 } |
|
1624 |
|
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. */ |
|
1646 |
|
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 } |
|
1688 |
|
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. */ |
|
1710 |
|
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 } |
|
1758 |
|
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. */ |
|
1780 |
|
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 } |
|
1796 |
|
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. */ |
|
1812 |
|
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 } |
|
1845 |
|
1846 /* eof */ |