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 
31 template <typename Entry, typename Int>
32 Int 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.
Ap
int Ap[]
Column offsets.
Definition: klu2_simple.cpp:52
Ai
int Ai[]
Row values.
Definition: klu2_simple.cpp:54
Generated by   doxygen 1.8.5