1 /* ========================================================================== */
2 /* === UMF_assemble ========================================================= */
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 /* Degree update and numerical assembly. This is compiled twice (with and
12 * without FIXQ) for each real/complex int/UF_long version, for a total of 8
13 * versions.*/
14
15 #include "umf_internal.h"
16 #include "umf_mem_free_tail_block.h"
17
18 /* ========================================================================== */
19 /* === row_assemble ========================================================= */
20 /* ========================================================================== */
21
row_assemble(Int row,NumericType * Numeric,WorkType * Work)22 PRIVATE void row_assemble
23 (
24 Int row,
25 NumericType *Numeric,
26 WorkType *Work
27 )
28 {
29
30 Entry *S, *Fcblock, *Frow ;
31 Int tpi, e, *E, *Fcpos, *Frpos, *Row_degree, *Row_tuples, *Row_tlen, rdeg0,
32 f, nrows, ncols, *Rows, *Cols, col, ncolsleft, j ;
33 Tuple *tp, *tp1, *tp2, *tpend ;
34 Unit *Memory, *p ;
35 Element *ep ;
36
37 #ifndef FIXQ
38 Int *Col_degree ;
39 Col_degree = Numeric->Cperm ;
40 #endif
41
42 Row_tuples = Numeric->Uip ;
43 tpi = Row_tuples [row] ;
44 if (!tpi) return ;
45
46 Memory = Numeric->Memory ;
47 E = Work->E ;
48 Fcpos = Work->Fcpos ;
49 Frpos = Work->Frpos ;
50 Row_degree = Numeric->Rperm ;
51 Row_tlen = Numeric->Uilen ;
52 E = Work->E ;
53 Memory = Numeric->Memory ;
54 rdeg0 = Work->rdeg0 ;
55 Fcblock = Work->Fcblock ;
56
57 #ifndef NDEBUG
58 DEBUG6 (("SCAN2-row: " ID "\n", row)) ;
59 UMF_dump_rowcol (0, Numeric, Work, row, FALSE) ;
60 #endif
61
62 ASSERT (NON_PIVOTAL_ROW (row)) ;
63
64 tp = (Tuple *) (Memory + tpi) ;
65 tp1 = tp ;
66 tp2 = tp ;
67 tpend = tp + Row_tlen [row] ;
68 for ( ; tp < tpend ; tp++)
69 {
70 e = tp->e ;
71 ASSERT (e > 0 && e <= Work->nel) ;
72 if (!E [e]) continue ; /* element already deallocated */
73 f = tp->f ;
74 p = Memory + E [e] ;
75 ep = (Element *) p ;
76 p += UNITS (Element, 1) ;
77 Cols = (Int *) p ;
78 Rows = Cols + ep->ncols ;
79 if (Rows [f] == EMPTY) continue ; /* row already assembled */
80 ASSERT (row == Rows [f] && row >= 0 && row < Work->n_row) ;
81
82 if (ep->rdeg == rdeg0)
83 {
84 /* ------------------------------------------------------ */
85 /* this is an old Lson - assemble just one row */
86 /* ------------------------------------------------------ */
87
88 /* flag the row as assembled from the Lson */
89 Rows [f] = EMPTY ;
90
91 nrows = ep->nrows ;
92 ncols = ep->ncols ;
93
94 p += UNITS (Int, ncols + nrows) ;
95 S = ((Entry *) p) + f ;
96
97 DEBUG6 (("Old LSON: " ID "\n", e)) ;
98 #ifndef NDEBUG
99 UMF_dump_element (Numeric, Work, e, FALSE) ;
100 #endif
101
102 ncolsleft = ep->ncolsleft ;
103
104 Frow = Fcblock + Frpos [row] ;
105 DEBUG6 (("LSON found (in scan2-row): " ID "\n", e)) ;
106
107 Row_degree [row] -= ncolsleft ;
108
109 if (ncols == ncolsleft)
110 {
111 /* -------------------------------------------------- */
112 /* no columns assembled out this Lson yet */
113 /* -------------------------------------------------- */
114
115 #pragma ivdep
116 for (j = 0 ; j < ncols ; j++)
117 {
118 col = Cols [j] ;
119 ASSERT (col >= 0 && col < Work->n_col) ;
120 #ifndef FIXQ
121 Col_degree [col] -- ;
122 #endif
123 /* Frow [Fcpos [col]] += *S ; */
124 ASSEMBLE (Frow [Fcpos [col]], *S) ;
125 S += nrows ;
126 }
127
128 }
129 else
130 {
131 /* -------------------------------------------------- */
132 /* some columns have been assembled out of this Lson */
133 /* -------------------------------------------------- */
134
135 #pragma ivdep
136 for (j = 0 ; j < ncols ; j++)
137 {
138 col = Cols [j] ;
139 if (col >= 0)
140 {
141 ASSERT (col < Work->n_col) ;
142 #ifndef FIXQ
143 Col_degree [col] -- ;
144 #endif
145 /* Frow [Fcpos [col]] += *S ; */
146 ASSEMBLE (Frow [Fcpos [col]], *S) ;
147 }
148 S += nrows ;
149 }
150
151 }
152 ep->nrowsleft-- ;
153 ASSERT (ep->nrowsleft > 0) ;
154 }
155 else
156 {
157 *tp2++ = *tp ; /* leave the tuple in the list */
158 }
159 }
160 Row_tlen [row] = tp2 - tp1 ;
161
162 #ifndef NDEBUG
163 DEBUG7 (("row assembled in scan2-row: " ID "\n", row)) ;
164 UMF_dump_rowcol (0, Numeric, Work, row, FALSE) ;
165 DEBUG7 (("Current frontal matrix: (scan 1b)\n")) ;
166 UMF_dump_current_front (Numeric, Work, TRUE) ;
167 #endif
168 }
169
170 /* ========================================================================== */
171 /* === col_assemble ========================================================= */
172 /* ========================================================================== */
173
col_assemble(Int col,NumericType * Numeric,WorkType * Work)174 PRIVATE void col_assemble
175 (
176 Int col,
177 NumericType *Numeric,
178 WorkType *Work
179 )
180 {
181
182 Entry *S, *Fcblock, *Fcol ;
183 Int tpi, e, *E, *Fcpos, *Frpos, *Row_degree, *Col_tuples, *Col_tlen, cdeg0,
184 f, nrows, ncols, *Rows, *Cols, row, nrowsleft, i ;
185 Tuple *tp, *tp1, *tp2, *tpend ;
186 Unit *Memory, *p ;
187 Element *ep ;
188
189 #if !defined (FIXQ) || !defined (NDEBUG)
190 Int *Col_degree ;
191 Col_degree = Numeric->Cperm ;
192 #endif
193
194 Col_tuples = Numeric->Lip ;
195 tpi = Col_tuples [col] ;
196 if (!tpi) return ;
197
198 Memory = Numeric->Memory ;
199 E = Work->E ;
200 Fcpos = Work->Fcpos ;
201 Frpos = Work->Frpos ;
202 Row_degree = Numeric->Rperm ;
203 Col_tlen = Numeric->Lilen ;
204 E = Work->E ;
205 Memory = Numeric->Memory ;
206 cdeg0 = Work->cdeg0 ;
207 Fcblock = Work->Fcblock ;
208
209 DEBUG6 (("SCAN2-col: " ID "\n", col)) ;
210 #ifndef NDEBUG
211 UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
212 #endif
213
214 ASSERT (NON_PIVOTAL_COL (col)) ;
215 tp = (Tuple *) (Memory + tpi) ;
216 tp1 = tp ;
217 tp2 = tp ;
218 tpend = tp + Col_tlen [col] ;
219 for ( ; tp < tpend ; tp++)
220 {
221 e = tp->e ;
222 ASSERT (e > 0 && e <= Work->nel) ;
223 if (!E [e]) continue ; /* element already deallocated */
224 f = tp->f ;
225 p = Memory + E [e] ;
226 ep = (Element *) p ;
227 p += UNITS (Element, 1) ;
228 Cols = (Int *) p ;
229
230 if (Cols [f] == EMPTY) continue ; /* col already assembled */
231 ASSERT (col == Cols [f] && col >= 0 && col < Work->n_col) ;
232
233 if (ep->cdeg == cdeg0)
234 {
235 /* ------------------------------------------------------ */
236 /* this is an old Uson - assemble just one col */
237 /* ------------------------------------------------------ */
238
239 /* flag the col as assembled from the Uson */
240 Cols [f] = EMPTY ;
241
242 nrows = ep->nrows ;
243 ncols = ep->ncols ;
244 Rows = Cols + ncols ;
245 p += UNITS (Int, ncols + nrows) ;
246 S = ((Entry *) p) + f * nrows ;
247
248 DEBUG6 (("Old USON: " ID "\n", e)) ;
249 #ifndef NDEBUG
250 UMF_dump_element (Numeric, Work, e, FALSE) ;
251 #endif
252
253 nrowsleft = ep->nrowsleft ;
254
255 Fcol = Fcblock + Fcpos [col] ;
256 DEBUG6 (("USON found (in scan2-col): " ID "\n", e)) ;
257 #ifndef FIXQ
258 Col_degree [col] -= nrowsleft ;
259 #endif
260 if (nrows == nrowsleft)
261 {
262 /* -------------------------------------------------- */
263 /* no rows assembled out of this Uson yet */
264 /* -------------------------------------------------- */
265
266 #pragma ivdep
267 for (i = 0 ; i < nrows ; i++)
268 {
269 row = Rows [i] ;
270 ASSERT (row >= 0 && row < Work->n_row) ;
271 Row_degree [row]-- ;
272 /* Fcol [Frpos [row]] += S [i] ; */
273 ASSEMBLE (Fcol [Frpos [row]], S [i]) ;
274 }
275 }
276 else
277 {
278 /* -------------------------------------------------- */
279 /* some rows have been assembled out of this Uson */
280 /* -------------------------------------------------- */
281
282 #pragma ivdep
283 for (i = 0 ; i < nrows ; i++)
284 {
285 row = Rows [i] ;
286 if (row >= 0)
287 {
288 ASSERT (row < Work->n_row) ;
289 Row_degree [row]-- ;
290 /* Fcol [Frpos [row]] += S [i] ; */
291 ASSEMBLE (Fcol [Frpos [row]], S [i]) ;
292 }
293 }
294 }
295 ep->ncolsleft-- ;
296 ASSERT (ep->ncolsleft > 0) ;
297 }
298 else
299 {
300 *tp2++ = *tp ; /* leave the tuple in the list */
301 }
302 }
303 Col_tlen [col] = tp2 - tp1 ;
304
305 #ifndef NDEBUG
306 DEBUG7 (("Column assembled in scan2-col: " ID "\n", col)) ;
307 UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
308 DEBUG7 (("Current frontal matrix: after scan2-col\n")) ;
309 UMF_dump_current_front (Numeric, Work, TRUE) ;
310 #endif
311
312 }
313
314
315 /* ========================================================================== */
316 /* === UMF_assemble / UMF_assemble_fixq ===================================== */
317 /* ========================================================================== */
318
319 #ifndef FIXQ
UMF_assemble(NumericType * Numeric,WorkType * Work)320 GLOBAL void UMF_assemble
321 #else
322 GLOBAL void UMF_assemble_fixq
323 #endif
324 (
325 NumericType *Numeric,
326 WorkType *Work
327 )
328 {
329 /* ---------------------------------------------------------------------- */
330 /* local variables */
331 /* ---------------------------------------------------------------------- */
332
333 Int e, i, row, col, i2, nrows, ncols, f, tpi, extcdeg, extrdeg, rdeg0,
334 cdeg0, son_list, next, nrows_to_assemble,
335 ncols_to_assemble, ngetrows, j, j2,
336 nrowsleft, /* number of rows remaining in S */
337 ncolsleft, /* number of columns remaining in S */
338 prior_Lson, prior_Uson, *E, *Cols, *Rows, *Wm, *Woo,
339 *Row_tuples, *Row_degree, *Row_tlen,
340 *Col_tuples, *Col_tlen ;
341 Unit *Memory, *p ;
342 Element *ep ;
343 Tuple *tp, *tp1, *tp2, *tpend ;
344 Entry
345 *S, /* a pointer into the contribution block of a son */
346 *Fcblock, /* current contribution block */
347 *Fcol ; /* a column of FC */
348 Int *Frpos,
349 *Fcpos,
350 fnrows, /* number of rows in contribution block in F */
351 fncols ; /* number of columns in contribution block in F */
352
353 #if !defined (FIXQ) || !defined (NDEBUG)
354 Int *Col_degree ;
355 #endif
356
357 #ifndef NDEBUG
358 Int n_row, n_col ;
359 n_row = Work->n_row ;
360 n_col = Work->n_col ;
361 DEBUG3 (("::Assemble SCANS 1-4\n")) ;
362 UMF_dump_current_front (Numeric, Work, TRUE) ;
363 #endif
364
365 #if !defined (FIXQ) || !defined (NDEBUG)
366 Col_degree = Numeric->Cperm ; /* not updated if FIXQ is true */
367 #endif
368
369 /* ---------------------------------------------------------------------- */
370 /* get parameters */
371 /* ---------------------------------------------------------------------- */
372
373 fncols = Work->fncols ;
374 fnrows = Work->fnrows ;
375 Fcpos = Work->Fcpos ;
376 Frpos = Work->Frpos ;
377 Row_degree = Numeric->Rperm ;
378 Row_tuples = Numeric->Uip ;
379 Row_tlen = Numeric->Uilen ;
380 Col_tuples = Numeric->Lip ;
381 Col_tlen = Numeric->Lilen ;
382 E = Work->E ;
383 Memory = Numeric->Memory ;
384 Wm = Work->Wm ;
385 Woo = Work->Woo ;
386 rdeg0 = Work->rdeg0 ;
387 cdeg0 = Work->cdeg0 ;
388
389 #ifndef NDEBUG
390 DEBUG6 (("============================================\n")) ;
391 DEBUG6 (("Degree update, assembly.\n")) ;
392 DEBUG6 (("pivot row pattern: fncols=" ID "\n", fncols)) ;
393 for (j = 0 ; j < fncols ; j++)
394 {
395 col = Work->Fcols [j] ;
396 DEBUG6 ((ID " ", col)) ;
397 ASSERT (Fcpos [col] == j * Work->fnr_curr) ;
398 ASSERT (NON_PIVOTAL_COL (col)) ;
399 }
400 ASSERT (Fcpos [Work->pivcol] >= 0) ;
401 DEBUG6 (("pivcol: " ID " pos " ID " fnr_curr " ID " fncols " ID "\n",
402 Work->pivcol, Fcpos [Work->pivcol], Work->fnr_curr, fncols)) ;
403 ASSERT (Fcpos [Work->pivcol] < fncols * Work->fnr_curr) ;
404 DEBUG6 (("\npivot col pattern: fnrows=" ID "\n", fnrows)) ;
405 for (i = 0 ; i < fnrows ; i++)
406 {
407 row = Work->Frows [i] ;
408 DEBUG6 ((ID " ", row)) ;
409 ASSERT (Frpos [row] == i) ;
410 ASSERT (NON_PIVOTAL_ROW (row)) ;
411 }
412 DEBUG6 (("\n")) ;
413 ASSERT (Frpos [Work->pivrow] >= 0) ;
414 ASSERT (Frpos [Work->pivrow] < fnrows) ;
415 ASSERT (Work->Flublock == (Entry *) (Numeric->Memory + E [0])) ;
416 ASSERT (Work->Fcblock == Work->Flublock + Work->nb *
417 (Work->nb + Work->fnr_curr + Work->fnc_curr)) ;
418 #endif
419
420 Fcblock = Work->Fcblock ;
421
422 /* ---------------------------------------------------------------------- */
423 /* determine the largest actual frontal matrix size (for Info only) */
424 /* ---------------------------------------------------------------------- */
425
426 ASSERT (fnrows == Work->fnrows_new + 1) ;
427 ASSERT (fncols == Work->fncols_new + 1) ;
428
429 Numeric->maxnrows = MAX (Numeric->maxnrows, fnrows) ;
430 Numeric->maxncols = MAX (Numeric->maxncols, fncols) ;
431
432 /* this is safe from integer overflow, since the current frontal matrix
433 * is already allocated. */
434 Numeric->maxfrsize = MAX (Numeric->maxfrsize, fnrows * fncols) ;
435
436 /* ---------------------------------------------------------------------- */
437 /* assemble from prior elements into the current frontal matrix */
438 /* ---------------------------------------------------------------------- */
439
440 DEBUG2 (("New assemble start [prior_element:" ID "\n", Work->prior_element)) ;
441
442 /* Currently no rows or columns are marked. No elements are scanned, */
443 /* that is, (ep->next == EMPTY) is true for all elements */
444
445 son_list = 0 ; /* start creating son_list [ */
446
447 /* ---------------------------------------------------------------------- */
448 /* determine if most recent element is Lson or Uson of current front */
449 /* ---------------------------------------------------------------------- */
450
451 if (!Work->do_extend)
452 {
453 prior_Uson = ( Work->pivcol_in_front && !Work->pivrow_in_front) ;
454 prior_Lson = (!Work->pivcol_in_front && Work->pivrow_in_front) ;
455 if (prior_Uson || prior_Lson)
456 {
457 e = Work->prior_element ;
458 if (e != EMPTY)
459 {
460 ASSERT (E [e]) ;
461 p = Memory + E [e] ;
462 ep = (Element *) p ;
463 ep->next = son_list ;
464 son_list = e ;
465 #ifndef NDEBUG
466 DEBUG2 (("e " ID " is Prior son " ID " " ID "\n",
467 e, prior_Uson, prior_Lson)) ;
468 UMF_dump_element (Numeric, Work, e, FALSE) ;
469 #endif
470 ASSERT (E [e]) ;
471 }
472 }
473 }
474 Work->prior_element = EMPTY ;
475
476 /* ---------------------------------------------------------------------- */
477 /* SCAN1-row: scan the element lists of each new row in the pivot col */
478 /* and compute the external column degree for each frontal */
479 /* ---------------------------------------------------------------------- */
480
481 for (i2 = Work->fscan_row ; i2 < fnrows ; i2++)
482 {
483 /* Get a row */
484 row = Work->NewRows [i2] ;
485 if (row < 0) row = FLIP (row) ;
486 ASSERT (row >= 0 && row < n_row) ;
487
488 DEBUG6 (("SCAN1-row: " ID "\n", row)) ;
489 #ifndef NDEBUG
490 UMF_dump_rowcol (0, Numeric, Work, row, FALSE) ;
491 #endif
492
493 ASSERT (NON_PIVOTAL_ROW (row)) ;
494 tpi = Row_tuples [row] ;
495 if (!tpi) continue ;
496 tp = (Tuple *) (Memory + tpi) ;
497 tp1 = tp ;
498 tp2 = tp ;
499 tpend = tp + Row_tlen [row] ;
500 for ( ; tp < tpend ; tp++)
501 {
502 e = tp->e ;
503 ASSERT (e > 0 && e <= Work->nel) ;
504 if (!E [e]) continue ; /* element already deallocated */
505 f = tp->f ;
506 p = Memory + E [e] ;
507 ep = (Element *) p ;
508 p += UNITS (Element, 1) ;
509 Rows = ((Int *) p) + ep->ncols ;
510 if (Rows [f] == EMPTY) continue ; /* row already assembled */
511 ASSERT (row == Rows [f]) ;
512
513 if (ep->cdeg < cdeg0)
514 {
515 /* first time seen in scan1-row */
516 ep->cdeg = ep->nrowsleft + cdeg0 ;
517 DEBUG6 (("e " ID " First seen: cdeg: " ID " ", e, ep->cdeg-cdeg0)) ;
518 ASSERT (ep->ncolsleft > 0 && ep->nrowsleft > 0) ;
519 }
520
521 ep->cdeg-- ; /* decrement external column degree */
522 DEBUG6 (("e " ID " New ext col deg: " ID "\n", e, ep->cdeg - cdeg0)) ;
523
524 /* this element is not yet in the new son list */
525 if (ep->cdeg == cdeg0 && ep->next == EMPTY)
526 {
527 /* A new LUson or Uson has been found */
528 ep->next = son_list ;
529 son_list = e ;
530 }
531
532 ASSERT (ep->cdeg >= cdeg0) ;
533 *tp2++ = *tp ; /* leave the tuple in the list */
534 }
535 Row_tlen [row] = tp2 - tp1 ;
536 }
537
538 /* ---------------------------------------------------------------------- */
539 /* SCAN1-col: scan the element lists of each new col in the pivot row */
540 /* and compute the external row degree for each frontal */
541 /* ---------------------------------------------------------------------- */
542
543 for (j2 = Work->fscan_col ; j2 < fncols ; j2++)
544 {
545 /* Get a column */
546 col = Work->NewCols [j2] ;
547 if (col < 0) col = FLIP (col) ;
548 ASSERT (col >= 0 && col < n_col) ;
549
550 DEBUG6 (("SCAN 1-col: " ID "\n", col)) ;
551 #ifndef NDEBUG
552 UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
553 #endif
554
555 ASSERT (NON_PIVOTAL_COL (col)) ;
556 tpi = Col_tuples [col] ;
557 if (!tpi) continue ;
558 tp = (Tuple *) (Memory + tpi) ;
559 tp1 = tp ;
560 tp2 = tp ;
561 tpend = tp + Col_tlen [col] ;
562 for ( ; tp < tpend ; tp++)
563 {
564 e = tp->e ;
565 ASSERT (e > 0 && e <= Work->nel) ;
566 if (!E [e]) continue ; /* element already deallocated */
567 f = tp->f ;
568 p = Memory + E [e] ;
569 ep = (Element *) p ;
570 p += UNITS (Element, 1) ;
571 Cols = (Int *) p ;
572 if (Cols [f] == EMPTY) continue ; /* column already assembled */
573 ASSERT (col == Cols [f]) ;
574
575 if (ep->rdeg < rdeg0)
576 {
577 /* first time seen in scan1-col */
578 ep->rdeg = ep->ncolsleft + rdeg0 ;
579 DEBUG6 (("e " ID " First seen: rdeg: " ID " ", e, ep->rdeg-rdeg0)) ;
580 ASSERT (ep->ncolsleft > 0 && ep->nrowsleft > 0) ;
581 }
582
583 ep->rdeg-- ; /* decrement external row degree */
584 DEBUG6 (("e " ID " New ext row degree: " ID "\n", e, ep->rdeg-rdeg0)) ;
585
586 if (ep->rdeg == rdeg0 && ep->next == EMPTY)
587 {
588 /* A new LUson or Lson has been found */
589 ep->next = son_list ;
590 son_list = e ;
591 }
592
593 ASSERT (ep->rdeg >= rdeg0) ;
594 *tp2++ = *tp ; /* leave the tuple in the list */
595 }
596 Col_tlen [col] = tp2 - tp1 ;
597 }
598
599 /* ---------------------------------------------------------------------- */
600 /* assemble new sons via full scans */
601 /* ---------------------------------------------------------------------- */
602
603 next = EMPTY ;
604
605 for (e = son_list ; e > 0 ; e = next)
606 {
607 ASSERT (e > 0 && e <= Work->nel && E [e]) ;
608 p = Memory + E [e] ;
609 DEBUG2 (("New son: " ID "\n", e)) ;
610 #ifndef NDEBUG
611 UMF_dump_element (Numeric, Work, e, FALSE) ;
612 #endif
613 GET_ELEMENT (ep, p, Cols, Rows, ncols, nrows, S) ;
614 nrowsleft = ep->nrowsleft ;
615 ncolsleft = ep->ncolsleft ;
616 next = ep->next ;
617 ep->next = EMPTY ;
618
619 extrdeg = (ep->rdeg < rdeg0) ? ncolsleft : (ep->rdeg - rdeg0) ;
620 extcdeg = (ep->cdeg < cdeg0) ? nrowsleft : (ep->cdeg - cdeg0) ;
621 ncols_to_assemble = ncolsleft - extrdeg ;
622 nrows_to_assemble = nrowsleft - extcdeg ;
623 DEBUG2 (("extrdeg " ID " extcdeg " ID "\n", extrdeg, extcdeg)) ;
624
625 if (extrdeg == 0 && extcdeg == 0)
626 {
627
628 /* -------------------------------------------------------------- */
629 /* this is an LUson - assemble an entire contribution block */
630 /* -------------------------------------------------------------- */
631
632 DEBUG6 (("LUson found: " ID "\n", e)) ;
633
634 if (nrows == nrowsleft)
635 {
636 /* ---------------------------------------------------------- */
637 /* no rows assembled out of this LUson yet */
638 /* ---------------------------------------------------------- */
639
640 /* compute the compressed column offset vector*/
641 /* [ use Wm [0..nrows-1] for offsets */
642 #pragma ivdep
643 for (i = 0 ; i < nrows ; i++)
644 {
645 row = Rows [i] ;
646 Row_degree [row] -= ncolsleft ;
647 Wm [i] = Frpos [row] ;
648 }
649
650 if (ncols == ncolsleft)
651 {
652 /* ------------------------------------------------------ */
653 /* no rows or cols assembled out of LUson yet */
654 /* ------------------------------------------------------ */
655
656 for (j = 0 ; j < ncols ; j++)
657 {
658 col = Cols [j] ;
659 #ifndef FIXQ
660 Col_degree [col] -= nrowsleft ;
661 #endif
662 Fcol = Fcblock + Fcpos [col] ;
663 #pragma ivdep
664 for (i = 0 ; i < nrows ; i++)
665 {
666 /* Fcol [Wm [i]] += S [i] ; */
667 ASSEMBLE (Fcol [Wm [i]], S [i]) ;
668 }
669 S += nrows ;
670 }
671
672
673 }
674 else
675 {
676 /* ------------------------------------------------------ */
677 /* only cols have been assembled out of LUson */
678 /* ------------------------------------------------------ */
679
680 for (j = 0 ; j < ncols ; j++)
681 {
682 col = Cols [j] ;
683 if (col >= 0)
684 {
685 #ifndef FIXQ
686 Col_degree [col] -= nrowsleft ;
687 #endif
688 Fcol = Fcblock + Fcpos [col] ;
689 #pragma ivdep
690 for (i = 0 ; i < nrows ; i++)
691 {
692 /* Fcol [Wm [i]] += S [i] ; */
693 ASSEMBLE (Fcol [Wm [i]], S [i]) ;
694 }
695 }
696 S += nrows ;
697 }
698
699 }
700 /* ] done using Wm [0..nrows-1] for offsets */
701 }
702 else
703 {
704 /* ---------------------------------------------------------- */
705 /* some rows have been assembled out of this LUson */
706 /* ---------------------------------------------------------- */
707
708 /* compute the compressed column offset vector*/
709 /* [ use Woo,Wm [0..nrowsleft-1] for offsets */
710 ngetrows = 0 ;
711 for (i = 0 ; i < nrows ; i++)
712 {
713 row = Rows [i] ;
714 if (row >= 0)
715 {
716 Row_degree [row] -= ncolsleft ;
717 Woo [ngetrows] = i ;
718 Wm [ngetrows++] = Frpos [row] ;
719 }
720 }
721 ASSERT (ngetrows == nrowsleft) ;
722
723 if (ncols == ncolsleft)
724 {
725 /* ------------------------------------------------------ */
726 /* only rows have been assembled out of this LUson */
727 /* ------------------------------------------------------ */
728
729 for (j = 0 ; j < ncols ; j++)
730 {
731 col = Cols [j] ;
732 #ifndef FIXQ
733 Col_degree [col] -= nrowsleft ;
734 #endif
735 Fcol = Fcblock + Fcpos [col] ;
736 #pragma ivdep
737 for (i = 0 ; i < nrowsleft ; i++)
738 {
739 /* Fcol [Wm [i]] += S [Woo [i]] ; */
740 ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
741 }
742 S += nrows ;
743 }
744
745 }
746 else
747 {
748 /* ------------------------------------------------------ */
749 /* both rows and columns have been assembled out of LUson */
750 /* ------------------------------------------------------ */
751
752 for (j = 0 ; j < ncols ; j++)
753 {
754 col = Cols [j] ;
755 if (col >= 0)
756 {
757 #ifndef FIXQ
758 Col_degree [col] -= nrowsleft ;
759 #endif
760 Fcol = Fcblock + Fcpos [col] ;
761 #pragma ivdep
762 for (i = 0 ; i < nrowsleft ; i++)
763 {
764 /* Fcol [Wm [i]] += S [Woo [i]] ; */
765 ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
766 }
767 }
768 S += nrows ;
769 }
770
771 }
772 /* ] done using Woo,Wm [0..nrowsleft-1] */
773 }
774
775 /* deallocate the element: remove from ordered list */
776 UMF_mem_free_tail_block (Numeric, E [e]) ;
777 E [e] = 0 ;
778
779 }
780 else if (extcdeg == 0)
781 {
782
783 /* -------------------------------------------------------------- */
784 /* this is a Uson - assemble all possible columns */
785 /* -------------------------------------------------------------- */
786
787 DEBUG6 (("New USON: " ID "\n", e)) ;
788 ASSERT (extrdeg > 0) ;
789
790 DEBUG6 (("New uson " ID " cols to do " ID "\n", e, ncols_to_assemble)) ;
791
792 if (ncols_to_assemble > 0)
793 {
794
795 Int skip = FALSE ;
796 if (ncols_to_assemble * 16 < ncols && nrows == 1)
797 {
798 /* this is a tall and thin frontal matrix consisting of
799 * only one column (most likely an original column). Do
800 * not assemble it. It cannot be the pivot column, since
801 * the pivot column element would be an LU son, not an Lson,
802 * of the current frontal matrix. */
803 ASSERT (nrowsleft == 1) ;
804 ASSERT (Rows [0] >= 0 && Rows [0] < Work->n_row) ;
805 skip = TRUE ;
806 Work->any_skip = TRUE ;
807 }
808
809 if (!skip)
810 {
811
812 if (nrows == nrowsleft)
813 {
814 /* -------------------------------------------------- */
815 /* no rows have been assembled out of this Uson yet */
816 /* -------------------------------------------------- */
817
818 /* compute the compressed column offset vector */
819 /* [ use Wm [0..nrows-1] for offsets */
820 #pragma ivdep
821 for (i = 0 ; i < nrows ; i++)
822 {
823 row = Rows [i] ;
824 ASSERT (row >= 0 && row < n_row) ;
825 Row_degree [row] -= ncols_to_assemble ;
826 Wm [i] = Frpos [row] ;
827 }
828
829 for (j = 0 ; j < ncols ; j++)
830 {
831 col = Cols [j] ;
832 if ((col >= 0) && (Fcpos [col] >= 0))
833 {
834 #ifndef FIXQ
835 Col_degree [col] -= nrowsleft ;
836 #endif
837 Fcol = Fcblock + Fcpos [col] ;
838 #pragma ivdep
839 for (i = 0 ; i < nrows ; i++)
840 {
841 /* Fcol [Wm [i]] += S [i] ; */
842 ASSEMBLE (Fcol [Wm [i]], S [i]) ;
843 }
844 /* flag the column as assembled from Uson */
845 Cols [j] = EMPTY ;
846 }
847 S += nrows ;
848 }
849
850
851 /* ] done using Wm [0..nrows-1] for offsets */
852 }
853 else
854 {
855 /* -------------------------------------------------- */
856 /* some rows have been assembled out of this Uson */
857 /* -------------------------------------------------- */
858
859 /* compute the compressed column offset vector*/
860 /* [ use Woo,Wm [0..nrows-1] for offsets */
861 ngetrows = 0 ;
862 for (i = 0 ; i < nrows ; i++)
863 {
864 row = Rows [i] ;
865 if (row >= 0)
866 {
867 Row_degree [row] -= ncols_to_assemble ;
868 ASSERT (row < n_row && Frpos [row] >= 0) ;
869 Woo [ngetrows] = i ;
870 Wm [ngetrows++] = Frpos [row] ;
871 }
872 }
873 ASSERT (ngetrows == nrowsleft) ;
874
875 for (j = 0 ; j < ncols ; j++)
876 {
877 col = Cols [j] ;
878 if ((col >= 0) && (Fcpos [col] >= 0))
879 {
880 #ifndef FIXQ
881 Col_degree [col] -= nrowsleft ;
882 #endif
883 Fcol = Fcblock + Fcpos [col] ;
884 #pragma ivdep
885 for (i = 0 ; i < nrowsleft ; i++)
886 {
887 /* Fcol [Wm [i]] += S [Woo [i]] ; */
888 ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
889 }
890 /* flag the column as assembled from Uson */
891 Cols [j] = EMPTY ;
892 }
893 S += nrows ;
894 }
895
896 /* ] done using Woo,Wm */
897 }
898 ep->ncolsleft = extrdeg ;
899 }
900 }
901
902 }
903 else
904 {
905
906 /* -------------------------------------------------------------- */
907 /* this is an Lson - assemble all possible rows */
908 /* -------------------------------------------------------------- */
909
910 DEBUG6 (("New LSON: " ID "\n", e)) ;
911 ASSERT (extrdeg == 0 && extcdeg > 0) ;
912
913 DEBUG6 (("New lson " ID " rows to do " ID "\n", e, nrows_to_assemble)) ;
914
915 if (nrows_to_assemble > 0)
916 {
917
918 Int skip = FALSE ;
919 if (nrows_to_assemble * 16 < nrows && ncols == 1)
920 {
921 /* this is a tall and thin frontal matrix consisting of
922 * only one column (most likely an original column). Do
923 * not assemble it. It cannot be the pivot column, since
924 * the pivot column element would be an LU son, not an Lson,
925 * of the current frontal matrix. */
926 ASSERT (ncolsleft == 1) ;
927 ASSERT (Cols [0] >= 0 && Cols [0] < Work->n_col) ;
928 Work->any_skip = TRUE ;
929 skip = TRUE ;
930 }
931
932 if (!skip)
933 {
934
935 /* compute the compressed column offset vector */
936 /* [ use Woo,Wm [0..nrows-1] for offsets */
937 ngetrows = 0 ;
938 for (i = 0 ; i < nrows ; i++)
939 {
940 row = Rows [i] ;
941 if ((row >= 0) && (Frpos [row] >= 0))
942 {
943 ASSERT (row < n_row) ;
944 Row_degree [row] -= ncolsleft ;
945 Woo [ngetrows] = i ;
946 Wm [ngetrows++] = Frpos [row] ;
947 /* flag the row as assembled from the Lson */
948 Rows [i] = EMPTY ;
949 }
950 }
951 ASSERT (nrowsleft - ngetrows == extcdeg) ;
952 ASSERT (ngetrows == nrows_to_assemble) ;
953
954 if (ncols == ncolsleft)
955 {
956 /* -------------------------------------------------- */
957 /* no columns assembled out this Lson yet */
958 /* -------------------------------------------------- */
959
960 for (j = 0 ; j < ncols ; j++)
961 {
962 col = Cols [j] ;
963 ASSERT (col >= 0 && col < n_col) ;
964 #ifndef FIXQ
965 Col_degree [col] -= nrows_to_assemble ;
966 #endif
967 Fcol = Fcblock + Fcpos [col] ;
968 #pragma ivdep
969 for (i = 0 ; i < nrows_to_assemble ; i++)
970 {
971 /* Fcol [Wm [i]] += S [Woo [i]] ; */
972 ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
973 }
974 S += nrows ;
975 }
976
977
978 }
979 else
980 {
981 /* -------------------------------------------------- */
982 /* some columns have been assembled out of this Lson */
983 /* -------------------------------------------------- */
984
985 for (j = 0 ; j < ncols ; j++)
986 {
987 col = Cols [j] ;
988 ASSERT (col < n_col) ;
989 if (col >= 0)
990 {
991 #ifndef FIXQ
992 Col_degree [col] -= nrows_to_assemble ;
993 #endif
994 Fcol = Fcblock + Fcpos [col] ;
995 #pragma ivdep
996 for (i = 0 ; i < nrows_to_assemble ; i++)
997 {
998 /* Fcol [Wm [i]] += S [Woo [i]] ; */
999 ASSEMBLE (Fcol [Wm [i]], S [Woo [i]]) ;
1000 }
1001 }
1002 S += nrows ;
1003 }
1004
1005 }
1006
1007 /* ] done using Woo,Wm */
1008
1009 ep->nrowsleft = extcdeg ;
1010 }
1011 }
1012 }
1013 }
1014
1015 /* Note that garbage collection, and build tuples */
1016 /* both destroy the son list. */
1017
1018 /* ] son_list now empty */
1019
1020 /* ---------------------------------------------------------------------- */
1021 /* If frontal matrix extended, assemble old L/Usons from new rows/cols */
1022 /* ---------------------------------------------------------------------- */
1023
1024 /* ---------------------------------------------------------------------- */
1025 /* SCAN2-row: assemble rows of old Lsons from the new rows */
1026 /* ---------------------------------------------------------------------- */
1027
1028 #ifndef NDEBUG
1029 DEBUG7 (("Current frontal matrix: (prior to scan2-row)\n")) ;
1030 UMF_dump_current_front (Numeric, Work, TRUE) ;
1031 #endif
1032
1033 /* rescan the pivot row */
1034 if (Work->any_skip)
1035 {
1036 row_assemble (Work->pivrow, Numeric, Work) ;
1037 }
1038
1039 if (Work->do_scan2row)
1040 {
1041 for (i2 = Work->fscan_row ; i2 < fnrows ; i2++)
1042 {
1043 /* Get a row */
1044 row = Work->NewRows [i2] ;
1045 if (row < 0) row = FLIP (row) ;
1046 ASSERT (row >= 0 && row < n_row) ;
1047 if (!(row == Work->pivrow && Work->any_skip))
1048 {
1049 /* assemble it */
1050 row_assemble (row, Numeric, Work) ;
1051 }
1052 }
1053 }
1054
1055 /* ---------------------------------------------------------------------- */
1056 /* SCAN2-col: assemble columns of old Usons from the new columns */
1057 /* ---------------------------------------------------------------------- */
1058
1059 #ifndef NDEBUG
1060 DEBUG7 (("Current frontal matrix: (prior to scan2-col)\n")) ;
1061 UMF_dump_current_front (Numeric, Work, TRUE) ;
1062 #endif
1063
1064 /* rescan the pivot col */
1065 if (Work->any_skip)
1066 {
1067 col_assemble (Work->pivcol, Numeric, Work) ;
1068 }
1069
1070 if (Work->do_scan2col)
1071 {
1072
1073 for (j2 = Work->fscan_col ; j2 < fncols ; j2++)
1074 {
1075 /* Get a column */
1076 col = Work->NewCols [j2] ;
1077 if (col < 0) col = FLIP (col) ;
1078 ASSERT (col >= 0 && col < n_col) ;
1079 if (!(col == Work->pivcol && Work->any_skip))
1080 {
1081 /* assemble it */
1082 col_assemble (col, Numeric, Work) ;
1083 }
1084 }
1085 }
1086
1087 /* ---------------------------------------------------------------------- */
1088 /* done. the remainder of this routine is used only when in debug mode */
1089 /* ---------------------------------------------------------------------- */
1090
1091 #ifndef NDEBUG
1092
1093 /* ---------------------------------------------------------------------- */
1094 /* when debugging: make sure the assembly did everything that it could */
1095 /* ---------------------------------------------------------------------- */
1096
1097 DEBUG3 (("::Assemble done\n")) ;
1098
1099 for (i2 = 0 ; i2 < fnrows ; i2++)
1100 {
1101 /* Get a row */
1102 row = Work->Frows [i2] ;
1103 ASSERT (row >= 0 && row < n_row) ;
1104
1105 DEBUG6 (("DEBUG SCAN 1: " ID "\n", row)) ;
1106 UMF_dump_rowcol (0, Numeric, Work, row, TRUE) ;
1107
1108 ASSERT (NON_PIVOTAL_ROW (row)) ;
1109 tpi = Row_tuples [row] ;
1110 if (!tpi) continue ;
1111 tp = (Tuple *) (Memory + tpi) ;
1112 tpend = tp + Row_tlen [row] ;
1113 for ( ; tp < tpend ; tp++)
1114 {
1115 e = tp->e ;
1116 ASSERT (e > 0 && e <= Work->nel) ;
1117 if (!E [e]) continue ; /* element already deallocated */
1118 f = tp->f ;
1119 p = Memory + E [e] ;
1120 ep = (Element *) p ;
1121 p += UNITS (Element, 1) ;
1122 Cols = (Int *) p ;
1123 Rows = ((Int *) p) + ep->ncols ;
1124 if (Rows [f] == EMPTY) continue ; /* row already assembled */
1125 ASSERT (row == Rows [f]) ;
1126 extrdeg = (ep->rdeg < rdeg0) ? ep->ncolsleft : (ep->rdeg - rdeg0) ;
1127 extcdeg = (ep->cdeg < cdeg0) ? ep->nrowsleft : (ep->cdeg - cdeg0) ;
1128 DEBUG6 ((
1129 "e " ID " After assembly ext row deg: " ID " ext col degree " ID "\n",
1130 e, extrdeg, extcdeg)) ;
1131
1132 if (Work->any_skip)
1133 {
1134 /* no Lsons in any row, except for very tall and thin ones */
1135 ASSERT (extrdeg >= 0) ;
1136 if (extrdeg == 0)
1137 {
1138 /* this is an unassemble Lson */
1139 ASSERT (ep->ncols == 1) ;
1140 ASSERT (ep->ncolsleft == 1) ;
1141 col = Cols [0] ;
1142 ASSERT (col != Work->pivcol) ;
1143 }
1144 }
1145 else
1146 {
1147 /* no Lsons in any row */
1148 ASSERT (extrdeg > 0) ;
1149 /* Uson external row degree is = number of cols left */
1150 ASSERT (IMPLIES (extcdeg == 0, extrdeg == ep->ncolsleft)) ;
1151 }
1152 }
1153 }
1154
1155 /* ---------------------------------------------------------------------- */
1156
1157 for (j2 = 0 ; j2 < fncols ; j2++)
1158 {
1159 /* Get a column */
1160 col = Work->Fcols [j2] ;
1161 ASSERT (col >= 0 && col < n_col) ;
1162
1163 DEBUG6 (("DEBUG SCAN 2: " ID "\n", col)) ;
1164 #ifndef FIXQ
1165 UMF_dump_rowcol (1, Numeric, Work, col, TRUE) ;
1166 #else
1167 UMF_dump_rowcol (1, Numeric, Work, col, FALSE) ;
1168 #endif
1169
1170 ASSERT (NON_PIVOTAL_COL (col)) ;
1171 tpi = Col_tuples [col] ;
1172 if (!tpi) continue ;
1173 tp = (Tuple *) (Memory + tpi) ;
1174 tpend = tp + Col_tlen [col] ;
1175 for ( ; tp < tpend ; tp++)
1176 {
1177 e = tp->e ;
1178 ASSERT (e > 0 && e <= Work->nel) ;
1179 if (!E [e]) continue ; /* element already deallocated */
1180 f = tp->f ;
1181 p = Memory + E [e] ;
1182 ep = (Element *) p ;
1183 p += UNITS (Element, 1) ;
1184 Cols = (Int *) p ;
1185 Rows = ((Int *) p) + ep->ncols ;
1186 if (Cols [f] == EMPTY) continue ; /* column already assembled */
1187 ASSERT (col == Cols [f]) ;
1188 extrdeg = (ep->rdeg < rdeg0) ? ep->ncolsleft : (ep->rdeg - rdeg0) ;
1189 extcdeg = (ep->cdeg < cdeg0) ? ep->nrowsleft : (ep->cdeg - cdeg0) ;
1190 DEBUG6 (("e " ID " After assembly ext col deg: " ID "\n", e, extcdeg)) ;
1191
1192 if (Work->any_skip)
1193 {
1194 /* no Usons in any column, except for very tall and thin ones */
1195 ASSERT (extcdeg >= 0) ;
1196 if (extcdeg == 0)
1197 {
1198 /* this is an unassemble Uson */
1199 ASSERT (ep->nrows == 1) ;
1200 ASSERT (ep->nrowsleft == 1) ;
1201 row = Rows [0] ;
1202 ASSERT (row != Work->pivrow) ;
1203 }
1204 }
1205 else
1206 {
1207 /* no Usons in any column */
1208 ASSERT (extcdeg > 0) ;
1209 /* Lson external column degree is = number of rows left */
1210 ASSERT (IMPLIES (extrdeg == 0, extcdeg == ep->nrowsleft)) ;
1211 }
1212 }
1213 }
1214 #endif /* NDEBUG */
1215 }
1216