Amesos2 - Direct Sparse Solver Interfaces Version of the Day
klu2_extract.hpp
1/* ========================================================================== */
2/* === KLU_extract ========================================================== */
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/* Extract KLU factorization into conventional compressed-column matrices.
38 * If any output array is NULL, that part of the LU factorization is not
39 * extracted (this is not an error condition).
40 *
41 * nnz(L) = Numeric->lnz, nnz(U) = Numeric->unz, and nnz(F) = Numeric->Offp [n]
42 */
43
44#ifndef KLU2_EXTRACT_HPP
45#define KLU2_EXTRACT_HPP
46
47#include "klu2_internal.h"
48
49template <typename Entry, typename Int>
50Int KLU_extract /* returns TRUE if successful, FALSE otherwise */
51(
52 /* inputs: */
53 KLU_numeric<Entry, Int> *Numeric,
54 KLU_symbolic<Entry, Int> *Symbolic,
55
56 /* outputs, all of which must be allocated on input */
57
58 /* L */
59 Int *Lp, /* size n+1 */
60 Int *Li, /* size nnz(L) */
61 double *Lx, /* size nnz(L) */
62#ifdef COMPLEX
63 double *Lz, /* size nnz(L) for the complex case, ignored if real */
64#endif
65
66 /* U */
67 Int *Up, /* size n+1 */
68 Int *Ui, /* size nnz(U) */
69 double *Ux, /* size nnz(U) */
70#ifdef COMPLEX
71 double *Uz, /* size nnz(U) for the complex case, ignored if real */
72#endif
73
74 /* F */
75 Int *Fp, /* size n+1 */
76 Int *Fi, /* size nnz(F) */
77 double *Fx, /* size nnz(F) */
78#ifdef COMPLEX
79 double *Fz, /* size nnz(F) for the complex case, ignored if real */
80#endif
81
82 /* P, row permutation */
83 Int *P, /* size n */
84
85 /* Q, column permutation */
86 Int *Q, /* size n */
87
88 /* Rs, scale factors */
89 double *Rs, /* size n */
90
91 /* R, block boundaries */
92 Int *R, /* size nblocks+1 */
93
94 KLU_common<Entry> *Common
95)
96{
97 Int *Lip, *Llen, *Uip, *Ulen, *Li2, *Ui2 ;
98 Unit *LU ;
99 Entry *Lx2, *Ux2, *Ukk ;
100 Int i, k, block, nblocks, n, nz, k1, k2, nk, len, kk, p ;
101
102 if (Common == NULL)
103 {
104 return (FALSE) ;
105 }
106
107 if (Symbolic == NULL || Numeric == NULL)
108 {
109 Common->status = KLU_INVALID ;
110 return (FALSE) ;
111 }
112
113 Common->status = KLU_OK ;
114 n = Symbolic->n ;
115 nblocks = Symbolic->nblocks ;
116
117 /* ---------------------------------------------------------------------- */
118 /* extract scale factors */
119 /* ---------------------------------------------------------------------- */
120
121 if (Rs != NULL)
122 {
123 if (Numeric->Rs != NULL)
124 {
125 for (i = 0 ; i < n ; i++)
126 {
127 Rs [i] = Numeric->Rs [i] ;
128 }
129 }
130 else
131 {
132 /* no scaling */
133 for (i = 0 ; i < n ; i++)
134 {
135 Rs [i] = 1 ;
136 }
137 }
138 }
139
140 /* ---------------------------------------------------------------------- */
141 /* extract block boundaries */
142 /* ---------------------------------------------------------------------- */
143
144 if (R != NULL)
145 {
146 for (block = 0 ; block <= nblocks ; block++)
147 {
148 R [block] = Symbolic->R [block] ;
149 }
150 }
151
152 /* ---------------------------------------------------------------------- */
153 /* extract final row permutation */
154 /* ---------------------------------------------------------------------- */
155
156 if (P != NULL)
157 {
158 for (k = 0 ; k < n ; k++)
159 {
160 P [k] = Numeric->Pnum [k] ;
161 }
162 }
163
164 /* ---------------------------------------------------------------------- */
165 /* extract column permutation */
166 /* ---------------------------------------------------------------------- */
167
168 if (Q != NULL)
169 {
170 for (k = 0 ; k < n ; k++)
171 {
172 Q [k] = Symbolic->Q [k] ;
173 }
174 }
175
176 /* ---------------------------------------------------------------------- */
177 /* extract each block of L */
178 /* ---------------------------------------------------------------------- */
179
180 if (Lp != NULL && Li != NULL && Lx != NULL
181#ifdef COMPLEX
182 && Lz != NULL
183#endif
184 )
185 {
186 nz = 0 ;
187 for (block = 0 ; block < nblocks ; block++)
188 {
189 k1 = Symbolic->R [block] ;
190 k2 = Symbolic->R [block+1] ;
191 nk = k2 - k1 ;
192 if (nk == 1)
193 {
194 /* singleton block */
195 Lp [k1] = nz ;
196 Li [nz] = k1 ;
197 Lx [nz] = 1 ;
198#ifdef COMPLEX
199 Lz [nz] = 0 ;
200#endif
201 nz++ ;
202 }
203 else
204 {
205 /* non-singleton block */
206 LU = (Unit *) Numeric->LUbx [block] ;
207 Lip = Numeric->Lip + k1 ;
208 Llen = Numeric->Llen + k1 ;
209 for (kk = 0 ; kk < nk ; kk++)
210 {
211 Lp [k1+kk] = nz ;
212 /* add the unit diagonal entry */
213 Li [nz] = k1 + kk ;
214 Lx [nz] = 1 ;
215#ifdef COMPLEX
216 Lz [nz] = 0 ;
217#endif
218 nz++ ;
219 GET_POINTER (LU, Lip, Llen, Li2, Lx2, kk, len) ;
220 for (p = 0 ; p < len ; p++)
221 {
222 Li [nz] = k1 + Li2 [p] ;
223 Lx [nz] = REAL (Lx2 [p]) ;
224#ifdef COMPLEX
225 Lz [nz] = IMAG (Lx2 [p]) ;
226#endif
227 nz++ ;
228 }
229 }
230 }
231 }
232 Lp [n] = nz ;
233 ASSERT (nz == Numeric->lnz) ;
234 }
235
236 /* ---------------------------------------------------------------------- */
237 /* extract each block of U */
238 /* ---------------------------------------------------------------------- */
239
240 if (Up != NULL && Ui != NULL && Ux != NULL
241#ifdef COMPLEX
242 && Uz != NULL
243#endif
244 )
245 {
246 nz = 0 ;
247 for (block = 0 ; block < nblocks ; block++)
248 {
249 k1 = Symbolic->R [block] ;
250 k2 = Symbolic->R [block+1] ;
251 nk = k2 - k1 ;
252 Ukk = ((Entry *) Numeric->Udiag) + k1 ;
253 if (nk == 1)
254 {
255 /* singleton block */
256 Up [k1] = nz ;
257 Ui [nz] = k1 ;
258 Ux [nz] = REAL (Ukk [0]) ;
259#ifdef COMPLEX
260 Uz [nz] = IMAG (Ukk [0]) ;
261#endif
262 nz++ ;
263 }
264 else
265 {
266 /* non-singleton block */
267 LU = (Unit *) Numeric->LUbx [block] ;
268 Uip = Numeric->Uip + k1 ;
269 Ulen = Numeric->Ulen + k1 ;
270 for (kk = 0 ; kk < nk ; kk++)
271 {
272 Up [k1+kk] = nz ;
273 GET_POINTER (LU, Uip, Ulen, Ui2, Ux2, kk, len) ;
274 for (p = 0 ; p < len ; p++)
275 {
276 Ui [nz] = k1 + Ui2 [p] ;
277 Ux [nz] = REAL (Ux2 [p]) ;
278#ifdef COMPLEX
279 Uz [nz] = IMAG (Ux2 [p]) ;
280#endif
281 nz++ ;
282 }
283 /* add the diagonal entry */
284 Ui [nz] = k1 + kk ;
285 Ux [nz] = REAL (Ukk [kk]) ;
286#ifdef COMPLEX
287 Uz [nz] = IMAG (Ukk [kk]) ;
288#endif
289 nz++ ;
290 }
291 }
292 }
293 Up [n] = nz ;
294 ASSERT (nz == Numeric->unz) ;
295 }
296
297 /* ---------------------------------------------------------------------- */
298 /* extract the off-diagonal blocks, F */
299 /* ---------------------------------------------------------------------- */
300
301 if (Fp != NULL && Fi != NULL && Fx != NULL
302#ifdef COMPLEX
303 && Fz != NULL
304#endif
305 )
306 {
307 for (k = 0 ; k <= n ; k++)
308 {
309 Fp [k] = Numeric->Offp [k] ;
310 }
311 nz = Fp [n] ;
312 for (k = 0 ; k < nz ; k++)
313 {
314 Fi [k] = Numeric->Offi [k] ;
315 }
316 for (k = 0 ; k < nz ; k++)
317 {
318 Fx [k] = REAL (((Entry *) Numeric->Offx) [k]) ;
319#ifdef COMPLEX
320 Fz [k] = IMAG (((Entry *) Numeric->Offx) [k]) ;
321#endif
322 }
323 }
324
325 return (TRUE) ;
326}
327
328#endif