1 /* ========================================================================== */
2 /* === UMF_store_lu ========================================================= */
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 Store the LU factors. Called by the kernel.
12 Returns TRUE if successful, FALSE if out of memory.
13 */
14
15 #include "umf_internal.h"
16 #include "umf_store_lu.h"
17 #include "umf_mem_alloc_head_block.h"
18 #include "umf_get_memory.h"
19
20 /* ========================================================================== */
21
22 #ifdef DROP
UMF_store_lu_drop(NumericType * Numeric,WorkType * Work)23 GLOBAL Int UMF_store_lu_drop
24 #else
25 GLOBAL Int UMF_store_lu
26 #endif
27 (
28 NumericType *Numeric,
29 WorkType *Work
30 )
31 {
32 /* ---------------------------------------------------------------------- */
33 /* local variables */
34 /* ---------------------------------------------------------------------- */
35
36 Entry pivot_value ;
37 #ifdef DROP
38 double droptol ;
39 #endif
40 Entry *D, *Lval, *Uval, *Fl1, *Fl2, *Fu1, *Fu2,
41 *Flublock, *Flblock, *Fublock ;
42 Int i, k, fnr_curr, fnrows, fncols, row, col, pivrow, pivcol, *Frows,
43 *Fcols, *Lpattern, *Upattern, *Lpos, *Upos, llen, ulen, fnc_curr, fnpiv,
44 uilen, lnz, unz, nb, *Lilen,
45 *Uilen, *Lip, *Uip, *Li, *Ui, pivcol_position, newLchain, newUchain,
46 pivrow_position, p, size, lip, uip, lnzi, lnzx, unzx, lnz2i, lnz2x,
47 unz2i, unz2x, zero_pivot, *Pivrow, *Pivcol, kk,
48 Lnz [MAXNB] ;
49
50 #ifndef NDEBUG
51 Int *Col_degree, *Row_degree ;
52 #endif
53
54 #ifdef DROP
55 Int all_lnz, all_unz ;
56 droptol = Numeric->droptol ;
57 #endif
58
59 /* ---------------------------------------------------------------------- */
60 /* get parameters */
61 /* ---------------------------------------------------------------------- */
62
63 fnrows = Work->fnrows ;
64 fncols = Work->fncols ;
65 fnpiv = Work->fnpiv ;
66
67 Lpos = Numeric->Lpos ;
68 Upos = Numeric->Upos ;
69 Lilen = Numeric->Lilen ;
70 Uilen = Numeric->Uilen ;
71
72 Lip = Numeric->Lip ;
73 Uip = Numeric->Uip ;
74 D = Numeric->D ;
75
76 Flublock = Work->Flublock ;
77 Flblock = Work->Flblock ;
78 Fublock = Work->Fublock ;
79
80 fnr_curr = Work->fnr_curr ;
81 fnc_curr = Work->fnc_curr ;
82 Frows = Work->Frows ;
83 Fcols = Work->Fcols ;
84
85 #ifndef NDEBUG
86 Col_degree = Numeric->Cperm ; /* for NON_PIVOTAL_COL macro */
87 Row_degree = Numeric->Rperm ; /* for NON_PIVOTAL_ROW macro */
88 #endif
89
90 Lpattern = Work->Lpattern ;
91 llen = Work->llen ;
92 Upattern = Work->Upattern ;
93 ulen = Work->ulen ;
94
95 nb = Work->nb ;
96
97 #ifndef NDEBUG
98 DEBUG1 (("\nSTORE LU: fnrows "ID
99 " fncols "ID"\n", fnrows, fncols)) ;
100
101 DEBUG2 (("\nFrontal matrix, including all space:\n"
102 "fnr_curr "ID" fnc_curr "ID" nb "ID"\n"
103 "fnrows "ID" fncols "ID" fnpiv "ID"\n",
104 fnr_curr, fnc_curr, nb, fnrows, fncols, fnpiv)) ;
105
106 DEBUG2 (("\nJust the active part:\n")) ;
107 DEBUG7 (("C block: ")) ;
108 UMF_dump_dense (Work->Fcblock, fnr_curr, fnrows, fncols) ;
109 DEBUG7 (("L block: ")) ;
110 UMF_dump_dense (Work->Flblock, fnr_curr, fnrows, fnpiv);
111 DEBUG7 (("U' block: ")) ;
112 UMF_dump_dense (Work->Fublock, fnc_curr, fncols, fnpiv) ;
113 DEBUG7 (("LU block: ")) ;
114 UMF_dump_dense (Work->Flublock, nb, fnpiv, fnpiv) ;
115 DEBUG7 (("Current frontal matrix: (prior to store LU)\n")) ;
116 UMF_dump_current_front (Numeric, Work, TRUE) ;
117 #endif
118
119 Pivrow = Work->Pivrow ;
120 Pivcol = Work->Pivcol ;
121
122 /* ---------------------------------------------------------------------- */
123 /* store the columns of L */
124 /* ---------------------------------------------------------------------- */
125
126 for (kk = 0 ; kk < fnpiv ; kk++)
127 {
128
129 /* ------------------------------------------------------------------ */
130 /* one more pivot row and column is being stored into L and U */
131 /* ------------------------------------------------------------------ */
132
133 k = Work->npiv + kk ;
134
135 /* ------------------------------------------------------------------ */
136 /* find the kth pivot row and pivot column */
137 /* ------------------------------------------------------------------ */
138
139 pivrow = Pivrow [kk] ;
140 pivcol = Pivcol [kk] ;
141
142 #ifndef NDEBUG
143 ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
144 ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
145
146 DEBUGm4 ((
147 "\n -------------------------------------------------------------"
148 "Store LU: step " ID"\n", k)) ;
149 ASSERT (k < MIN (Work->n_row, Work->n_col)) ;
150 DEBUG2 (("Store column of L, k = "ID", llen "ID"\n", k, llen)) ;
151 for (i = 0 ; i < llen ; i++)
152 {
153 row = Lpattern [i] ;
154 ASSERT (row >= 0 && row < Work->n_row) ;
155 DEBUG2 ((" Lpattern["ID"] "ID" Lpos "ID, i, row, Lpos [row])) ;
156 if (row == pivrow) DEBUG2 ((" <- pivot row")) ;
157 DEBUG2 (("\n")) ;
158 ASSERT (i == Lpos [row]) ;
159 }
160 #endif
161
162 /* ------------------------------------------------------------------ */
163 /* remove pivot row from L */
164 /* ------------------------------------------------------------------ */
165
166 /* remove pivot row index from current column of L */
167 /* if a new Lchain starts, then all entries are removed later */
168 DEBUG2 (("Removing pivrow from Lpattern, k = "ID"\n", k)) ;
169 ASSERT (!NON_PIVOTAL_ROW (pivrow)) ;
170 pivrow_position = Lpos [pivrow] ;
171 if (pivrow_position != EMPTY)
172 {
173 /* place the last entry in the column in the */
174 /* position of the pivot row index */
175 ASSERT (pivrow == Lpattern [pivrow_position]) ;
176 row = Lpattern [--llen] ;
177 /* ASSERT (NON_PIVOTAL_ROW (row)) ; */
178 Lpattern [pivrow_position] = row ;
179 Lpos [row] = pivrow_position ;
180 Lpos [pivrow] = EMPTY ;
181 }
182
183 /* ------------------------------------------------------------------ */
184 /* store the pivot value, for the diagonal matrix D */
185 /* ------------------------------------------------------------------ */
186
187 /* kk-th column of LU block */
188 Fl1 = Flublock + kk * nb ;
189
190 /* kk-th column of L in the L block */
191 Fl2 = Flblock + kk * fnr_curr ;
192
193 /* kk-th pivot in frontal matrix located in Flublock [kk, kk] */
194 pivot_value = Fl1 [kk] ;
195
196 D [k] = pivot_value ;
197 zero_pivot = IS_ZERO (pivot_value) ;
198
199 DEBUG4 (("Pivot D["ID"]=", k)) ;
200 EDEBUG4 (pivot_value) ;
201 DEBUG4 (("\n")) ;
202
203 /* ------------------------------------------------------------------ */
204 /* count nonzeros in kth column of L */
205 /* ------------------------------------------------------------------ */
206
207 lnz = 0 ;
208 lnz2i = 0 ;
209 lnz2x = llen ;
210
211 #ifdef DROP
212 all_lnz = 0 ;
213
214 for (i = kk + 1 ; i < fnpiv ; i++)
215 {
216 Entry x ;
217 double s ;
218 x = Fl1 [i] ;
219 if (IS_ZERO (x)) continue ;
220 all_lnz++ ;
221 APPROX_ABS (s, x) ;
222 if (s <= droptol) continue ;
223 lnz++ ;
224 if (Lpos [Pivrow [i]] == EMPTY) lnz2i++ ;
225 }
226
227 for (i = 0 ; i < fnrows ; i++)
228 {
229 Entry x ;
230 double s ;
231 x = Fl2 [i] ;
232 if (IS_ZERO (x)) continue ;
233 all_lnz++ ;
234 APPROX_ABS (s, x) ;
235 if (s <= droptol) continue ;
236 lnz++ ;
237 if (Lpos [Frows [i]] == EMPTY) lnz2i++ ;
238 }
239
240 #else
241
242 for (i = kk + 1 ; i < fnpiv ; i++)
243 {
244 if (IS_ZERO (Fl1 [i])) continue ;
245 lnz++ ;
246 if (Lpos [Pivrow [i]] == EMPTY) lnz2i++ ;
247 }
248
249 for (i = 0 ; i < fnrows ; i++)
250 {
251 if (IS_ZERO (Fl2 [i])) continue ;
252 lnz++ ;
253 if (Lpos [Frows [i]] == EMPTY) lnz2i++ ;
254 }
255
256 #endif
257
258 lnz2x += lnz2i ;
259
260 /* determine if we start a new Lchain or continue the old one */
261 if (llen == 0 || zero_pivot)
262 {
263 /* llen == 0 means there is no prior Lchain */
264 /* D [k] == 0 means the pivot column is empty */
265 newLchain = TRUE ;
266 }
267 else
268 {
269 newLchain =
270 /* storage for starting a new Lchain */
271 UNITS (Entry, lnz) + UNITS (Int, lnz)
272 <=
273 /* storage for continuing a prior Lchain */
274 UNITS (Entry, lnz2x) + UNITS (Int, lnz2i) ;
275 }
276
277 if (newLchain)
278 {
279 /* start a new chain for column k of L */
280 DEBUG2 (("Start new Lchain, k = "ID"\n", k)) ;
281
282 pivrow_position = EMPTY ;
283
284 /* clear the prior Lpattern */
285 for (i = 0 ; i < llen ; i++)
286 {
287 row = Lpattern [i] ;
288 Lpos [row] = EMPTY ;
289 }
290 llen = 0 ;
291
292 lnzi = lnz ;
293 lnzx = lnz ;
294 }
295 else
296 {
297 /* continue the prior Lchain */
298 DEBUG2 (("Continue Lchain, k = "ID"\n", k)) ;
299 lnzi = lnz2i ;
300 lnzx = lnz2x ;
301 }
302
303 /* ------------------------------------------------------------------ */
304 /* allocate space for the column of L */
305 /* ------------------------------------------------------------------ */
306
307 size = UNITS (Int, lnzi) + UNITS (Entry, lnzx) ;
308
309 #ifndef NDEBUG
310 UMF_allocfail = FALSE ;
311 if (UMF_gprob > 0)
312 {
313 double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
314 DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
315 UMF_allocfail = rrr < UMF_gprob ;
316 if (UMF_allocfail) DEBUGm2 (("Random garbage coll. (store LU)\n"));
317 }
318 #endif
319
320 p = UMF_mem_alloc_head_block (Numeric, size) ;
321 if (!p)
322 {
323 Int r2, c2 ;
324 /* Do garbage collection, realloc, and try again. */
325 /* Note that there are pivot rows/columns in current front. */
326 if (Work->do_grow)
327 {
328 /* full compaction of current frontal matrix, since
329 * UMF_grow_front will be called next anyway. */
330 r2 = fnrows ;
331 c2 = fncols ;
332 }
333 else
334 {
335 /* partial compaction. */
336 r2 = MAX (fnrows, Work->fnrows_new + 1) ;
337 c2 = MAX (fncols, Work->fncols_new + 1) ;
338 }
339 DEBUGm3 (("get_memory from umf_store_lu:\n")) ;
340 if (!UMF_get_memory (Numeric, Work, size, r2, c2, TRUE))
341 {
342 DEBUGm4 (("out of memory: store LU (1)\n")) ;
343 return (FALSE) ; /* out of memory */
344 }
345 p = UMF_mem_alloc_head_block (Numeric, size) ;
346 if (!p)
347 {
348 DEBUGm4 (("out of memory: store LU (2)\n")) ;
349 return (FALSE) ; /* out of memory */
350 }
351 /* garbage collection may have moved the current front */
352 fnc_curr = Work->fnc_curr ;
353 fnr_curr = Work->fnr_curr ;
354 Flublock = Work->Flublock ;
355 Flblock = Work->Flblock ;
356 Fublock = Work->Fublock ;
357 Fl1 = Flublock + kk * nb ;
358 Fl2 = Flblock + kk * fnr_curr ;
359 }
360
361 /* ------------------------------------------------------------------ */
362 /* store the column of L */
363 /* ------------------------------------------------------------------ */
364
365 lip = p ;
366
367 Li = (Int *) (Numeric->Memory + p) ;
368 p += UNITS (Int, lnzi) ;
369 Lval = (Entry *) (Numeric->Memory + p) ;
370 p += UNITS (Entry, lnzx) ;
371
372 for (i = 0 ; i < lnzx ; i++)
373 {
374 CLEAR (Lval [i]) ;
375 }
376
377 /* store the numerical entries */
378
379 if (newLchain)
380 {
381 /* flag the first column in the Lchain by negating Lip [k] */
382 lip = -lip ;
383
384 ASSERT (llen == 0) ;
385
386 #ifdef DROP
387
388 for (i = kk + 1 ; i < fnpiv ; i++)
389 {
390 Entry x ;
391 double s ;
392 Int row2, pos ;
393 x = Fl1 [i] ;
394 APPROX_ABS (s, x) ;
395 if (s <= droptol) continue ;
396 row2 = Pivrow [i] ;
397 pos = llen++ ;
398 Lpattern [pos] = row2 ;
399 Lpos [row2] = pos ;
400 Li [pos] = row2 ;
401 Lval [pos] = x ;
402 }
403
404 for (i = 0 ; i < fnrows ; i++)
405 {
406 Entry x ;
407 double s ;
408 Int row2, pos ;
409 x = Fl2 [i] ;
410 APPROX_ABS (s, x) ;
411 if (s <= droptol) continue ;
412 row2 = Frows [i] ;
413 pos = llen++ ;
414 Lpattern [pos] = row2 ;
415 Lpos [row2] = pos ;
416 Li [pos] = row2 ;
417 Lval [pos] = x ;
418 }
419
420 #else
421
422 for (i = kk + 1 ; i < fnpiv ; i++)
423 {
424 Entry x ;
425 Int row2, pos ;
426 x = Fl1 [i] ;
427 if (IS_ZERO (x)) continue ;
428 row2 = Pivrow [i] ;
429 pos = llen++ ;
430 Lpattern [pos] = row2 ;
431 Lpos [row2] = pos ;
432 Li [pos] = row2 ;
433 Lval [pos] = x ;
434 }
435
436 for (i = 0 ; i < fnrows ; i++)
437 {
438 Entry x ;
439 Int row2, pos ;
440 x = Fl2 [i] ;
441 if (IS_ZERO (x)) continue ;
442 row2 = Frows [i] ;
443 pos = llen++ ;
444 Lpattern [pos] = row2 ;
445 Lpos [row2] = pos ;
446 Li [pos] = row2 ;
447 Lval [pos] = x ;
448 }
449
450 #endif
451
452 }
453 else
454 {
455 ASSERT (llen > 0) ;
456
457 #ifdef DROP
458
459 for (i = kk + 1 ; i < fnpiv ; i++)
460 {
461 Entry x ;
462 double s ;
463 Int row2, pos ;
464 x = Fl1 [i] ;
465 APPROX_ABS (s, x) ;
466 if (s <= droptol) continue ;
467 row2 = Pivrow [i] ;
468 pos = Lpos [row2] ;
469 if (pos == EMPTY)
470 {
471 pos = llen++ ;
472 Lpattern [pos] = row2 ;
473 Lpos [row2] = pos ;
474 *Li++ = row2 ;
475 }
476 Lval [pos] = x ;
477 }
478
479 for (i = 0 ; i < fnrows ; i++)
480 {
481 Entry x ;
482 double s ;
483 Int row2, pos ;
484 x = Fl2 [i] ;
485 APPROX_ABS (s, x) ;
486 if (s <= droptol) continue ;
487 row2 = Frows [i] ;
488 pos = Lpos [row2] ;
489 if (pos == EMPTY)
490 {
491 pos = llen++ ;
492 Lpattern [pos] = row2 ;
493 Lpos [row2] = pos ;
494 *Li++ = row2 ;
495 }
496 Lval [pos] = x ;
497 }
498
499 #else
500
501 for (i = kk + 1 ; i < fnpiv ; i++)
502 {
503 Entry x ;
504 Int row2, pos ;
505 x = Fl1 [i] ;
506 if (IS_ZERO (x)) continue ;
507 row2 = Pivrow [i] ;
508 pos = Lpos [row2] ;
509 if (pos == EMPTY)
510 {
511 pos = llen++ ;
512 Lpattern [pos] = row2 ;
513 Lpos [row2] = pos ;
514 *Li++ = row2 ;
515 }
516 Lval [pos] = x ;
517 }
518
519 for (i = 0 ; i < fnrows ; i++)
520 {
521 Entry x ;
522 Int row2, pos ;
523 x = Fl2 [i] ;
524 if (IS_ZERO (x)) continue ;
525 row2 = Frows [i] ;
526 pos = Lpos [row2] ;
527 if (pos == EMPTY)
528 {
529 pos = llen++ ;
530 Lpattern [pos] = row2 ;
531 Lpos [row2] = pos ;
532 *Li++ = row2 ;
533 }
534 Lval [pos] = x ;
535 }
536
537 #endif
538
539 }
540 DEBUG4 (("llen "ID" lnzx "ID"\n", llen, lnzx)) ;
541 ASSERT (llen == lnzx) ;
542 ASSERT (lnz <= llen) ;
543 DEBUG4 (("lnz "ID" \n", lnz)) ;
544
545 #ifdef DROP
546
547 DEBUG4 (("all_lnz "ID" \n", all_lnz)) ;
548 ASSERT (lnz <= all_lnz) ;
549 Numeric->lnz += lnz ;
550 Numeric->all_lnz += all_lnz ;
551 Lnz [kk] = all_lnz ;
552
553 #else
554
555 Numeric->lnz += lnz ;
556 Numeric->all_lnz += lnz ;
557 Lnz [kk] = lnz ;
558 #endif
559
560 Numeric->nLentries += lnzx ;
561 Work->llen = llen ;
562 Numeric->isize += lnzi ;
563
564 /* ------------------------------------------------------------------ */
565 /* the pivot column is fully assembled and scaled, and is now the */
566 /* k-th column of L */
567 /* ------------------------------------------------------------------ */
568
569 Lpos [pivrow] = pivrow_position ; /* not aliased */
570 Lip [pivcol] = lip ; /* aliased with Col_tuples */
571 Lilen [pivcol] = lnzi ; /* aliased with Col_tlen */
572
573 }
574
575 /* ---------------------------------------------------------------------- */
576 /* store the rows of U */
577 /* ---------------------------------------------------------------------- */
578
579 for (kk = 0 ; kk < fnpiv ; kk++)
580 {
581
582 /* ------------------------------------------------------------------ */
583 /* one more pivot row and column is being stored into L and U */
584 /* ------------------------------------------------------------------ */
585
586 k = Work->npiv + kk ;
587
588 /* ------------------------------------------------------------------ */
589 /* find the kth pivot row and pivot column */
590 /* ------------------------------------------------------------------ */
591
592 pivrow = Pivrow [kk] ;
593 pivcol = Pivcol [kk] ;
594
595 #ifndef NDEBUG
596 ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
597 ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
598
599 DEBUG2 (("Store row of U, k = "ID", ulen "ID"\n", k, ulen)) ;
600 for (i = 0 ; i < ulen ; i++)
601 {
602 col = Upattern [i] ;
603 DEBUG2 ((" Upattern["ID"] "ID, i, col)) ;
604 if (col == pivcol) DEBUG2 ((" <- pivot col")) ;
605 DEBUG2 (("\n")) ;
606 ASSERT (col >= 0 && col < Work->n_col) ;
607 ASSERT (i == Upos [col]) ;
608 }
609 #endif
610
611 /* ------------------------------------------------------------------ */
612 /* get the pivot value, for the diagonal matrix D */
613 /* ------------------------------------------------------------------ */
614
615 zero_pivot = IS_ZERO (D [k]) ;
616
617 /* ------------------------------------------------------------------ */
618 /* count the nonzeros in the row of U */
619 /* ------------------------------------------------------------------ */
620
621 /* kk-th row of U in the LU block */
622 Fu1 = Flublock + kk ;
623
624 /* kk-th row of U in the U block */
625 Fu2 = Fublock + kk * fnc_curr ;
626
627 unz = 0 ;
628 unz2i = 0 ;
629 unz2x = ulen ;
630 DEBUG2 (("unz2x is "ID", lnzx "ID"\n", unz2x, lnzx)) ;
631
632 /* if row k does not end a Uchain, pivcol not included in ulen */
633 ASSERT (!NON_PIVOTAL_COL (pivcol)) ;
634 pivcol_position = Upos [pivcol] ;
635 if (pivcol_position != EMPTY)
636 {
637 unz2x-- ;
638 DEBUG2 (("(exclude pivcol) unz2x is now "ID"\n", unz2x)) ;
639 }
640
641 ASSERT (unz2x >= 0) ;
642
643 #ifdef DROP
644 all_unz = 0 ;
645
646 for (i = kk + 1 ; i < fnpiv ; i++)
647 {
648 Entry x ;
649 double s ;
650 x = Fu1 [i*nb] ;
651 if (IS_ZERO (x)) continue ;
652 all_unz++ ;
653 APPROX_ABS (s, x) ;
654 if (s <= droptol) continue ;
655 unz++ ;
656 if (Upos [Pivcol [i]] == EMPTY) unz2i++ ;
657 }
658
659 for (i = 0 ; i < fncols ; i++)
660 {
661 Entry x ;
662 double s ;
663 x = Fu2 [i] ;
664 if (IS_ZERO (x)) continue ;
665 all_unz++ ;
666 APPROX_ABS (s, x) ;
667 if (s <= droptol) continue ;
668 unz++ ;
669 if (Upos [Fcols [i]] == EMPTY) unz2i++ ;
670 }
671
672 #else
673
674 for (i = kk + 1 ; i < fnpiv ; i++)
675 {
676 if (IS_ZERO (Fu1 [i*nb])) continue ;
677 unz++ ;
678 if (Upos [Pivcol [i]] == EMPTY) unz2i++ ;
679 }
680
681 for (i = 0 ; i < fncols ; i++)
682 {
683 if (IS_ZERO (Fu2 [i])) continue ;
684 unz++ ;
685 if (Upos [Fcols [i]] == EMPTY) unz2i++ ;
686 }
687
688 #endif
689
690 unz2x += unz2i ;
691
692 ASSERT (IMPLIES (k == 0, ulen == 0)) ;
693
694 /* determine if we start a new Uchain or continue the old one */
695 if (ulen == 0 || zero_pivot)
696 {
697 /* ulen == 0 means there is no prior Uchain */
698 /* D [k] == 0 means the matrix is singular (pivot row might */
699 /* not be empty, however, but start a new Uchain to prune zero */
700 /* entries for the deg > 0 test in UMF_u*solve) */
701 newUchain = TRUE ;
702 }
703 else
704 {
705 newUchain =
706 /* approximate storage for starting a new Uchain */
707 UNITS (Entry, unz) + UNITS (Int, unz)
708 <=
709 /* approximate storage for continuing a prior Uchain */
710 UNITS (Entry, unz2x) + UNITS (Int, unz2i) ;
711
712 /* this would be exact, except for the Int to Unit rounding, */
713 /* because the Upattern is stored only at the end of the Uchain */
714 }
715
716 /* ------------------------------------------------------------------ */
717 /* allocate space for the row of U */
718 /* ------------------------------------------------------------------ */
719
720 size = 0 ;
721 if (newUchain)
722 {
723 /* store the pattern of the last row in the prior Uchain */
724 size += UNITS (Int, ulen) ;
725 unzx = unz ;
726 }
727 else
728 {
729 unzx = unz2x ;
730 }
731 size += UNITS (Entry, unzx) ;
732
733 #ifndef NDEBUG
734 UMF_allocfail = FALSE ;
735 if (UMF_gprob > 0)
736 {
737 double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
738 DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
739 UMF_allocfail = rrr < UMF_gprob ;
740 if (UMF_allocfail) DEBUGm2 (("Random garbage coll. (store LU)\n"));
741 }
742 #endif
743
744 p = UMF_mem_alloc_head_block (Numeric, size) ;
745 if (!p)
746 {
747 Int r2, c2 ;
748 /* Do garbage collection, realloc, and try again. */
749 /* Note that there are pivot rows/columns in current front. */
750 if (Work->do_grow)
751 {
752 /* full compaction of current frontal matrix, since
753 * UMF_grow_front will be called next anyway. */
754 r2 = fnrows ;
755 c2 = fncols ;
756 }
757 else
758 {
759 /* partial compaction. */
760 r2 = MAX (fnrows, Work->fnrows_new + 1) ;
761 c2 = MAX (fncols, Work->fncols_new + 1) ;
762 }
763 DEBUGm3 (("get_memory from umf_store_lu:\n")) ;
764 if (!UMF_get_memory (Numeric, Work, size, r2, c2, TRUE))
765 {
766 /* :: get memory, column of L :: */
767 DEBUGm4 (("out of memory: store LU (1)\n")) ;
768 return (FALSE) ; /* out of memory */
769 }
770 p = UMF_mem_alloc_head_block (Numeric, size) ;
771 if (!p)
772 {
773 /* :: out of memory, column of U :: */
774 DEBUGm4 (("out of memory: store LU (2)\n")) ;
775 return (FALSE) ; /* out of memory */
776 }
777 /* garbage collection may have moved the current front */
778 fnc_curr = Work->fnc_curr ;
779 fnr_curr = Work->fnr_curr ;
780 Flublock = Work->Flublock ;
781 Flblock = Work->Flblock ;
782 Fublock = Work->Fublock ;
783 Fu1 = Flublock + kk ;
784 Fu2 = Fublock + kk * fnc_curr ;
785 }
786
787 /* ------------------------------------------------------------------ */
788 /* store the row of U */
789 /* ------------------------------------------------------------------ */
790
791 uip = p ;
792
793 if (newUchain)
794 {
795 /* starting a new Uchain - flag this by negating Uip [k] */
796 uip = -uip ;
797 DEBUG2 (("Start new Uchain, k = "ID"\n", k)) ;
798
799 pivcol_position = EMPTY ;
800
801 /* end the prior Uchain */
802 /* save the current Upattern, and then */
803 /* clear it and start a new Upattern */
804 DEBUG2 (("Ending prior chain, k-1 = "ID"\n", k-1)) ;
805 uilen = ulen ;
806 Ui = (Int *) (Numeric->Memory + p) ;
807 Numeric->isize += ulen ;
808 p += UNITS (Int, ulen) ;
809 for (i = 0 ; i < ulen ; i++)
810 {
811 col = Upattern [i] ;
812 ASSERT (col >= 0 && col < Work->n_col) ;
813 Upos [col] = EMPTY ;
814 Ui [i] = col ;
815 }
816
817 ulen = 0 ;
818
819 }
820 else
821 {
822 /* continue the prior Uchain */
823 DEBUG2 (("Continue Uchain, k = "ID"\n", k)) ;
824 ASSERT (k > 0) ;
825
826 /* remove pivot col index from current row of U */
827 /* if a new Uchain starts, then all entries are removed later */
828 DEBUG2 (("Removing pivcol from Upattern, k = "ID"\n", k)) ;
829
830 if (pivcol_position != EMPTY)
831 {
832 /* place the last entry in the row in the */
833 /* position of the pivot col index */
834 ASSERT (pivcol == Upattern [pivcol_position]) ;
835 col = Upattern [--ulen] ;
836 ASSERT (col >= 0 && col < Work->n_col) ;
837 Upattern [pivcol_position] = col ;
838 Upos [col] = pivcol_position ;
839 Upos [pivcol] = EMPTY ;
840 }
841
842 /* this row continues the Uchain. Keep track of how much */
843 /* to trim from the k-th length to get the length of the */
844 /* (k-1)st row of U */
845 uilen = unz2i ;
846
847 }
848
849 Uval = (Entry *) (Numeric->Memory + p) ;
850 /* p += UNITS (Entry, unzx), no need to increment p */
851
852 for (i = 0 ; i < unzx ; i++)
853 {
854 CLEAR (Uval [i]) ;
855 }
856
857 if (newUchain)
858 {
859 ASSERT (ulen == 0) ;
860
861 #ifdef DROP
862
863 for (i = kk + 1 ; i < fnpiv ; i++)
864 {
865 Entry x ;
866 double s ;
867 Int col2, pos ;
868 x = Fu1 [i*nb] ;
869 APPROX_ABS (s, x) ;
870 if (s <= droptol) continue ;
871 col2 = Pivcol [i] ;
872 pos = ulen++ ;
873 Upattern [pos] = col2 ;
874 Upos [col2] = pos ;
875 Uval [pos] = x ;
876 }
877
878 for (i = 0 ; i < fncols ; i++)
879 {
880 Entry x ;
881 double s ;
882 Int col2, pos ;
883 x = Fu2 [i] ;
884 APPROX_ABS (s, x) ;
885 if (s <= droptol) continue ;
886 col2 = Fcols [i] ;
887 pos = ulen++ ;
888 Upattern [pos] = col2 ;
889 Upos [col2] = pos ;
890 Uval [pos] = x ;
891 }
892
893 #else
894
895 for (i = kk + 1 ; i < fnpiv ; i++)
896 {
897 Entry x ;
898 Int col2, pos ;
899 x = Fu1 [i*nb] ;
900 if (IS_ZERO (x)) continue ;
901 col2 = Pivcol [i] ;
902 pos = ulen++ ;
903 Upattern [pos] = col2 ;
904 Upos [col2] = pos ;
905 Uval [pos] = x ;
906 }
907
908 for (i = 0 ; i < fncols ; i++)
909 {
910 Entry x ;
911 Int col2, pos ;
912 x = Fu2 [i] ;
913 if (IS_ZERO (x)) continue ;
914 col2 = Fcols [i] ;
915 pos = ulen++ ;
916 Upattern [pos] = col2 ;
917 Upos [col2] = pos ;
918 Uval [pos] = x ;
919 }
920
921 #endif
922
923 }
924 else
925 {
926
927 ASSERT (ulen > 0) ;
928
929 /* store the numerical entries and find new nonzeros */
930
931 #ifdef DROP
932
933 for (i = kk + 1 ; i < fnpiv ; i++)
934 {
935 Entry x ;
936 double s ;
937 Int col2, pos ;
938 x = Fu1 [i*nb] ;
939 APPROX_ABS (s, x) ;
940 if (s <= droptol) continue ;
941 col2 = Pivcol [i] ;
942 pos = Upos [col2] ;
943 if (pos == EMPTY)
944 {
945 pos = ulen++ ;
946 Upattern [pos] = col2 ;
947 Upos [col2] = pos ;
948 }
949 Uval [pos] = x ;
950 }
951
952 for (i = 0 ; i < fncols ; i++)
953 {
954 Entry x ;
955 double s ;
956 Int col2, pos ;
957 x = Fu2 [i] ;
958 APPROX_ABS (s, x) ;
959 if (s <= droptol) continue ;
960 col2 = Fcols [i] ;
961 pos = Upos [col2] ;
962 if (pos == EMPTY)
963 {
964 pos = ulen++ ;
965 Upattern [pos] = col2 ;
966 Upos [col2] = pos ;
967 }
968 Uval [pos] = x ;
969 }
970
971 #else
972
973 for (i = kk + 1 ; i < fnpiv ; i++)
974 {
975 Entry x ;
976 Int col2, pos ;
977 x = Fu1 [i*nb] ;
978 if (IS_ZERO (x)) continue ;
979 col2 = Pivcol [i] ;
980 pos = Upos [col2] ;
981 if (pos == EMPTY)
982 {
983 pos = ulen++ ;
984 Upattern [pos] = col2 ;
985 Upos [col2] = pos ;
986 }
987 Uval [pos] = x ;
988 }
989
990 for (i = 0 ; i < fncols ; i++)
991 {
992 Entry x ;
993 Int col2, pos ;
994 x = Fu2 [i] ;
995 if (IS_ZERO (x)) continue ;
996 col2 = Fcols [i] ;
997 pos = Upos [col2] ;
998 if (pos == EMPTY)
999 {
1000 pos = ulen++ ;
1001 Upattern [pos] = col2 ;
1002 Upos [col2] = pos ;
1003 }
1004 Uval [pos] = x ;
1005 }
1006
1007 #endif
1008
1009 }
1010
1011 ASSERT (ulen == unzx) ;
1012 ASSERT (unz <= ulen) ;
1013 DEBUG4 (("unz "ID" \n", unz)) ;
1014
1015 #ifdef DROP
1016
1017 DEBUG4 (("all_unz "ID" \n", all_unz)) ;
1018 ASSERT (unz <= all_unz) ;
1019 Numeric->unz += unz ;
1020 Numeric->all_unz += all_unz ;
1021 /* count the "true" flops, based on LU pattern only */
1022 Numeric->flops += DIV_FLOPS * Lnz [kk] /* scale pivot column */
1023 + MULTSUB_FLOPS * (Lnz [kk] * all_unz) ; /* outer product */
1024
1025 #else
1026
1027 Numeric->unz += unz ;
1028 Numeric->all_unz += unz ;
1029 /* count the "true" flops, based on LU pattern only */
1030 Numeric->flops += DIV_FLOPS * Lnz [kk] /* scale pivot column */
1031 + MULTSUB_FLOPS * (Lnz [kk] * unz) ; /* outer product */
1032 #endif
1033
1034 Numeric->nUentries += unzx ;
1035 Work->ulen = ulen ;
1036 DEBUG1 (("Work->ulen = "ID" at end of pivot step, k: "ID"\n", ulen, k));
1037
1038 /* ------------------------------------------------------------------ */
1039 /* the pivot row is the k-th row of U */
1040 /* ------------------------------------------------------------------ */
1041
1042 Upos [pivcol] = pivcol_position ; /* not aliased */
1043 Uip [pivrow] = uip ; /* aliased with Row_tuples */
1044 Uilen [pivrow] = uilen ; /* aliased with Row_tlen */
1045
1046 }
1047
1048 /* ---------------------------------------------------------------------- */
1049 /* no more pivots in frontal working array */
1050 /* ---------------------------------------------------------------------- */
1051
1052 Work->npiv += fnpiv ;
1053 Work->fnpiv = 0 ;
1054 Work->fnzeros = 0 ;
1055 return (TRUE) ;
1056 }
1057