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