1 /* ========================================================================== */
2 /* === UMF_kernel_init ====================================================== */
3 /* ========================================================================== */
4
5 /* -------------------------------------------------------------------------- */
6 /* Copyright (c) 2005-2012 by Timothy A. Davis, http://www.suitesparse.com. */
7 /* All Rights Reserved. See ../Doc/License.txt for License. */
8 /* -------------------------------------------------------------------------- */
9
10 /*
11 Initialize the kernel: scale the matrix, load the initial elements, and
12 build the tuple lists.
13
14 Returns TRUE if successful, FALSE if out of memory or if the pattern has
15 changed since UMFPACK_*symbolic. UMFPACK_numeric allocates at least enough
16 space for UMF_kernel_init to succeed; otherwise it does not call
17 UMF_kernel_init. So an out-of-memory condition means that the pattern must
18 have gotten larger.
19 */
20
21 #include "umf_internal.h"
22 #include "umf_kernel_init.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) != AMD_OK)
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
489 #ifndef NDEBUG
490 for (i = 0 ; i < nn ; i++)
491 {
492 Diagonal_map [i] = EMPTY ;
493 Diagonal_imap [i] = EMPTY ;
494 }
495 #endif
496
497 for (k = 0 ; k < nn ; k++)
498 {
499 newrow = Symbolic->Diagonal_map [k] ;
500 Diagonal_map [k] = newrow ;
501 Diagonal_imap [newrow] = k ;
502 }
503
504 #ifndef NDEBUG
505 for (i = 0 ; i < nn ; i++)
506 {
507 ASSERT (Diagonal_map [i] != EMPTY) ;
508 ASSERT (Diagonal_imap [i] != EMPTY) ;
509 }
510 #endif
511 }
512
513 /* ---------------------------------------------------------------------- */
514 /* allocate O (n_row) workspace at the tail end of Memory */
515 /* ---------------------------------------------------------------------- */
516
517 rpi = UMF_mem_alloc_tail_block (Numeric, UNITS (Int *, n_row+1)) ;
518 rpx = UMF_mem_alloc_tail_block (Numeric, UNITS (Entry *, n_row+1)) ;
519 if (!rpi || !rpx)
520 {
521 /* :: pattern change (out of memory for Rpx, Rpx) :: */
522 /* out of memory, which can only mean that the pattern has changed */
523 return (FALSE) ; /* pattern changed */
524 }
525 Rpi = (Int **) (Memory + rpx) ;
526 Rpx = (Entry **) (Memory + rpi) ;
527
528 /* ---------------------------------------------------------------------- */
529 /* allocate the LU factors for the columns of the singletons */
530 /* ---------------------------------------------------------------------- */
531
532 DEBUG1 (("Allocating singletons:\n")) ;
533 for (k = 0 ; k < n1 ; k++)
534 {
535 lnz = Cdeg [k] - 1 ;
536 unz = Rdeg [k] - 1 ;
537
538 DEBUG1 (("Singleton k "ID" pivrow "ID" pivcol "ID" cdeg "ID" rdeg "
539 ID"\n", k, Rperm_init [k], Cperm_init [k], Cdeg [k], Rdeg [k])) ;
540 ASSERT (unz >= 0 && lnz >= 0 && (lnz == 0 || unz == 0)) ;
541 DEBUG1 ((" lnz "ID" unz "ID"\n", lnz, unz)) ;
542
543 size = UNITS (Int, lnz) + UNITS (Entry, lnz)
544 + UNITS (Int, unz) + UNITS (Entry, unz) ;
545 p = UMF_mem_alloc_head_block (Numeric, size) ;
546 DEBUG1 (("Kernel init head usage: "ID"\n", Numeric->ihead)) ;
547 if (!p)
548 {
549 /* :: pattern change (out of memory for singletons) :: */
550 DEBUG0 (("Pattern has gotten larger - kernel init failed\n")) ;
551 return (FALSE) ; /* pattern changed */
552 }
553
554 Numeric->all_lnz += lnz ;
555 Numeric->all_unz += unz ;
556
557 /* allocate the column of L */
558 lip = p ;
559 p += UNITS (Int, lnz) ;
560 p += UNITS (Entry, lnz) ;
561
562 /* allocate the row of U */
563 uip = p ;
564 Rpi [k] = (Int *) (Memory + p) ;
565 p += UNITS (Int, unz) ;
566 Rpx [k] = (Entry *) (Memory + p) ;
567 /* p += UNITS (Entry, unz) ; (not needed) */
568
569 /* a single column of L (no Lchains) */
570 Lip [k] = lip ;
571 Lilen [k] = lnz ;
572
573 /* a single row of L (no Uchains) */
574 Uip [k] = uip ;
575 Uilen [k] = unz ;
576
577 Wp [k] = unz ;
578
579 /* save row and column inverse permutation */
580 k1 = ONES_COMPLEMENT (k) ;
581 Rperm [k] = k1 ; /* aliased with Row_degree */
582 Cperm [k] = k1 ; /* aliased with Col_degree */
583 }
584
585 /* ---------------------------------------------------------------------- */
586 /* current frontal matrix is empty */
587 /* ---------------------------------------------------------------------- */
588
589 e = 0 ;
590 E [e] = 0 ;
591 Work->Flublock = (Entry *) NULL ;
592 Work->Flblock = (Entry *) NULL ;
593 Work->Fublock = (Entry *) NULL ;
594 Work->Fcblock = (Entry *) NULL ;
595
596 /* ---------------------------------------------------------------------- */
597 /* allocate the column elements */
598 /* ---------------------------------------------------------------------- */
599
600 Esize = Symbolic->Esize ;
601 empty_elements = FALSE ;
602 for (k = n1 ; k < n_col - nempty_col ; k++)
603 {
604 e = k - n1 + 1 ;
605 ASSERT (e < Work->elen) ;
606 esize = Esize ? Esize [k-n1] : Cdeg [k] ;
607 if (esize > 0)
608 {
609 /* allocate an element for this column */
610 E [e] = UMF_mem_alloc_element (Numeric, esize, 1, &Rows, &Cols, &C,
611 &size, &ep) ;
612 if (E [e] <= 0)
613 {
614 /* :: pattern change (out of memory for column elements) :: */
615 return (FALSE) ; /* pattern has changed */
616 }
617 Cols [0] = k ;
618 DEBUG0 (("Got column element e "ID" esize "ID"\n", e, esize)) ;
619 }
620 else
621 {
622 /* all rows in this column are dense, or empty */
623 E [e] = 0 ;
624 empty_elements = TRUE ;
625 DEBUG0 (("column element e is empty "ID"\n", e)) ;
626 }
627 }
628 DEBUG0 (("e "ID" n_col "ID" nempty_col "ID" n1 "ID"\n", e, n_col,
629 nempty_col, n1)) ;
630 ASSERT (e == n_col - nempty_col - n1) ;
631
632 /* ---------------------------------------------------------------------- */
633 /* allocate the row elements for dense rows of A (if any) */
634 /* ---------------------------------------------------------------------- */
635
636 if (Esize)
637 {
638 for (k = n1 ; k < n_row - nempty_row ; k++)
639 {
640 rdeg = Rdeg [k] ;
641 if (rdeg > dense_row_threshold)
642 {
643 /* allocate an element for this dense row */
644 e++ ;
645 ASSERT (e < Work->elen) ;
646 E [e] = UMF_mem_alloc_element (Numeric, 1, rdeg, &Rows, &Cols,
647 &C, &size, &ep) ;
648 if (E [e] <= 0)
649 {
650 /* :: pattern change (out of memory for row elements) :: */
651 return (FALSE) ; /* pattern has changed */
652 }
653 Rows [0] = k ;
654 Rpi [k] = Cols ;
655 Rpx [k] = C ;
656 Wp [k] = rdeg ;
657 DEBUG0 (("Got row element e "ID" rdeg "ID"\n", e, rdeg)) ;
658 }
659 }
660 }
661
662 /* elements are currently in the range 0 to e */
663 Work->nel = e ;
664
665 /* ---------------------------------------------------------------------- */
666 /* create the first n1 columns of L and U */
667 /* ---------------------------------------------------------------------- */
668
669 for (k = 0 ; k < n1 ; k++)
670 {
671 pivcol = Cperm_init [k] ;
672 p2 = Ap [pivcol+1] ;
673
674 /* get the kth column of L */
675 p = Lip [k] ;
676 Li = (Int *) (Memory + p) ;
677 lilen = Lilen [k] ;
678 p += UNITS (Int, lilen) ;
679 Lval = (Entry *) (Memory + p) ;
680
681 llen = 0 ;
682
683 for (pa = Ap [pivcol] ; pa < p2 ; pa++)
684 {
685 oldrow = Ai [pa] ;
686 newrow = Frpos [oldrow] ;
687 ASSIGN (x, Ax, Az, pa, split) ;
688
689 /* scale the value using the scale factors, Rs */
690 if (do_scale)
691 {
692 #ifndef NRECIPROCAL
693 if (do_recip)
694 {
695 SCALE (x, Rs [oldrow]) ;
696 }
697 else
698 #endif
699 {
700 SCALE_DIV (x, Rs [oldrow]) ;
701 }
702 }
703
704 if (newrow == k)
705 {
706 /* this is the pivot entry itself */
707 ASSERT (oldrow == Rperm_init [k]) ;
708 D [k] = x ;
709 }
710 else if (newrow < k)
711 {
712 /* this entry goes in a row of U */
713 DEBUG1 (("Singleton row of U: k "ID" newrow "ID"\n",
714 k, newrow)) ;
715 if (--(Wp [newrow]) < 0)
716 {
717 /* :: pattern change (singleton row too lengthy) :: */
718 DEBUGm4 (("bad U singleton row (too lengthy)\n")) ;
719 return (FALSE) ; /* pattern changed */
720 }
721 *(Rpi [newrow]++) = k ;
722 *(Rpx [newrow]++) = x ;
723 }
724 else
725 {
726 /* this entry goes in a column of L */
727 DEBUG1 (("Singleton col of L: k "ID" newrow "ID"\n",
728 k, newrow)) ;
729 if (llen >= lilen)
730 {
731 DEBUGm4 (("bad L singleton col (too lengthy)\n")) ;
732 return (FALSE) ; /* pattern changed */
733 }
734 Li [llen] = newrow ;
735 Lval [llen] = x ;
736 llen++ ;
737 }
738 }
739
740 if (llen != lilen)
741 {
742 /* :: pattern change (singleton column too lengthy) :: */
743 DEBUGm4 (("bad L singleton col (too short)\n")) ;
744 return (FALSE) ; /* pattern changed */
745 }
746
747 /* scale the column of L */
748 if (llen > 0)
749 {
750 pivot_value = D [k] ;
751 UMF_scale (llen, pivot_value, Lval) ;
752 }
753
754 }
755
756 /* ---------------------------------------------------------------------- */
757 /* allocate the elements and copy the columns of A */
758 /* ---------------------------------------------------------------------- */
759
760 /* also apply the row and column pre-ordering. */
761 for (k = n1 ; k < n_col ; k++)
762 {
763 /* The newcol is k, which is what the name of the column is in the
764 * UMFPACK kernel. The user's name for the column is oldcol. */
765 oldcol = Cperm_init [k] ;
766
767 ASSERT (oldcol >= 0 && oldcol < n_col) ;
768
769 p2 = Ap [oldcol+1] ;
770
771 cdeg = Cdeg [k] ;
772 ASSERT (cdeg >= 0) ;
773
774 /* if fixQ: set Col_degree to 0 for the NON_PIVOTAL_COL macro */
775 Col_degree [k] = fixQ ? 0 : cdeg ;
776
777 /* get the element for this column (if any) */
778 e = k - n1 + 1 ;
779 if (k < n_col - nempty_col)
780 {
781 esize = Esize ? Esize [k-n1] : cdeg ;
782 if (E [e])
783 {
784 Int ncols, nrows ;
785 Unit *pp ;
786 pp = Memory + E [e] ;
787 GET_ELEMENT (ep, pp, Cols, Rows, ncols, nrows, C) ;
788 ASSERT (ncols == 1) ;
789 ASSERT (nrows == esize) ;
790 ASSERT (Cols [0] == k) ;
791 }
792 }
793 else
794 {
795 ASSERT (cdeg == 0) ;
796 esize = 0 ;
797 }
798
799 clen = 0 ;
800
801 for (pa = Ap [oldcol] ; pa < p2 ; pa++)
802 {
803 oldrow = Ai [pa] ;
804 newrow = Frpos [oldrow] ;
805 ASSIGN (x, Ax, Az, pa, split) ;
806
807 /* scale the value using the scale factors, Rs */
808 if (do_scale)
809 {
810 #ifndef NRECIPROCAL
811 if (do_recip)
812 {
813 /* multiply by the reciprocal */
814 SCALE (x, Rs [oldrow]) ;
815 }
816 else
817 #endif
818 {
819 /* divide instead */
820 SCALE_DIV (x, Rs [oldrow]) ;
821 }
822 }
823
824 rdeg = Rdeg [newrow] ;
825 if (newrow < n1 || rdeg > dense_row_threshold)
826 {
827 /* this entry goes in a row of U or into a dense row */
828 DEBUG1 (("Singleton/dense row of U: k "ID" newrow "ID"\n",
829 k, newrow)) ;
830 if (--(Wp [newrow]) < 0)
831 {
832 DEBUGm4 (("bad row of U or A (too lengthy)\n")) ;
833 return (FALSE) ; /* pattern changed */
834 }
835 *(Rpi [newrow]++) = k ;
836 *(Rpx [newrow]++) = x ;
837 }
838 else
839 {
840 /* this entry goes in an initial element */
841 DEBUG1 (("In element k "ID" e "ID" newrow "ID"\n",
842 k, e, newrow)) ;
843 if (clen >= esize)
844 {
845 DEBUGm4 (("bad A column (too lengthy)\n")) ;
846 return (FALSE) ; /* pattern changed */
847 }
848 ASSERT (E [e]) ;
849 ASSERT (k < n_col - nempty_col) ;
850 Rows [clen] = newrow ;
851 C [clen] = x ;
852 clen++ ;
853 #ifndef NDEBUG
854 if (Diagonal_map && (newrow == Diagonal_map [k]))
855 {
856 DEBUG0 (("Diagonal: old: row "ID" col "ID" : "
857 "new: row "ID" col "ID" : ",
858 oldrow, oldcol, newrow, k)) ;
859 EDEBUGk (0, x) ;
860 }
861 #endif
862 }
863 }
864
865 if (clen != esize)
866 {
867 /* :: pattern change (singleton column too short) :: */
868 DEBUGm4 (("bad A column (too short)\n")) ;
869 return (FALSE) ; /* pattern changed */
870 }
871 }
872
873 /* ---------------------------------------------------------------------- */
874 /* free the Rpi and Rpx workspace at the tail end of memory */
875 /* ---------------------------------------------------------------------- */
876
877 UMF_mem_free_tail_block (Numeric, rpi) ;
878 UMF_mem_free_tail_block (Numeric, rpx) ;
879
880 /* ---------------------------------------------------------------------- */
881 /* prune zeros and small entries from the singleton rows and columns */
882 /* ---------------------------------------------------------------------- */
883
884 if (n1 > 0)
885 {
886 pnew = Lip [0] ;
887 ASSERT (pnew == 1) ;
888 for (k = 0 ; k < n1 ; k++)
889 {
890 DEBUGm4 (("\nPrune singleton L col "ID"\n", k)) ;
891 pnew = packsp (pnew, &Lip [k], &Lilen [k], drop, droptol, Memory) ;
892 Numeric->lnz += Lilen [k] ;
893 DEBUGm4 (("\nPrune singleton U row "ID"\n", k)) ;
894 pnew = packsp (pnew, &Uip [k], &Uilen [k], drop, droptol, Memory) ;
895 Numeric->unz += Uilen [k] ;
896 }
897 /* free the unused space at the head of memory */
898 Numeric->ihead = pnew ;
899 }
900
901 /* ---------------------------------------------------------------------- */
902 /* initialize row degrees */
903 /* ---------------------------------------------------------------------- */
904
905 for (k = 0 ; k < n1 ; k++)
906 {
907 if (Wp [k] != 0)
908 {
909 /* :: pattern change (singleton row too short) :: */
910 DEBUGm4 (("bad U singleton row (too short)\n")) ;
911 return (FALSE) ; /* pattern changed */
912 }
913 }
914
915 for (k = n1 ; k < n_row ; k++)
916 {
917 DEBUG1 (("Initial row degree k "ID" oldrow "ID" Rdeg "ID"\n",
918 k, Rperm_init [k], Rdeg [k])) ;
919 rdeg = Rdeg [k] ;
920 Row_degree [k] = rdeg ;
921 if (rdeg > dense_row_threshold && Wp [k] != 0)
922 {
923 /* :: pattern change (dense row too short) :: */
924 DEBUGm4 (("bad dense row (too short)\n")) ;
925 return (FALSE) ; /* pattern changed */
926 }
927 }
928
929 #ifndef NDEBUG
930 if (prefer_diagonal)
931 {
932 Entry aij ;
933 Int *InvCperm, newcol ;
934 UMF_dump_diagonal_map (Diagonal_map, Diagonal_imap, nn) ;
935 InvCperm = (Int *) malloc (n_col * sizeof (Int)) ;
936 ASSERT (InvCperm != (Int *) NULL) ;
937 for (newcol = 0 ; newcol < n_col ; newcol++)
938 {
939 oldcol = Cperm_init [newcol] ;
940 InvCperm [oldcol] = newcol ;
941 }
942 DEBUGm3 (("Diagonal of P2*A:\n")) ;
943 for (oldcol = 0 ; oldcol < n_col ; oldcol++)
944 {
945 newcol = InvCperm [oldcol] ;
946 for (p = Ap [oldcol] ; p < Ap [oldcol+1] ; p++)
947 {
948 oldrow = Ai [p] ;
949 newrow = Frpos [oldrow] ;
950 ASSIGN (aij, Ax, Az, p, split) ;
951 if (newrow == Diagonal_map [newcol])
952 {
953 DEBUG0 (("old row "ID" col "ID" new row "ID" col "ID,
954 oldrow, oldcol, newrow, newcol)) ;
955 EDEBUGk (0, aij) ;
956 DEBUG0 ((" scaled ")) ;
957 if (do_scale)
958 {
959 #ifndef NRECIPROCAL
960 if (do_recip)
961 {
962 SCALE (aij, Rs [oldrow]) ;
963 }
964 else
965 #endif
966 {
967 SCALE_DIV (aij, Rs [oldrow]) ;
968 }
969 }
970 EDEBUGk (0, aij) ;
971 DEBUG0 (("\n")) ;
972 }
973 }
974 }
975 free (InvCperm) ;
976 }
977 #endif
978
979 Col_degree [n_col] = 0 ;
980
981 /* ---------------------------------------------------------------------- */
982 /* pack the element name space */
983 /* ---------------------------------------------------------------------- */
984
985 if (empty_elements)
986 {
987 Int e2 = 0 ;
988 DEBUG0 (("\n\n============= Packing element space\n")) ;
989 for (e = 1 ; e <= Work->nel ; e++)
990 {
991 if (E [e])
992 {
993 e2++ ;
994 E [e2] = E [e] ;
995 }
996 }
997 Work->nel = e2 ;
998 }
999
1000 #ifndef NDEBUG
1001 DEBUG0 (("Number of initial elements: "ID"\n", Work->nel)) ;
1002 for (e = 0 ; e <= Work->nel ; e++) UMF_dump_element (Numeric, Work,e,TRUE) ;
1003 #endif
1004
1005 for (e = Work->nel + 1 ; e < Work->elen ; e++)
1006 {
1007 E [e] = 0 ;
1008 }
1009
1010 /* Frpos no longer needed */
1011 for (row = 0 ; row <= n_row ; row++)
1012 {
1013 Frpos [row] = EMPTY ;
1014 }
1015
1016 /* clear Wp */
1017 for (i = 0 ; i <= nn ; i++)
1018 {
1019 Wp [i] = EMPTY ;
1020 }
1021
1022 DEBUG1 (("Kernel init head usage: "ID"\n", Numeric->ihead)) ;
1023
1024 /* ---------------------------------------------------------------------- */
1025 /* build the tuple lists */
1026 /* ---------------------------------------------------------------------- */
1027
1028 /* if the memory usage changes, then the pattern has changed */
1029
1030 (void) UMF_tuple_lengths (Numeric, Work, &unused) ;
1031 if (!UMF_build_tuples (Numeric, Work))
1032 {
1033 /* :: pattern change (out of memory in umf_build_tuples) :: */
1034 /* We ran out of memory, which can only mean that */
1035 /* the pattern (Ap and or Ai) has changed (gotten larger). */
1036 DEBUG0 (("Pattern has gotten larger - build tuples failed\n")) ;
1037 return (FALSE) ; /* pattern changed */
1038 }
1039
1040 Numeric->init_usage = Numeric->max_usage ;
1041
1042 /* ---------------------------------------------------------------------- */
1043 /* construct the row merge sets */
1044 /* ---------------------------------------------------------------------- */
1045
1046 for (i = 0 ; i <= Symbolic->nfr ; i++)
1047 {
1048 Work->Front_new1strow [i] = Symbolic->Front_1strow [i] ;
1049 }
1050
1051 #ifndef NDEBUG
1052 UMF_dump_rowmerge (Numeric, Symbolic, Work) ;
1053 DEBUG6 (("Column form of original matrix:\n")) ;
1054 UMF_dump_col_matrix (Ax,
1055 #ifdef COMPLEX
1056 Az,
1057 #endif
1058 Ai, Ap, n_row, n_col, nz) ;
1059 UMF_dump_memory (Numeric) ;
1060 UMF_dump_matrix (Numeric, Work, FALSE) ;
1061 DEBUG0 (("kernel init done...\n")) ;
1062 #endif
1063
1064 return (TRUE) ;
1065 }
1066