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