1 /* ========================================================================== */
2 /* === UMF_kernel_init ====================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis. CISE Dept, */
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 /*
12 Initialize the kernel: scale the matrix, load the initial elements, and
13 build the tuple lists.
14
15 Returns TRUE if successful, FALSE if out of memory or if the pattern has
16 changed since UMFPACK_*symbolic. UMFPACK_numeric allocates at least enough
17 space for UMF_kernel_init to succeed; otherwise it does not call
18 UMF_kernel_init. So an out-of-memory condition means that the pattern must
19 have gotten larger.
20 */
21
22 #include "umf_internal.h"
23 #include "umf_tuple_lengths.h"
24 #include "umf_build_tuples.h"
25 #include "umf_mem_init_memoryspace.h"
26 #include "umf_mem_alloc_element.h"
27 #include "umf_mem_alloc_head_block.h"
28 #include "umf_mem_alloc_tail_block.h"
29 #include "umf_mem_free_tail_block.h"
30 #include "umf_scale.h"
31
32 /* ========================================================================== */
33 /* === packsp =============================================================== */
34 /* ========================================================================== */
35
36 /* remove zero or small entries from a column of L or a row of U */
37
packsp(Int pnew,Int * p_p,Int * p_len,Int drop,double droptol,Unit * Memory)38 PRIVATE Int packsp /* returns new value of pnew */
39 (
40 Int pnew, /* index into Memory of next free space */
41 Int *p_p, /* ptr to index of old pattern in Memory on input,
42 new pattern on output */
43 Int *p_len, /* ptr to length of old pattern on input,
44 new pattern on output */
45 Int drop, /* TRUE if small nonzero entries are to be dropped */
46 double droptol, /* the drop tolerance */
47 Unit *Memory /* contains the sparse vector on input and output */
48 )
49 {
50 Entry x ;
51 double s ;
52 Entry *Bx, *Bx2 ;
53 Int p, i, len, len_new, *Bi, *Bi2 ;
54
55 /* get the pointers to the sparse vector, and its length */
56 p = *p_p ;
57 len = *p_len ;
58 Bi = (Int *) (Memory + p) ; p += UNITS (Int, len) ;
59 Bx = (Entry *) (Memory + p) ; p += UNITS (Entry, len) ;
60 DEBUGm4 ((" p "ID" len "ID" pnew "ID"\n", p, len, pnew)) ;
61
62 /* the vector resides in Bi [0..len-1] and Bx [0..len-1] */
63
64 /* first, compact the vector in place */
65 len_new = 0 ;
66 for (p = 0 ; p < len ; p++)
67 {
68 i = Bi [p] ;
69 x = Bx [p] ;
70 DEBUGm4 ((" old vector: i "ID" value: ", i)) ;
71 EDEBUGk (-4, x) ;
72 DEBUGm4 (("\n")) ;
73 ASSERT (i >= 0) ;
74 /* skip if zero or below drop tolerance */
75 if (IS_ZERO (x)) continue ;
76 if (drop)
77 {
78 APPROX_ABS (s, x) ;
79 if (s <= droptol) continue ;
80 }
81 /* store the value back into the vector */
82 if (len_new != p)
83 {
84 Bi [len_new] = i ;
85 Bx [len_new] = x ;
86 }
87 len_new++ ;
88 }
89 ASSERT (len_new <= len) ;
90
91 /* the vector is now in Bi [0..len_new-1] and Bx [0..len_new-1] */
92
93 #ifndef NDEBUG
94 for (p = 0 ; p < len_new ; p++)
95 {
96 DEBUGm4 ((" new vector: i "ID" value: ", Bi [p])) ;
97 EDEBUGk (-4, Bx [p]) ;
98 DEBUGm4 (("\n")) ;
99 ASSERT (Bi [p] >= 0) ;
100 }
101 #endif
102
103 /* allocate new space for the compacted vector */
104 *p_p = pnew ;
105 *p_len = len_new ;
106 Bi2 = (Int *) (Memory + pnew) ; pnew += UNITS (Int, len_new) ;
107 Bx2 = (Entry *) (Memory + pnew) ; pnew += UNITS (Entry, len_new) ;
108 DEBUGm4 ((" pnew "ID" len_new "ID"\n", pnew, len_new)) ;
109
110 /* shift the vector upwards, into its new space */
111 for (p = 0 ; p < len_new ; p++)
112 {
113 Bi2 [p] = Bi [p] ;
114 }
115 for (p = 0 ; p < len_new ; p++)
116 {
117 Bx2 [p] = Bx [p] ;
118 }
119
120 #ifndef NDEBUG
121 for (p = 0 ; p < len_new ; p++)
122 {
123 DEBUGm4 ((" packed vec: i "ID" value: ", Bi2 [p])) ;
124 EDEBUGk (-4, Bx2 [p]) ;
125 DEBUGm4 (("\n")) ;
126 ASSERT (Bi2 [p] >= 0) ;
127 }
128 #endif
129
130 /* return the pointer to the space just after the new vector */
131 return (pnew) ;
132 }
133
134
135 /* ========================================================================== */
136 /* === UMF_kernel_init ====================================================== */
137 /* ========================================================================== */
138
UMF_kernel_init(const Int Ap[],const Int Ai[],const double Ax[],const double Az[],NumericType * Numeric,WorkType * Work,SymbolicType * Symbolic)139 GLOBAL Int UMF_kernel_init
140 (
141 const Int Ap [ ], /* user's input matrix (not modified) */
142 const Int Ai [ ],
143 const double Ax [ ],
144 #ifdef COMPLEX
145 const double Az [ ],
146 #endif
147 NumericType *Numeric,
148 WorkType *Work,
149 SymbolicType *Symbolic
150 )
151 {
152 /* ---------------------------------------------------------------------- */
153 /* local variables */
154 /* ---------------------------------------------------------------------- */
155
156 Entry x, pivot_value ;
157 double unused = 0, rsmin, rsmax, rs, droptol ;
158 Entry *D, *C, *Lval, **Rpx ;
159 double *Rs ;
160 Int row, k, oldcol, size, e, p1, p2, p, nz, *Rows, *Cols, *E, i, *Upos,
161 *Lpos, n_row, n_col, *Wp, *Cperm_init, *Frpos, *Fcpos, *Row_degree, nn,
162 *Row_tlen, *Col_degree, *Col_tlen, oldrow, newrow, ilast, *Wrp,
163 *Rperm_init, col, n_inner, prefer_diagonal, *Diagonal_map, nempty,
164 *Diagonal_imap, fixQ, rdeg, cdeg, nempty_col, *Esize, esize, pnew,
165 *Lip, *Uip, *Lilen, *Uilen, llen, pa, *Cdeg, *Rdeg, n1, clen, do_scale,
166 lnz, unz, lip, uip, k1, *Rperm, *Cperm, pivcol, *Li, lilen, drop,
167 **Rpi, nempty_row, dense_row_threshold, empty_elements, rpi, rpx ;
168 Element *ep ;
169 Unit *Memory ;
170 #ifdef COMPLEX
171 Int split = SPLIT (Az) ;
172 #endif
173 #ifndef NRECIPROCAL
174 Int do_recip = FALSE ;
175 #endif
176
177 /* ---------------------------------------------------------------------- */
178 /* get parameters */
179 /* ---------------------------------------------------------------------- */
180
181 DEBUG0 (("KERNEL INIT\n")) ;
182
183 n_row = Symbolic->n_row ;
184 n_col = Symbolic->n_col ;
185 nn = MAX (n_row, n_col) ;
186 n_inner = MIN (n_row, n_col) ;
187 nempty_col = Symbolic->nempty_col ;
188 nempty_row = Symbolic->nempty_row ;
189 nempty = MIN (nempty_row, nempty_col) ;
190 ASSERT (n_row > 0 && n_col > 0) ;
191 Cperm_init = Symbolic->Cperm_init ;
192 Rperm_init = Symbolic->Rperm_init ;
193 Cdeg = Symbolic->Cdeg ;
194 Rdeg = Symbolic->Rdeg ;
195 n1 = Symbolic->n1 ;
196 dense_row_threshold = Symbolic->dense_row_threshold ;
197 DEBUG0 (("Singletons: "ID"\n", n1)) ;
198 Work->nforced = 0 ;
199 Work->ndiscard = 0 ;
200 Work->noff_diagonal = 0 ;
201
202 nz = Ap [n_col] ;
203 if (nz < 0 || Ap [0] != 0 || nz != Symbolic->nz)
204 {
205 DEBUGm4 (("nz or Ap [0] bad\n")) ;
206 return (FALSE) ; /* pattern changed */
207 }
208
209 prefer_diagonal = Symbolic->prefer_diagonal ;
210 Diagonal_map = Work->Diagonal_map ;
211 Diagonal_imap = Work->Diagonal_imap ;
212
213 /* ---------------------------------------------------------------------- */
214 /* initialize the Numeric->Memory space for LU, elements, and tuples */
215 /* ---------------------------------------------------------------------- */
216
217 UMF_mem_init_memoryspace (Numeric) ;
218 DEBUG1 (("Kernel init head usage, before allocs: "ID"\n", Numeric->ihead)) ;
219
220 /* ---------------------------------------------------------------------- */
221 /* initialize the Work and Numeric objects */
222 /* ---------------------------------------------------------------------- */
223
224 /* current front is empty */
225 Work->fnpiv = 0 ;
226 Work->fncols = 0 ;
227 Work->fnrows = 0 ;
228 Work->fncols_max = 0 ;
229 Work->fnrows_max = 0 ;
230 Work->fnzeros = 0 ;
231 Work->fcurr_size = 0 ;
232 Work->fnr_curr = 0 ;
233 Work->fnc_curr = 0 ;
234
235 Work->nz = nz ;
236 Work->prior_element = EMPTY ;
237 Work->ulen = 0 ;
238 Work->llen = 0 ;
239 Work->npiv = n1 ;
240 Work->frontid = 0 ;
241 Work->nextcand = n1 ;
242
243 Memory = Numeric->Memory ;
244 Rperm = Numeric->Rperm ;
245 Cperm = Numeric->Cperm ;
246 Row_degree = Numeric->Rperm ;
247 Col_degree = Numeric->Cperm ;
248 /* Row_tuples = Numeric->Uip ; not needed */
249 Row_tlen = Numeric->Uilen ;
250 /* Col_tuples = Numeric->Lip ; not needed */
251 Col_tlen = Numeric->Lilen ;
252
253 Lip = Numeric->Lip ;
254 Uip = Numeric->Uip ;
255 Lilen = Numeric->Lilen ;
256 Uilen = Numeric->Uilen ;
257
258 Frpos = Work->Frpos ;
259 Fcpos = Work->Fcpos ;
260 Wp = Work->Wp ;
261 Wrp = Work->Wrp ;
262
263 D = Numeric->D ;
264 Upos = Numeric->Upos ;
265 Lpos = Numeric->Lpos ;
266 for (k = 0 ; k < n_inner ; k++)
267 {
268 CLEAR (D [k]) ;
269 }
270
271 Rs = Numeric->Rs ;
272
273 for (row = 0 ; row <= n_row ; row++)
274 {
275 Lpos [row] = EMPTY ;
276 /* Row_tuples [row] = 0 ; set in UMF_build_tuples */
277 /* Row_degree [row] = 0 ; initialized below */
278 Row_tlen [row] = 0 ;
279 /* Frpos [row] = EMPTY ; do this later */
280 }
281
282 for (col = 0 ; col <= n_col ; col++)
283 {
284 Upos [col] = EMPTY ;
285 /* Col_tuples [col] = 0 ; set in UMF_build_tuples */
286 /* Col_degree [col] = 0 ; initialized below */
287 Col_tlen [col] = 0 ;
288 Fcpos [col] = EMPTY ;
289 Wrp [col] = 0 ;
290 }
291 Work->Wrpflag = 1 ;
292
293 /* When cleared, Wp [0..nn] is < 0 */
294 for (i = 0 ; i <= nn ; i++)
295 {
296 Wp [i] = EMPTY ;
297 }
298 /* In col search, Wp [row] is set to a position, which is >= 0. */
299
300 /* When cleared, Wrp [0..n_col] is < Wrpflag */
301 /* In row search, Wrp [col] is set to Wrpflag. */
302
303 /* no need to initialize Wm, Wio, Woi, and Woo */
304
305 /* clear the external degree counters */
306 Work->cdeg0 = 1 ;
307 Work->rdeg0 = 1 ;
308
309 fixQ = Symbolic->fixQ ;
310
311 E = Work->E ;
312
313 Numeric->n_row = n_row ;
314 Numeric->n_col = n_col ;
315 Numeric->npiv = 0 ;
316 Numeric->nnzpiv = 0 ;
317 Numeric->min_udiag = 0.0 ;
318 Numeric->max_udiag = 0.0 ;
319 Numeric->rcond = 0.0 ;
320 Numeric->isize = 0 ;
321 Numeric->nLentries = 0 ;
322 Numeric->nUentries = 0 ;
323 Numeric->lnz = 0 ;
324 Numeric->unz = 0 ;
325 Numeric->all_lnz = 0 ;
326 Numeric->all_unz = 0 ;
327 Numeric->maxfrsize = 0 ;
328 Numeric->maxnrows = 0 ;
329 Numeric->maxncols = 0 ;
330 Numeric->flops = 0. ;
331 Numeric->n1 = n1 ;
332 droptol = Numeric->droptol ;
333 drop = (droptol > 0) ;
334
335 /* ---------------------------------------------------------------------- */
336 /* compute the scale factors, if requested, and check the input matrix */
337 /* ---------------------------------------------------------------------- */
338
339 /* UMFPACK_SCALE_SUM: Rs [i] = sum of the absolute values in row i.
340 * UMFPACK_SCALE_MAX: Rs [i] = max of the absolute values in row i.
341 *
342 * If A is complex, an approximate abs is used (|xreal| + |ximag|).
343 *
344 * If min (Rs [0..n_row]) >= RECIPROCAL_TOLERANCE, then the scale
345 * factors are inverted, and the rows of A are multiplied by the scale
346 * factors. Otherwise, the rows are divided by the scale factors. If
347 * NRECIPROCAL is defined, then the rows are always divided by the scale
348 * factors.
349 *
350 * For MATLAB (either built-in routine or mexFunction), or for gcc,
351 * the rows are always divided by the scale factors.
352 */
353
354 do_scale = (Numeric->scale != UMFPACK_SCALE_NONE) ;
355
356 if (do_scale)
357 {
358 int do_max = Numeric->scale == UMFPACK_SCALE_MAX ;
359 for (row = 0 ; row < n_row ; row++)
360 {
361 Rs [row] = 0.0 ;
362 }
363 for (col = 0 ; col < n_col ; col++)
364 {
365 ilast = EMPTY ;
366 p1 = Ap [col] ;
367 p2 = Ap [col+1] ;
368 if (p1 > p2)
369 {
370 /* invalid matrix */
371 DEBUGm4 (("invalid matrix (Ap)\n")) ;
372 return (FALSE) ;
373 }
374 for (p = p1 ; p < p2 ; p++)
375 {
376 Entry aij ;
377 double value ;
378 row = Ai [p] ;
379 if (row <= ilast || row >= n_row)
380 {
381 /* invalid matrix, columns must be sorted, no duplicates */
382 DEBUGm4 (("invalid matrix (Ai)\n")) ;
383 return (FALSE) ;
384 }
385 ASSIGN (aij, Ax, Az, p, split) ;
386 APPROX_ABS (value, aij) ;
387 rs = Rs [row] ;
388 if (!SCALAR_IS_NAN (rs))
389 {
390 if (SCALAR_IS_NAN (value))
391 {
392 /* if any entry in the row is NaN, then the scale factor
393 * is NaN too (for now) and then set to 1.0 below */
394 Rs [row] = value ;
395 }
396 else if (do_max)
397 {
398 Rs [row] = MAX (rs, value) ;
399 }
400 else
401 {
402 Rs [row] += value ;
403 }
404 }
405 DEBUG4 (("i "ID" j "ID" value %g, Rs[i]: %g\n",
406 row, col, value, Rs[row])) ;
407 ilast = row ;
408 }
409 }
410 DEBUG2 (("Rs[0] = %30.20e\n", Rs [0])) ;
411 for (row = 0 ; row < n_row ; row++)
412 {
413 rs = Rs [row] ;
414 if (SCALAR_IS_ZERO (rs) || SCALAR_IS_NAN (rs))
415 {
416 /* don't scale a completely zero row, or one with NaN's */
417 Rs [row] = 1.0 ;
418 }
419 }
420 rsmin = Rs [0] ;
421 rsmax = Rs [0] ;
422 for (row = 0 ; row < n_row ; row++)
423 {
424 DEBUG2 (("sum %30.20e ", Rs [row])) ;
425 rsmin = MIN (rsmin, Rs [row]) ;
426 rsmax = MAX (rsmax, Rs [row]) ;
427 DEBUG2 (("Rs["ID"] = %30.20e\n", row, Rs [row])) ;
428 }
429 #ifndef NRECIPROCAL
430 /* multiply by the reciprocal if Rs is not too small */
431 do_recip = (rsmin >= RECIPROCAL_TOLERANCE) ;
432 if (do_recip)
433 {
434 /* invert the scale factors */
435 for (row = 0 ; row < n_row ; row++)
436 {
437 Rs [row] = 1.0 / Rs [row] ;
438 }
439 }
440 #endif
441 }
442 else
443 {
444 /* no scaling, rsmin and rsmax not computed */
445 rsmin = -1 ;
446 rsmax = -1 ;
447 #ifndef NRECIPROCAL
448 do_recip = FALSE ;
449 #endif
450 /* check the input matrix */
451 if (!AMD_valid (n_row, n_col, Ap, Ai))
452 {
453 /* matrix is invalid */
454 return (FALSE) ;
455 }
456 }
457
458 Numeric->rsmin = rsmin ;
459 Numeric->rsmax = rsmax ;
460 #ifndef NRECIPROCAL
461 Numeric->do_recip = do_recip ;
462 #else
463 Numeric->do_recip = FALSE ;
464 #endif
465
466 /* ---------------------------------------------------------------------- */
467 /* construct the inverse row Rperm_init permutation (use Frpos as temp) */
468 /* ---------------------------------------------------------------------- */
469
470 DEBUG3 (("\n\n===================LOAD_MATRIX:\n")) ;
471
472 for (newrow = 0 ; newrow < n_row ; newrow++)
473 {
474 oldrow = Rperm_init [newrow] ;
475 ASSERT (oldrow >= 0 && oldrow < n_row) ;
476 Frpos [oldrow] = newrow ;
477 }
478
479 /* ---------------------------------------------------------------------- */
480 /* construct the diagonal imap if doing symmetric pivoting */
481 /* ---------------------------------------------------------------------- */
482
483 if (prefer_diagonal)
484 {
485 ASSERT (n_row == n_col) ;
486 ASSERT (nempty_col == Symbolic->nempty_row) ;
487 ASSERT (nempty_col == nempty) ;
488 for (i = 0 ; i < nn ; i++)
489 {
490 Diagonal_map [i] = EMPTY ;
491 Diagonal_imap [i] = EMPTY ;
492 }
493 for (k = n1 ; k < nn - nempty ; k++)
494 {
495 newrow = Symbolic->Diagonal_map [k] ;
496 Diagonal_map [k] = newrow ;
497 Diagonal_imap [newrow] = k ;
498 }
499 }
500
501 /* ---------------------------------------------------------------------- */
502 /* allocate O (n_row) workspace at the tail end of Memory */
503 /* ---------------------------------------------------------------------- */
504
505 rpi = UMF_mem_alloc_tail_block (Numeric, UNITS (Int *, n_row+1)) ;
506 rpx = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry *, n_row+1)) ;
507 if (!rpi || !rpx)
508 {
509 /* :: pattern change (out of memory for Rpx, Rpx) :: */
510 /* out of memory, which can only mean that the pattern has changed */
511 return (FALSE) ; /* pattern changed */
512 }
513 Rpi = (Int **) (Memory + rpx) ;
514 Rpx = (Entry **) (Memory + rpi) ;
515
516 /* ---------------------------------------------------------------------- */
517 /* allocate the LU factors for the columns of the singletons */
518 /* ---------------------------------------------------------------------- */
519
520 DEBUG1 (("Allocating singletons:\n")) ;
521 for (k = 0 ; k < n1 ; k++)
522 {
523 lnz = Cdeg [k] - 1 ;
524 unz = Rdeg [k] - 1 ;
525
526 DEBUG1 (("Singleton k "ID" pivrow "ID" pivcol "ID" cdeg "ID" rdeg "
527 ID"\n", k, Rperm_init [k], Cperm_init [k], Cdeg [k], Rdeg [k])) ;
528 ASSERT (unz >= 0 && lnz >= 0 && (lnz == 0 || unz == 0)) ;
529 DEBUG1 ((" lnz "ID" unz "ID"\n", lnz, unz)) ;
530
531 size = UNITS (Int, lnz) + UNITS (Entry, lnz)
532 + UNITS (Int, unz) + UNITS (Entry, unz) ;
533 p = UMF_mem_alloc_head_block (Numeric, size) ;
534 DEBUG1 (("Kernel init head usage: "ID"\n", Numeric->ihead)) ;
535 if (!p)
536 {
537 /* :: pattern change (out of memory for singletons) :: */
538 DEBUG0 (("Pattern has gotten larger - kernel init failed\n")) ;
539 return (FALSE) ; /* pattern changed */
540 }
541
542 Numeric->all_lnz += lnz ;
543 Numeric->all_unz += unz ;
544
545 /* allocate the column of L */
546 lip = p ;
547 p += UNITS (Int, lnz) ;
548 p += UNITS (Entry, lnz) ;
549
550 /* allocate the row of U */
551 uip = p ;
552 Rpi [k] = (Int *) (Memory + p) ;
553 p += UNITS (Int, unz) ;
554 Rpx [k] = (Entry *) (Memory + p) ;
555 /* p += UNITS (Entry, unz) ; (not needed) */
556
557 /* a single column of L (no Lchains) */
558 Lip [k] = lip ;
559 Lilen [k] = lnz ;
560
561 /* a single row of L (no Uchains) */
562 Uip [k] = uip ;
563 Uilen [k] = unz ;
564
565 Wp [k] = unz ;
566
567 /* save row and column inverse permutation */
568 k1 = ONES_COMPLEMENT (k) ;
569 Rperm [k] = k1 ; /* aliased with Row_degree */
570 Cperm [k] = k1 ; /* aliased with Col_degree */
571 }
572
573 /* ---------------------------------------------------------------------- */
574 /* current frontal matrix is empty */
575 /* ---------------------------------------------------------------------- */
576
577 e = 0 ;
578 E [e] = 0 ;
579 Work->Flublock = (Entry *) NULL ;
580 Work->Flblock = (Entry *) NULL ;
581 Work->Fublock = (Entry *) NULL ;
582 Work->Fcblock = (Entry *) NULL ;
583
584 /* ---------------------------------------------------------------------- */
585 /* allocate the column elements */
586 /* ---------------------------------------------------------------------- */
587
588 Esize = Symbolic->Esize ;
589 empty_elements = FALSE ;
590 for (k = n1 ; k < n_col - nempty_col ; k++)
591 {
592 e = k - n1 + 1 ;
593 ASSERT (e < Work->elen) ;
594 esize = Esize ? Esize [k-n1] : Cdeg [k] ;
595 if (esize > 0)
596 {
597 /* allocate an element for this column */
598 E [e] = UMF_mem_alloc_element (Numeric, esize, 1, &Rows, &Cols, &C,
599 &size, &ep) ;
600 if (E [e] <= 0)
601 {
602 /* :: pattern change (out of memory for column elements) :: */
603 return (FALSE) ; /* pattern has changed */
604 }
605 Cols [0] = k ;
606 DEBUG0 (("Got column element e "ID" esize "ID"\n", e, esize)) ;
607 }
608 else
609 {
610 /* all rows in this column are dense, or empty */
611 E [e] = 0 ;
612 empty_elements = TRUE ;
613 DEBUG0 (("column element e is empty "ID"\n", e)) ;
614 }
615 }
616 DEBUG0 (("e "ID" n_col "ID" nempty_col "ID" n1 "ID"\n", e, n_col,
617 nempty_col, n1)) ;
618 ASSERT (e == n_col - nempty_col - n1) ;
619
620 /* ---------------------------------------------------------------------- */
621 /* allocate the row elements for dense rows of A (if any) */
622 /* ---------------------------------------------------------------------- */
623
624 if (Esize)
625 {
626 for (k = n1 ; k < n_row - nempty_row ; k++)
627 {
628 rdeg = Rdeg [k] ;
629 if (rdeg > dense_row_threshold)
630 {
631 /* allocate an element for this dense row */
632 e++ ;
633 ASSERT (e < Work->elen) ;
634 E [e] = UMF_mem_alloc_element (Numeric, 1, rdeg, &Rows, &Cols,
635 &C, &size, &ep) ;
636 if (E [e] <= 0)
637 {
638 /* :: pattern change (out of memory for row elements) :: */
639 return (FALSE) ; /* pattern has changed */
640 }
641 Rows [0] = k ;
642 Rpi [k] = Cols ;
643 Rpx [k] = C ;
644 Wp [k] = rdeg ;
645 DEBUG0 (("Got row element e "ID" rdeg "ID"\n", e, rdeg)) ;
646 }
647 }
648 }
649
650 /* elements are currently in the range 0 to e */
651 Work->nel = e ;
652
653 /* ---------------------------------------------------------------------- */
654 /* create the first n1 columns of L and U */
655 /* ---------------------------------------------------------------------- */
656
657 for (k = 0 ; k < n1 ; k++)
658 {
659 pivcol = Cperm_init [k] ;
660 p2 = Ap [pivcol+1] ;
661
662 /* get the kth column of L */
663 p = Lip [k] ;
664 Li = (Int *) (Memory + p) ;
665 lilen = Lilen [k] ;
666 p += UNITS (Int, lilen) ;
667 Lval = (Entry *) (Memory + p) ;
668
669 llen = 0 ;
670
671 for (pa = Ap [pivcol] ; pa < p2 ; pa++)
672 {
673 oldrow = Ai [pa] ;
674 newrow = Frpos [oldrow] ;
675 ASSIGN (x, Ax, Az, pa, split) ;
676
677 /* scale the value using the scale factors, Rs */
678 if (do_scale)
679 {
680 #ifndef NRECIPROCAL
681 if (do_recip)
682 {
683 SCALE (x, Rs [oldrow]) ;
684 }
685 else
686 #endif
687 {
688 SCALE_DIV (x, Rs [oldrow]) ;
689 }
690 }
691
692 if (newrow == k)
693 {
694 /* this is the pivot entry itself */
695 ASSERT (oldrow == Rperm_init [k]) ;
696 D [k] = x ;
697 }
698 else if (newrow < k)
699 {
700 /* this entry goes in a row of U */
701 DEBUG1 (("Singleton row of U: k "ID" newrow "ID"\n",
702 k, newrow)) ;
703 if (--(Wp [newrow]) < 0)
704 {
705 /* :: pattern change (singleton row too long) :: */
706 DEBUGm4 (("bad U singleton row (too long)\n")) ;
707 return (FALSE) ; /* pattern changed */
708 }
709 *(Rpi [newrow]++) = k ;
710 *(Rpx [newrow]++) = x ;
711 }
712 else
713 {
714 /* this entry goes in a column of L */
715 DEBUG1 (("Singleton col of L: k "ID" newrow "ID"\n",
716 k, newrow)) ;
717 if (llen >= lilen)
718 {
719 DEBUGm4 (("bad L singleton col (too long)\n")) ;
720 return (FALSE) ; /* pattern changed */
721 }
722 Li [llen] = newrow ;
723 Lval [llen] = x ;
724 llen++ ;
725 }
726 }
727
728 if (llen != lilen)
729 {
730 /* :: pattern change (singleton column too long) :: */
731 DEBUGm4 (("bad L singleton col (too short)\n")) ;
732 return (FALSE) ; /* pattern changed */
733 }
734
735 /* scale the column of L */
736 if (llen > 0)
737 {
738 pivot_value = D [k] ;
739 UMF_scale (llen, pivot_value, Lval) ;
740 }
741
742 }
743
744 /* ---------------------------------------------------------------------- */
745 /* allocate the elements and copy the columns of A */
746 /* ---------------------------------------------------------------------- */
747
748 /* also apply the row and column pre-ordering. */
749 for (k = n1 ; k < n_col ; k++)
750 {
751 /* The newcol is k, which is what the name of the column is in the
752 * UMFPACK kernel. The user's name for the column is oldcol. */
753 oldcol = Cperm_init [k] ;
754
755 ASSERT (oldcol >= 0 && oldcol < n_col) ;
756
757 p2 = Ap [oldcol+1] ;
758
759 cdeg = Cdeg [k] ;
760 ASSERT (cdeg >= 0) ;
761 ASSERT (IMPLIES (
762 (Symbolic->ordering != UMFPACK_ORDERING_GIVEN) && n1 > 0,
763 cdeg > 1 || cdeg == 0)) ;
764
765 /* if fixQ: set Col_degree to 0 for the NON_PIVOTAL_COL macro */
766 Col_degree [k] = fixQ ? 0 : cdeg ;
767
768 /* get the element for this column (if any) */
769 e = k - n1 + 1 ;
770 if (k < n_col - nempty_col)
771 {
772 esize = Esize ? Esize [k-n1] : cdeg ;
773 if (E [e])
774 {
775 Int ncols, nrows ;
776 Unit *pp ;
777 pp = Memory + E [e] ;
778 GET_ELEMENT (ep, pp, Cols, Rows, ncols, nrows, C) ;
779 ASSERT (ncols == 1) ;
780 ASSERT (nrows == esize) ;
781 ASSERT (Cols [0] == k) ;
782 }
783 }
784 else
785 {
786 ASSERT (cdeg == 0) ;
787 esize = 0 ;
788 }
789
790 clen = 0 ;
791
792 for (pa = Ap [oldcol] ; pa < p2 ; pa++)
793 {
794 oldrow = Ai [pa] ;
795 newrow = Frpos [oldrow] ;
796 ASSIGN (x, Ax, Az, pa, split) ;
797
798 /* scale the value using the scale factors, Rs */
799 if (do_scale)
800 {
801 #ifndef NRECIPROCAL
802 if (do_recip)
803 {
804 /* multiply by the reciprocal */
805 SCALE (x, Rs [oldrow]) ;
806 }
807 else
808 #endif
809 {
810 /* divide instead */
811 SCALE_DIV (x, Rs [oldrow]) ;
812 }
813 }
814
815 rdeg = Rdeg [newrow] ;
816 if (newrow < n1 || rdeg > dense_row_threshold)
817 {
818 /* this entry goes in a row of U or into a dense row */
819 DEBUG1 (("Singleton/dense row of U: k "ID" newrow "ID"\n",
820 k, newrow)) ;
821 if (--(Wp [newrow]) < 0)
822 {
823 DEBUGm4 (("bad row of U or A (too long)\n")) ;
824 return (FALSE) ; /* pattern changed */
825 }
826 *(Rpi [newrow]++) = k ;
827 *(Rpx [newrow]++) = x ;
828 }
829 else
830 {
831 /* this entry goes in an initial element */
832 DEBUG1 (("In element k "ID" e "ID" newrow "ID"\n",
833 k, e, newrow)) ;
834 if (clen >= esize)
835 {
836 DEBUGm4 (("bad A column (too long)\n")) ;
837 return (FALSE) ; /* pattern changed */
838 }
839 ASSERT (E [e]) ;
840 ASSERT (k < n_col - nempty_col) ;
841 Rows [clen] = newrow ;
842 C [clen] = x ;
843 clen++ ;
844 #ifndef NDEBUG
845 if (Diagonal_map && (newrow == Diagonal_map [k]))
846 {
847 DEBUG0 (("Diagonal: old: row "ID" col "ID" : "
848 "new: row "ID" col "ID" : ",
849 oldrow, oldcol, newrow, k)) ;
850 EDEBUGk (0, x) ;
851 }
852 #endif
853 }
854 }
855
856 if (clen != esize)
857 {
858 /* :: pattern change (singleton column too short) :: */
859 DEBUGm4 (("bad A column (too short)\n")) ;
860 return (FALSE) ; /* pattern changed */
861 }
862 }
863
864 /* ---------------------------------------------------------------------- */
865 /* free the Rpi and Rpx workspace at the tail end of memory */
866 /* ---------------------------------------------------------------------- */
867
868 UMF_mem_free_tail_block (Numeric, rpi) ;
869 UMF_mem_free_tail_block (Numeric, rpx) ;
870
871 /* ---------------------------------------------------------------------- */
872 /* prune zeros and small entries from the singleton rows and columns */
873 /* ---------------------------------------------------------------------- */
874
875 if (n1 > 0)
876 {
877 pnew = Lip [0] ;
878 ASSERT (pnew == 1) ;
879 for (k = 0 ; k < n1 ; k++)
880 {
881 DEBUGm4 (("\nPrune singleton L col "ID"\n", k)) ;
882 pnew = packsp (pnew, &Lip [k], &Lilen [k], drop, droptol, Memory) ;
883 Numeric->lnz += Lilen [k] ;
884 DEBUGm4 (("\nPrune singleton U row "ID"\n", k)) ;
885 pnew = packsp (pnew, &Uip [k], &Uilen [k], drop, droptol, Memory) ;
886 Numeric->unz += Uilen [k] ;
887 }
888 /* free the unused space at the head of memory */
889 Numeric->ihead = pnew ;
890 }
891
892 /* ---------------------------------------------------------------------- */
893 /* initialize row degrees */
894 /* ---------------------------------------------------------------------- */
895
896 for (k = 0 ; k < n1 ; k++)
897 {
898 if (Wp [k] != 0)
899 {
900 /* :: pattern change (singleton row too short) :: */
901 DEBUGm4 (("bad U singleton row (too short)\n")) ;
902 return (FALSE) ; /* pattern changed */
903 }
904 }
905
906 for (k = n1 ; k < n_row ; k++)
907 {
908 DEBUG1 (("Initial row degree k "ID" oldrow "ID" Rdeg "ID"\n",
909 k, Rperm_init [k], Rdeg [k])) ;
910 rdeg = Rdeg [k] ;
911 Row_degree [k] = rdeg ;
912 if (rdeg > dense_row_threshold && Wp [k] != 0)
913 {
914 /* :: pattern change (dense row too short) :: */
915 DEBUGm4 (("bad dense row (too short)\n")) ;
916 return (FALSE) ; /* pattern changed */
917 }
918 }
919
920 #ifndef NDEBUG
921 if (prefer_diagonal)
922 {
923 Entry aij ;
924 Int *InvCperm, newcol ;
925 UMF_dump_diagonal_map (Diagonal_map, Diagonal_imap, n1, nn, nempty) ;
926 InvCperm = (Int *) malloc (n_col * sizeof (Int)) ;
927 ASSERT (InvCperm != (Int *) NULL) ;
928 for (newcol = 0 ; newcol < n_col ; newcol++)
929 {
930 oldcol = Cperm_init [newcol] ;
931 InvCperm [oldcol] = newcol ;
932 }
933 DEBUGm3 (("Diagonal of P2*A:\n")) ;
934 for (oldcol = 0 ; oldcol < n_col ; oldcol++)
935 {
936 newcol = InvCperm [oldcol] ;
937 for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
938 {
939 oldrow = Ai [p] ;
940 newrow = Frpos [oldrow] ;
941 ASSIGN (aij, Ax, Az, p, split) ;
942 if (newrow == Diagonal_map [newcol])
943 {
944 DEBUG0 (("old row "ID" col "ID" new row "ID" col "ID,
945 oldrow, oldcol, newrow, newcol)) ;
946 EDEBUGk (0, aij) ;
947 DEBUG0 ((" scaled ")) ;
948 if (do_scale)
949 {
950 #ifndef NRECIPROCAL
951 if (do_recip)
952 {
953 SCALE (aij, Rs [oldrow]) ;
954 }
955 else
956 #endif
957 {
958 SCALE_DIV (aij, Rs [oldrow]) ;
959 }
960 }
961 EDEBUGk (0, aij) ;
962 DEBUG0 (("\n")) ;
963 }
964 }
965 }
966 free (InvCperm) ;
967 }
968 #endif
969
970 Col_degree [n_col] = 0 ;
971
972 /* ---------------------------------------------------------------------- */
973 /* pack the element name space */
974 /* ---------------------------------------------------------------------- */
975
976 if (empty_elements)
977 {
978 Int e2 = 0 ;
979 DEBUG0 (("\n\n============= Packing element space\n")) ;
980 for (e = 1 ; e <= Work->nel ; e++)
981 {
982 if (E [e])
983 {
984 e2++ ;
985 E [e2] = E [e] ;
986 }
987 }
988 Work->nel = e2 ;
989 }
990
991 #ifndef NDEBUG
992 DEBUG0 (("Number of initial elements: "ID"\n", Work->nel)) ;
993 for (e = 0 ; e <= Work->nel ; e++) UMF_dump_element (Numeric, Work,e,TRUE) ;
994 #endif
995
996 for (e = Work->nel + 1 ; e < Work->elen ; e++)
997 {
998 E [e] = 0 ;
999 }
1000
1001 /* Frpos no longer needed */
1002 for (row = 0 ; row <= n_row ; row++)
1003 {
1004 Frpos [row] = EMPTY ;
1005 }
1006
1007 /* clear Wp */
1008 for (i = 0 ; i <= nn ; i++)
1009 {
1010 Wp [i] = EMPTY ;
1011 }
1012
1013 DEBUG1 (("Kernel init head usage: "ID"\n", Numeric->ihead)) ;
1014
1015 /* ---------------------------------------------------------------------- */
1016 /* build the tuple lists */
1017 /* ---------------------------------------------------------------------- */
1018
1019 /* if the memory usage changes, then the pattern has changed */
1020
1021 (void) UMF_tuple_lengths (Numeric, Work, &unused) ;
1022 if (!UMF_build_tuples (Numeric, Work))
1023 {
1024 /* :: pattern change (out of memory in umf_build_tuples) :: */
1025 /* We ran out of memory, which can only mean that */
1026 /* the pattern (Ap and or Ai) has changed (gotten larger). */
1027 DEBUG0 (("Pattern has gotten larger - build tuples failed\n")) ;
1028 return (FALSE) ; /* pattern changed */
1029 }
1030
1031 Numeric->init_usage = Numeric->max_usage ;
1032
1033 /* ---------------------------------------------------------------------- */
1034 /* construct the row merge sets */
1035 /* ---------------------------------------------------------------------- */
1036
1037 for (i = 0 ; i <= Symbolic->nfr ; i++)
1038 {
1039 Work->Front_new1strow [i] = Symbolic->Front_1strow [i] ;
1040 }
1041
1042 #ifndef NDEBUG
1043 UMF_dump_rowmerge (Numeric, Symbolic, Work) ;
1044 DEBUG6 (("Column form of original matrix:\n")) ;
1045 UMF_dump_col_matrix (Ax,
1046 #ifdef COMPLEX
1047 Az,
1048 #endif
1049 Ai, Ap, n_row, n_col, nz) ;
1050 UMF_dump_memory (Numeric) ;
1051 UMF_dump_matrix (Numeric, Work, FALSE) ;
1052 DEBUG0 (("kernel init done...\n")) ;
1053 #endif
1054
1055 return (TRUE) ;
1056 }
1057