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