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 
55 template <typename Entry, typename Int>
56 Int 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 = EMPTY ;
97  Common->singular_col = 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 == 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 == 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
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.
Generated by   doxygen 1.8.5