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