1 /* ========================================================================== */
2 /* === KLU_refactor ========================================================= */
3 /* ========================================================================== */
4
5 /* Factor the matrix, after ordering and analyzing it with KLU_analyze, and
6 * factoring it once with KLU_factor. This routine cannot do any numerical
7 * pivoting. The pattern of the input matrix (Ap, Ai) must be identical to
8 * the pattern given to KLU_factor.
9 */
10
11 #include "klu_internal.h"
12
13
14 /* ========================================================================== */
15 /* === KLU_refactor ========================================================= */
16 /* ========================================================================== */
17
KLU_refactor(Int Ap[],Int Ai[],double Ax[],KLU_symbolic * Symbolic,KLU_numeric * Numeric,KLU_common * Common)18 Int KLU_refactor /* returns TRUE if successful, FALSE otherwise */
19 (
20 /* inputs, not modified */
21 Int Ap [ ], /* size n+1, column pointers */
22 Int Ai [ ], /* size nz, row indices */
23 double Ax [ ],
24 KLU_symbolic *Symbolic,
25
26 /* input/output */
27 KLU_numeric *Numeric,
28 KLU_common *Common
29 )
30 {
31 Entry ukk, ujk, s ;
32 Entry *Offx, *Lx, *Ux, *X, *Az, *Udiag ;
33 double *Rs ;
34 Int *Q, *R, *Pnum, *Ui, *Li, *Pinv, *Lip, *Uip, *Llen, *Ulen ;
35 Unit **LUbx ;
36 Unit *LU ;
37 Int k1, k2, nk, k, block, oldcol, pend, oldrow, n, p, newrow, scale,
38 nblocks, poff, i, j, up, ulen, llen, maxblock, nzoff ;
39
40 /* ---------------------------------------------------------------------- */
41 /* check inputs */
42 /* ---------------------------------------------------------------------- */
43
44 if (Common == NULL)
45 {
46 return (FALSE) ;
47 }
48 Common->status = KLU_OK ;
49
50 if (Numeric == NULL)
51 {
52 /* invalid Numeric object */
53 Common->status = KLU_INVALID ;
54 return (FALSE) ;
55 }
56
57 Common->numerical_rank = EMPTY ;
58 Common->singular_col = EMPTY ;
59
60 Az = (Entry *) Ax ;
61
62 /* ---------------------------------------------------------------------- */
63 /* get the contents of the Symbolic object */
64 /* ---------------------------------------------------------------------- */
65
66 n = Symbolic->n ;
67 Q = Symbolic->Q ;
68 R = Symbolic->R ;
69 nblocks = Symbolic->nblocks ;
70 maxblock = Symbolic->maxblock ;
71
72 /* ---------------------------------------------------------------------- */
73 /* get the contents of the Numeric object */
74 /* ---------------------------------------------------------------------- */
75
76 Pnum = Numeric->Pnum ;
77 Offx = (Entry *) Numeric->Offx ;
78
79 LUbx = (Unit **) Numeric->LUbx ;
80
81 scale = Common->scale ;
82 if (scale > 0)
83 {
84 /* factorization was not scaled, but refactorization is scaled */
85 if (Numeric->Rs == NULL)
86 {
87 Numeric->Rs = KLU_malloc (n, sizeof (double), Common) ;
88 if (Common->status < KLU_OK)
89 {
90 Common->status = KLU_OUT_OF_MEMORY ;
91 return (FALSE) ;
92 }
93 }
94 }
95 else
96 {
97 /* no scaling for refactorization; ensure Numeric->Rs is freed. This
98 * does nothing if Numeric->Rs is already NULL. */
99 Numeric->Rs = KLU_free (Numeric->Rs, n, sizeof (double), Common) ;
100 }
101 Rs = Numeric->Rs ;
102
103 Pinv = Numeric->Pinv ;
104 X = (Entry *) Numeric->Xwork ;
105 Common->nrealloc = 0 ;
106 Udiag = Numeric->Udiag ;
107 nzoff = Symbolic->nzoff ;
108
109 /* ---------------------------------------------------------------------- */
110 /* check the input matrix compute the row scale factors, Rs */
111 /* ---------------------------------------------------------------------- */
112
113 /* do no scale, or check the input matrix, if scale < 0 */
114 if (scale >= 0)
115 {
116 /* check for out-of-range indices, but do not check for duplicates */
117 if (!KLU_scale (scale, n, Ap, Ai, Ax, Rs, NULL, Common))
118 {
119 return (FALSE) ;
120 }
121 }
122
123 /* ---------------------------------------------------------------------- */
124 /* clear workspace X */
125 /* ---------------------------------------------------------------------- */
126
127 for (k = 0 ; k < maxblock ; k++)
128 {
129 /* X [k] = 0 */
130 CLEAR (X [k]) ;
131 }
132
133 poff = 0 ;
134
135 /* ---------------------------------------------------------------------- */
136 /* factor each block */
137 /* ---------------------------------------------------------------------- */
138
139 if (scale <= 0)
140 {
141
142 /* ------------------------------------------------------------------ */
143 /* no scaling */
144 /* ------------------------------------------------------------------ */
145
146 for (block = 0 ; block < nblocks ; block++)
147 {
148
149 /* -------------------------------------------------------------- */
150 /* the block is from rows/columns k1 to k2-1 */
151 /* -------------------------------------------------------------- */
152
153 k1 = R [block] ;
154 k2 = R [block+1] ;
155 nk = k2 - k1 ;
156
157 if (nk == 1)
158 {
159
160 /* ---------------------------------------------------------- */
161 /* singleton case */
162 /* ---------------------------------------------------------- */
163
164 oldcol = Q [k1] ;
165 pend = Ap [oldcol+1] ;
166 CLEAR (s) ;
167 for (p = Ap [oldcol] ; p < pend ; p++)
168 {
169 newrow = Pinv [Ai [p]] - k1 ;
170 if (newrow < 0 && poff < nzoff)
171 {
172 /* entry in off-diagonal block */
173 Offx [poff] = Az [p] ;
174 poff++ ;
175 }
176 else
177 {
178 /* singleton */
179 s = Az [p] ;
180 }
181 }
182 Udiag [k1] = s ;
183
184 }
185 else
186 {
187
188 /* ---------------------------------------------------------- */
189 /* construct and factor the kth block */
190 /* ---------------------------------------------------------- */
191
192 Lip = Numeric->Lip + k1 ;
193 Llen = Numeric->Llen + k1 ;
194 Uip = Numeric->Uip + k1 ;
195 Ulen = Numeric->Ulen + k1 ;
196 LU = LUbx [block] ;
197
198 for (k = 0 ; k < nk ; k++)
199 {
200
201 /* ------------------------------------------------------ */
202 /* scatter kth column of the block into workspace X */
203 /* ------------------------------------------------------ */
204
205 oldcol = Q [k+k1] ;
206 pend = Ap [oldcol+1] ;
207 for (p = Ap [oldcol] ; p < pend ; p++)
208 {
209 newrow = Pinv [Ai [p]] - k1 ;
210 if (newrow < 0 && poff < nzoff)
211 {
212 /* entry in off-diagonal block */
213 Offx [poff] = Az [p] ;
214 poff++ ;
215 }
216 else
217 {
218 /* (newrow,k) is an entry in the block */
219 X [newrow] = Az [p] ;
220 }
221 }
222
223 /* ------------------------------------------------------ */
224 /* compute kth column of U, and update kth column of A */
225 /* ------------------------------------------------------ */
226
227 GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ;
228 for (up = 0 ; up < ulen ; up++)
229 {
230 j = Ui [up] ;
231 ujk = X [j] ;
232 /* X [j] = 0 */
233 CLEAR (X [j]) ;
234 Ux [up] = ujk ;
235 GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ;
236 for (p = 0 ; p < llen ; p++)
237 {
238 /* X [Li [p]] -= Lx [p] * ujk */
239 MULT_SUB (X [Li [p]], Lx [p], ujk) ;
240 }
241 }
242 /* get the diagonal entry of U */
243 ukk = X [k] ;
244 /* X [k] = 0 */
245 CLEAR (X [k]) ;
246 /* singular case */
247 if (IS_ZERO (ukk))
248 {
249 /* matrix is numerically singular */
250 Common->status = KLU_SINGULAR ;
251 if (Common->numerical_rank == EMPTY)
252 {
253 Common->numerical_rank = k+k1 ;
254 Common->singular_col = Q [k+k1] ;
255 }
256 if (Common->halt_if_singular)
257 {
258 /* do not continue the factorization */
259 return (FALSE) ;
260 }
261 }
262 Udiag [k+k1] = ukk ;
263 /* gather and divide by pivot to get kth column of L */
264 GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ;
265 for (p = 0 ; p < llen ; p++)
266 {
267 i = Li [p] ;
268 DIV (Lx [p], X [i], ukk) ;
269 CLEAR (X [i]) ;
270 }
271
272 }
273 }
274 }
275
276 }
277 else
278 {
279
280 /* ------------------------------------------------------------------ */
281 /* scaling */
282 /* ------------------------------------------------------------------ */
283
284 for (block = 0 ; block < nblocks ; block++)
285 {
286
287 /* -------------------------------------------------------------- */
288 /* the block is from rows/columns k1 to k2-1 */
289 /* -------------------------------------------------------------- */
290
291 k1 = R [block] ;
292 k2 = R [block+1] ;
293 nk = k2 - k1 ;
294
295 if (nk == 1)
296 {
297
298 /* ---------------------------------------------------------- */
299 /* singleton case */
300 /* ---------------------------------------------------------- */
301
302 oldcol = Q [k1] ;
303 pend = Ap [oldcol+1] ;
304 CLEAR (s) ;
305 for (p = Ap [oldcol] ; p < pend ; p++)
306 {
307 oldrow = Ai [p] ;
308 newrow = Pinv [oldrow] - k1 ;
309 if (newrow < 0 && poff < nzoff)
310 {
311 /* entry in off-diagonal block */
312 /* Offx [poff] = Az [p] / Rs [oldrow] */
313 SCALE_DIV_ASSIGN (Offx [poff], Az [p], Rs [oldrow]) ;
314 poff++ ;
315 }
316 else
317 {
318 /* singleton */
319 /* s = Az [p] / Rs [oldrow] */
320 SCALE_DIV_ASSIGN (s, Az [p], Rs [oldrow]) ;
321 }
322 }
323 Udiag [k1] = s ;
324
325 }
326 else
327 {
328
329 /* ---------------------------------------------------------- */
330 /* construct and factor the kth block */
331 /* ---------------------------------------------------------- */
332
333 Lip = Numeric->Lip + k1 ;
334 Llen = Numeric->Llen + k1 ;
335 Uip = Numeric->Uip + k1 ;
336 Ulen = Numeric->Ulen + k1 ;
337 LU = LUbx [block] ;
338
339 for (k = 0 ; k < nk ; k++)
340 {
341
342 /* ------------------------------------------------------ */
343 /* scatter kth column of the block into workspace X */
344 /* ------------------------------------------------------ */
345
346 oldcol = Q [k+k1] ;
347 pend = Ap [oldcol+1] ;
348 for (p = Ap [oldcol] ; p < pend ; p++)
349 {
350 oldrow = Ai [p] ;
351 newrow = Pinv [oldrow] - k1 ;
352 if (newrow < 0 && poff < nzoff)
353 {
354 /* entry in off-diagonal part */
355 /* Offx [poff] = Az [p] / Rs [oldrow] */
356 SCALE_DIV_ASSIGN (Offx [poff], Az [p], Rs [oldrow]);
357 poff++ ;
358 }
359 else
360 {
361 /* (newrow,k) is an entry in the block */
362 /* X [newrow] = Az [p] / Rs [oldrow] */
363 SCALE_DIV_ASSIGN (X [newrow], Az [p], Rs [oldrow]) ;
364 }
365 }
366
367 /* ------------------------------------------------------ */
368 /* compute kth column of U, and update kth column of A */
369 /* ------------------------------------------------------ */
370
371 GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ;
372 for (up = 0 ; up < ulen ; up++)
373 {
374 j = Ui [up] ;
375 ujk = X [j] ;
376 /* X [j] = 0 */
377 CLEAR (X [j]) ;
378 Ux [up] = ujk ;
379 GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ;
380 for (p = 0 ; p < llen ; p++)
381 {
382 /* X [Li [p]] -= Lx [p] * ujk */
383 MULT_SUB (X [Li [p]], Lx [p], ujk) ;
384 }
385 }
386 /* get the diagonal entry of U */
387 ukk = X [k] ;
388 /* X [k] = 0 */
389 CLEAR (X [k]) ;
390 /* singular case */
391 if (IS_ZERO (ukk))
392 {
393 /* matrix is numerically singular */
394 Common->status = KLU_SINGULAR ;
395 if (Common->numerical_rank == EMPTY)
396 {
397 Common->numerical_rank = k+k1 ;
398 Common->singular_col = Q [k+k1] ;
399 }
400 if (Common->halt_if_singular)
401 {
402 /* do not continue the factorization */
403 return (FALSE) ;
404 }
405 }
406 Udiag [k+k1] = ukk ;
407 /* gather and divide by pivot to get kth column of L */
408 GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ;
409 for (p = 0 ; p < llen ; p++)
410 {
411 i = Li [p] ;
412 DIV (Lx [p], X [i], ukk) ;
413 CLEAR (X [i]) ;
414 }
415 }
416 }
417 }
418 }
419
420 /* ---------------------------------------------------------------------- */
421 /* permute scale factors Rs according to pivotal row order */
422 /* ---------------------------------------------------------------------- */
423
424 if (scale > 0)
425 {
426 for (k = 0 ; k < n ; k++)
427 {
428 REAL (X [k]) = Rs [Pnum [k]] ;
429 }
430 for (k = 0 ; k < n ; k++)
431 {
432 Rs [k] = REAL (X [k]) ;
433 }
434 }
435
436 #ifndef NDEBUG
437 ASSERT (Numeric->Offp [n] == poff) ;
438 ASSERT (Symbolic->nzoff == poff) ;
439 PRINTF (("\n------------------- Off diagonal entries, new:\n")) ;
440 ASSERT (KLU_valid (n, Numeric->Offp, Numeric->Offi, Offx)) ;
441 if (Common->status == KLU_OK)
442 {
443 PRINTF (("\n ########### KLU_BTF_REFACTOR done, nblocks %d\n",nblocks));
444 for (block = 0 ; block < nblocks ; block++)
445 {
446 k1 = R [block] ;
447 k2 = R [block+1] ;
448 nk = k2 - k1 ;
449 PRINTF ((
450 "\n================KLU_refactor output: k1 %d k2 %d nk %d\n",
451 k1, k2, nk)) ;
452 if (nk == 1)
453 {
454 PRINTF (("singleton ")) ;
455 PRINT_ENTRY (Udiag [k1]) ;
456 }
457 else
458 {
459 Lip = Numeric->Lip + k1 ;
460 Llen = Numeric->Llen + k1 ;
461 LU = (Unit *) Numeric->LUbx [block] ;
462 PRINTF (("\n---- L block %d\n", block)) ;
463 ASSERT (KLU_valid_LU (nk, TRUE, Lip, Llen, LU)) ;
464 Uip = Numeric->Uip + k1 ;
465 Ulen = Numeric->Ulen + k1 ;
466 PRINTF (("\n---- U block %d\n", block)) ;
467 ASSERT (KLU_valid_LU (nk, FALSE, Uip, Ulen, LU)) ;
468 }
469 }
470 }
471 #endif
472
473 return (TRUE) ;
474 }
475