1 /* ========================================================================== */
2 /* === UMF_kernel_wrapup ==================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Copyright (c) Timothy A. Davis, CISE, */
7 /* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
8 /* web: http://www.cise.ufl.edu/research/sparse/umfpack */
9 /* -------------------------------------------------------------------------- */
10
11 /* The matrix is factorized. Finish the LU data structure. */
12
13 #include "umf_internal.h"
14 #include "umf_kernel_wrapup.h"
15
UMF_kernel_wrapup(NumericType * Numeric,SymbolicType * Symbolic,WorkType * Work)16 GLOBAL void UMF_kernel_wrapup
17 (
18 NumericType *Numeric,
19 SymbolicType *Symbolic,
20 WorkType *Work
21 )
22 {
23
24 /* ---------------------------------------------------------------------- */
25 /* local variables */
26 /* ---------------------------------------------------------------------- */
27
28 Entry pivot_value ;
29 double d ;
30 Entry *D ;
31 Int i, k, col, row, llen, ulen, *ip, *Rperm, *Cperm, *Lilen, npiv, lp,
32 *Uilen, *Lip, *Uip, *Cperm_init, up, pivrow, pivcol, *Lpos, *Upos, *Wr,
33 *Wc, *Wp, *Frpos, *Fcpos, *Row_degree, *Col_degree, *Rperm_init,
34 n_row, n_col, n_inner, zero_pivot, nan_pivot, n1 ;
35
36 #ifndef NDEBUG
37 UMF_dump_matrix (Numeric, Work, FALSE) ;
38 #endif
39
40 DEBUG0 (("Kernel complete, Starting Kernel wrapup\n")) ;
41 n_row = Symbolic->n_row ;
42 n_col = Symbolic->n_col ;
43 n_inner = MIN (n_row, n_col) ;
44 Rperm = Numeric->Rperm ;
45 Cperm = Numeric->Cperm ;
46 Lilen = Numeric->Lilen ;
47 Uilen = Numeric->Uilen ;
48 Upos = Numeric->Upos ;
49 Lpos = Numeric->Lpos ;
50 Lip = Numeric->Lip ;
51 Uip = Numeric->Uip ;
52 D = Numeric->D ;
53
54 npiv = Work->npiv ;
55 Numeric->npiv = npiv ;
56 Numeric->ulen = Work->ulen ;
57
58 ASSERT (n_row == Numeric->n_row) ;
59 ASSERT (n_col == Symbolic->n_col) ;
60 DEBUG0 (("Wrap-up: npiv "ID" ulen "ID"\n", npiv, Numeric->ulen)) ;
61 ASSERT (npiv <= n_inner) ;
62
63 /* this will be nonzero only if matrix is singular or rectangular */
64 ASSERT (IMPLIES (npiv == n_col, Work->ulen == 0)) ;
65
66 /* ---------------------------------------------------------------------- */
67 /* find the smallest and largest entries in D */
68 /* ---------------------------------------------------------------------- */
69
70 for (k = 0 ; k < npiv ; k++)
71 {
72 pivot_value = D [k] ;
73 ABS (d, pivot_value) ;
74 zero_pivot = SCALAR_IS_ZERO (d) ;
75 nan_pivot = SCALAR_IS_NAN (d) ;
76
77 if (!zero_pivot)
78 {
79 /* the pivot is nonzero, but might be Inf or NaN */
80 Numeric->nnzpiv++ ;
81 }
82
83 if (k == 0)
84 {
85 Numeric->min_udiag = d ;
86 Numeric->max_udiag = d ;
87 }
88 else
89 {
90 /* min (abs (diag (U))) behaves as follows: If any entry is zero,
91 then the result is zero (regardless of the presence of NaN's).
92 Otherwise, if any entry is NaN, then the result is NaN.
93 Otherwise, the result is the smallest absolute value on the
94 diagonal of U.
95 */
96
97 if (SCALAR_IS_NONZERO (Numeric->min_udiag))
98 {
99 if (zero_pivot || nan_pivot)
100 {
101 Numeric->min_udiag = d ;
102 }
103 else if (!SCALAR_IS_NAN (Numeric->min_udiag))
104 {
105 /* d and min_udiag are both non-NaN */
106 Numeric->min_udiag = MIN (Numeric->min_udiag, d) ;
107 }
108 }
109
110 /*
111 max (abs (diag (U))) behaves as follows: If any entry is NaN
112 then the result is NaN. Otherise, the result is the largest
113 absolute value on the diagonal of U.
114 */
115
116 if (nan_pivot)
117 {
118 Numeric->max_udiag = d ;
119 }
120 else if (!SCALAR_IS_NAN (Numeric->max_udiag))
121 {
122 /* d and max_udiag are both non-NaN */
123 Numeric->max_udiag = MAX (Numeric->max_udiag, d) ;
124 }
125 }
126 }
127
128 /* ---------------------------------------------------------------------- */
129 /* check if matrix is singular or rectangular */
130 /* ---------------------------------------------------------------------- */
131
132 Col_degree = Cperm ; /* for NON_PIVOTAL_COL macro */
133 Row_degree = Rperm ; /* for NON_PIVOTAL_ROW macro */
134
135 if (npiv < n_row)
136 {
137 /* finalize the row permutation */
138 k = npiv ;
139 DEBUGm3 (("Singular pivot rows "ID" to "ID"\n", k, n_row-1)) ;
140 for (row = 0 ; row < n_row ; row++)
141 {
142 if (NON_PIVOTAL_ROW (row))
143 {
144 Rperm [row] = ONES_COMPLEMENT (k) ;
145 DEBUGm3 (("Singular row "ID" is k: "ID" pivot row\n", row, k)) ;
146 ASSERT (!NON_PIVOTAL_ROW (row)) ;
147 Lpos [row] = EMPTY ;
148 Uip [row] = EMPTY ;
149 Uilen [row] = 0 ;
150 k++ ;
151 }
152 }
153 ASSERT (k == n_row) ;
154 }
155
156 if (npiv < n_col)
157 {
158 /* finalize the col permutation */
159 k = npiv ;
160 DEBUGm3 (("Singular pivot cols "ID" to "ID"\n", k, n_col-1)) ;
161 for (col = 0 ; col < n_col ; col++)
162 {
163 if (NON_PIVOTAL_COL (col))
164 {
165 Cperm [col] = ONES_COMPLEMENT (k) ;
166 DEBUGm3 (("Singular col "ID" is k: "ID" pivot row\n", col, k)) ;
167 ASSERT (!NON_PIVOTAL_COL (col)) ;
168 Upos [col] = EMPTY ;
169 Lip [col] = EMPTY ;
170 Lilen [col] = 0 ;
171 k++ ;
172 }
173 }
174 ASSERT (k == n_col) ;
175 }
176
177 if (npiv < n_inner)
178 {
179 /* finalize the diagonal of U */
180 DEBUGm3 (("Diag of U is zero, "ID" to "ID"\n", npiv, n_inner-1)) ;
181 for (k = npiv ; k < n_inner ; k++)
182 {
183 CLEAR (D [k]) ;
184 }
185 }
186
187 /* save the pattern of the last row of U */
188 if (Numeric->ulen > 0)
189 {
190 DEBUGm3 (("Last row of U is not empty\n")) ;
191 Numeric->Upattern = Work->Upattern ;
192 Work->Upattern = (Int *) NULL ;
193 }
194
195 DEBUG2 (("Nnzpiv: "ID" npiv "ID"\n", Numeric->nnzpiv, npiv)) ;
196 ASSERT (Numeric->nnzpiv <= npiv) ;
197 if (Numeric->nnzpiv < n_inner && !SCALAR_IS_NAN (Numeric->min_udiag))
198 {
199 /* the rest of the diagonal is zero, so min_udiag becomes 0,
200 * unless it is already NaN. */
201 Numeric->min_udiag = 0.0 ;
202 }
203
204 /* ---------------------------------------------------------------------- */
205 /* size n_row, n_col workspaces that can be used here: */
206 /* ---------------------------------------------------------------------- */
207
208 Frpos = Work->Frpos ; /* of size n_row+1 */
209 Fcpos = Work->Fcpos ; /* of size n_col+1 */
210 Wp = Work->Wp ; /* of size MAX(n_row,n_col)+1 */
211 /* Work->Upattern ; cannot be used (in Numeric) */
212 Wr = Work->Lpattern ; /* of size n_row+1 */
213 Wc = Work->Wrp ; /* of size n_col+1 or bigger */
214
215 /* ---------------------------------------------------------------------- */
216 /* construct Rperm from inverse permutations */
217 /* ---------------------------------------------------------------------- */
218
219 /* use Frpos for temporary copy of inverse row permutation [ */
220
221 for (pivrow = 0 ; pivrow < n_row ; pivrow++)
222 {
223 k = Rperm [pivrow] ;
224 ASSERT (k < 0) ;
225 k = ONES_COMPLEMENT (k) ;
226 ASSERT (k >= 0 && k < n_row) ;
227 Wp [k] = pivrow ;
228 Frpos [pivrow] = k ;
229 }
230 for (k = 0 ; k < n_row ; k++)
231 {
232 Rperm [k] = Wp [k] ;
233 }
234
235 /* ---------------------------------------------------------------------- */
236 /* construct Cperm from inverse permutation */
237 /* ---------------------------------------------------------------------- */
238
239 /* use Fcpos for temporary copy of inverse column permutation [ */
240
241 for (pivcol = 0 ; pivcol < n_col ; pivcol++)
242 {
243 k = Cperm [pivcol] ;
244 ASSERT (k < 0) ;
245 k = ONES_COMPLEMENT (k) ;
246 ASSERT (k >= 0 && k < n_col) ;
247 Wp [k] = pivcol ;
248 /* save a copy of the inverse column permutation in Fcpos */
249 Fcpos [pivcol] = k ;
250 }
251 for (k = 0 ; k < n_col ; k++)
252 {
253 Cperm [k] = Wp [k] ;
254 }
255
256 #ifndef NDEBUG
257 for (k = 0 ; k < n_col ; k++)
258 {
259 col = Cperm [k] ;
260 ASSERT (col >= 0 && col < n_col) ;
261 ASSERT (Fcpos [col] == k) ; /* col is the kth pivot */
262 }
263 for (k = 0 ; k < n_row ; k++)
264 {
265 row = Rperm [k] ;
266 ASSERT (row >= 0 && row < n_row) ;
267 ASSERT (Frpos [row] == k) ; /* row is the kth pivot */
268 }
269 #endif
270
271 #ifndef NDEBUG
272 UMF_dump_lu (Numeric) ;
273 #endif
274
275 /* ---------------------------------------------------------------------- */
276 /* permute Lpos, Upos, Lilen, Lip, Uilen, and Uip */
277 /* ---------------------------------------------------------------------- */
278
279 for (k = 0 ; k < npiv ; k++)
280 {
281 pivrow = Rperm [k] ;
282 Wr [k] = Uilen [pivrow] ;
283 Wp [k] = Uip [pivrow] ;
284 }
285
286 for (k = 0 ; k < npiv ; k++)
287 {
288 Uilen [k] = Wr [k] ;
289 Uip [k] = Wp [k] ;
290 }
291
292 for (k = 0 ; k < npiv ; k++)
293 {
294 pivrow = Rperm [k] ;
295 Wp [k] = Lpos [pivrow] ;
296 }
297
298 for (k = 0 ; k < npiv ; k++)
299 {
300 Lpos [k] = Wp [k] ;
301 }
302
303 for (k = 0 ; k < npiv ; k++)
304 {
305 pivcol = Cperm [k] ;
306 Wc [k] = Lilen [pivcol] ;
307 Wp [k] = Lip [pivcol] ;
308 }
309
310 for (k = 0 ; k < npiv ; k++)
311 {
312 Lilen [k] = Wc [k] ;
313 Lip [k] = Wp [k] ;
314 }
315
316 for (k = 0 ; k < npiv ; k++)
317 {
318 pivcol = Cperm [k] ;
319 Wp [k] = Upos [pivcol] ;
320 }
321
322 for (k = 0 ; k < npiv ; k++)
323 {
324 Upos [k] = Wp [k] ;
325 }
326
327 /* ---------------------------------------------------------------------- */
328 /* terminate the last Uchain and last Lchain */
329 /* ---------------------------------------------------------------------- */
330
331 Upos [npiv] = EMPTY ;
332 Lpos [npiv] = EMPTY ;
333 Uip [npiv] = EMPTY ;
334 Lip [npiv] = EMPTY ;
335 Uilen [npiv] = 0 ;
336 Lilen [npiv] = 0 ;
337
338 /* ---------------------------------------------------------------------- */
339 /* convert U to the new pivot order */
340 /* ---------------------------------------------------------------------- */
341
342 n1 = Symbolic->n1 ;
343
344 for (k = 0 ; k < n1 ; k++)
345 {
346 /* this is a singleton row of U */
347 ulen = Uilen [k] ;
348 DEBUG4 (("K "ID" New U. ulen "ID" Singleton 1\n", k, ulen)) ;
349 if (ulen > 0)
350 {
351 up = Uip [k] ;
352 ip = (Int *) (Numeric->Memory + up) ;
353 for (i = 0 ; i < ulen ; i++)
354 {
355 col = *ip ;
356 DEBUG4 ((" old col "ID" new col "ID"\n", col, Fcpos [col]));
357 ASSERT (col >= 0 && col < n_col) ;
358 *ip++ = Fcpos [col] ;
359 }
360 }
361 }
362
363 for (k = n1 ; k < npiv ; k++)
364 {
365 up = Uip [k] ;
366 if (up < 0)
367 {
368 /* this is the start of a new Uchain (with a pattern) */
369 ulen = Uilen [k] ;
370 DEBUG4 (("K "ID" New U. ulen "ID" End_Uchain 1\n", k, ulen)) ;
371 if (ulen > 0)
372 {
373 up = -up ;
374 ip = (Int *) (Numeric->Memory + up) ;
375 for (i = 0 ; i < ulen ; i++)
376 {
377 col = *ip ;
378 DEBUG4 ((" old col "ID" new col "ID"\n", col, Fcpos [col]));
379 ASSERT (col >= 0 && col < n_col) ;
380 *ip++ = Fcpos [col] ;
381 }
382 }
383 }
384 }
385
386 ulen = Numeric->ulen ;
387 if (ulen > 0)
388 {
389 /* convert last pivot row of U to the new pivot order */
390 DEBUG4 (("K "ID" (last)\n", k)) ;
391 for (i = 0 ; i < ulen ; i++)
392 {
393 col = Numeric->Upattern [i] ;
394 DEBUG4 ((" old col "ID" new col "ID"\n", col, Fcpos [col])) ;
395 Numeric->Upattern [i] = Fcpos [col] ;
396 }
397 }
398
399 /* Fcpos no longer needed ] */
400
401 /* ---------------------------------------------------------------------- */
402 /* convert L to the new pivot order */
403 /* ---------------------------------------------------------------------- */
404
405 for (k = 0 ; k < n1 ; k++)
406 {
407 llen = Lilen [k] ;
408 DEBUG4 (("K "ID" New L. llen "ID" Singleton col\n", k, llen)) ;
409 if (llen > 0)
410 {
411 lp = Lip [k] ;
412 ip = (Int *) (Numeric->Memory + lp) ;
413 for (i = 0 ; i < llen ; i++)
414 {
415 row = *ip ;
416 DEBUG4 ((" old row "ID" new row "ID"\n", row, Frpos [row])) ;
417 ASSERT (row >= 0 && row < n_row) ;
418 *ip++ = Frpos [row] ;
419 }
420 }
421 }
422
423 for (k = n1 ; k < npiv ; k++)
424 {
425 llen = Lilen [k] ;
426 DEBUG4 (("K "ID" New L. llen "ID" \n", k, llen)) ;
427 if (llen > 0)
428 {
429 lp = Lip [k] ;
430 if (lp < 0)
431 {
432 /* this starts a new Lchain */
433 lp = -lp ;
434 }
435 ip = (Int *) (Numeric->Memory + lp) ;
436 for (i = 0 ; i < llen ; i++)
437 {
438 row = *ip ;
439 DEBUG4 ((" old row "ID" new row "ID"\n", row, Frpos [row])) ;
440 ASSERT (row >= 0 && row < n_row) ;
441 *ip++ = Frpos [row] ;
442 }
443 }
444 }
445
446 /* Frpos no longer needed ] */
447
448 /* ---------------------------------------------------------------------- */
449 /* combine symbolic and numeric permutations */
450 /* ---------------------------------------------------------------------- */
451
452 Cperm_init = Symbolic->Cperm_init ;
453 Rperm_init = Symbolic->Rperm_init ;
454
455 for (k = 0 ; k < n_row ; k++)
456 {
457 Rperm [k] = Rperm_init [Rperm [k]] ;
458 }
459
460 for (k = 0 ; k < n_col ; k++)
461 {
462 Cperm [k] = Cperm_init [Cperm [k]] ;
463 }
464
465 /* Work object will be freed immediately upon return (to UMF_kernel */
466 /* and then to UMFPACK_numeric). */
467 }
468