1 /* ========================================================================== */
2 /* === UMF_local_search ===================================================== */
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 /*
12 Perform pivot search to find pivot row and pivot column.
13 The pivot column is selected from the candidate set. The candidate set
14 corresponds to a supercolumn from colamd or UMF_analyze. The pivot column
15 is then removed from that set. Constructs the pivot column pattern and
16 values. Called by umf_kernel. Returns UMFPACK_OK if successful, or
17 UMFPACK_WARNING_singular_matrix or UMFPACK_ERROR_different_pattern if not.
18 */
19
20 #include "umf_internal.h"
21 #include "umf_local_search.h"
22 #include "umf_row_search.h"
23 #include "umf_mem_free_tail_block.h"
24
25 /* relaxed amalgamation control parameters are fixed at compile time */
26 #define RELAX1 0.25
27 #define SYM_RELAX1 0.0
28 #define RELAX2 0.1
29 #define RELAX3 0.125
30
31 /* ========================================================================== */
32 /* === remove_candidate ===================================================== */
33 /* ========================================================================== */
34
35 /* Remove a column from the set of candidate pivot columns. */
36
remove_candidate(Int jj,WorkType * Work,SymbolicType * Symbolic)37 PRIVATE void remove_candidate (Int jj, WorkType *Work, SymbolicType *Symbolic)
38 {
39
40 #ifndef NDEBUG
41 Int j ;
42 DEBUGm2 (("pivot column Candidates before remove: nCand "ID" ncand "ID
43 " lo "ID" hi "ID" jj "ID"\n", Work->nCandidates, Work->ncand,
44 Work->lo, Work->hi, jj)) ;
45 for (j = 0 ; j < Work->nCandidates ; j++)
46 {
47 Int col = Work->Candidates [j] ;
48 DEBUGm2 ((ID" ", col));
49 ASSERT (col >= 0 && col < Work->n_col) ;
50 /* ASSERT (NON_PIVOTAL_COL (col)) ; */
51 ASSERT (col >= Work->lo && col <= Work->hi) ;
52 }
53 DEBUGm2 (("\n")) ;
54 #endif
55
56 if (Symbolic->fixQ)
57 {
58 DEBUGm2 (("FixQ\n")) ;
59 /* do not modify the column ordering */
60 ASSERT (Work->nCandidates == 1) ;
61 ASSERT (jj == 0) ;
62 if (Work->ncand > 1)
63 {
64 Work->Candidates [0] = Work->nextcand++ ;
65 }
66 else
67 {
68 Work->nCandidates = 0 ;
69 }
70 }
71 else
72 {
73 /* place the next candidate in the set */
74 if (Work->ncand > MAX_CANDIDATES)
75 {
76 Work->Candidates [jj] = Work->nextcand++ ;
77 }
78 else
79 {
80 ASSERT (Work->nCandidates == Work->ncand) ;
81 Work->Candidates [jj] = Work->Candidates [Work->ncand - 1] ;
82 Work->Candidates [Work->ncand - 1] = EMPTY ;
83 Work->nCandidates-- ;
84 }
85 }
86 Work->ncand-- ;
87
88 #ifndef NDEBUG
89 DEBUGm2 (("pivot column Candidates after remove: nCand "ID" ncand "ID
90 " lo "ID" hi "ID" jj "ID"\n", Work->nCandidates, Work->ncand, Work->lo,
91 Work->hi, jj)) ;
92 for (j = 0 ; j < Work->nCandidates ; j++)
93 {
94 Int col = Work->Candidates [j] ;
95 DEBUGm2 ((ID" ", col));
96 ASSERT (col >= 0 && col < Work->n_col) ;
97 /* ASSERT (NON_PIVOTAL_COL (col)) ; */
98 ASSERT (col >= Work->lo && col <= Work->hi) ;
99 }
100 DEBUGm2 (("\n")) ;
101 ASSERT (Work->ncand >= 0) ;
102 ASSERT (Work->nCandidates <= Work->ncand) ;
103 #endif
104 }
105
106 /* ========================================================================== */
107 /* === UMF_local_search ===================================================== */
108 /* ========================================================================== */
109
UMF_local_search(NumericType * Numeric,WorkType * Work,SymbolicType * Symbolic)110 GLOBAL Int UMF_local_search
111 (
112 NumericType *Numeric,
113 WorkType *Work,
114 SymbolicType *Symbolic
115 )
116 {
117 /* ---------------------------------------------------------------------- */
118 /* local variables */
119 /* ---------------------------------------------------------------------- */
120
121 double relax1 ;
122 Entry *Flblock, *Fublock, *Fs, *Fcblock, *C, *Wx, *Wy, *Fu, *Flublock,
123 *Flu ;
124 Int pos, nrows, *Cols, *Rows, e, f, status, max_cdeg, fnzeros, nb, j, col,
125 i, row, cdeg_in, rdeg [2][2], fnpiv, nothing [2], new_LUsize,
126 pivrow [2][2], pivcol [2], *Wp, *Fcpos, *Frpos, new_fnzeros, cdeg_out,
127 *Wm, *Wio, *Woi, *Woo, *Frows, *Fcols, fnrows, fncols, *E, deg, nr_in,
128 nc, thiscost, bestcost, nr_out, do_update, extra_cols, extra_rows,
129 extra_zeros, relaxed_front, do_extend, fnr_curr, fnc_curr, tpi,
130 *Col_tuples, *Col_degree, *Col_tlen, jj, jcand [2], freebie [2],
131 did_rowmerge, fnrows_new [2][2], fncols_new [2][2], search_pivcol_out,
132 *Diagonal_map, *Diagonal_imap, row2, col2 ;
133 Unit *Memory, *p ;
134 Tuple *tp, *tpend, *tp1, *tp2 ;
135 Element *ep ;
136
137 #ifndef NBLAS
138 Int blas_ok = TRUE ;
139 #else
140 #define blas_ok FALSE
141 #endif
142
143 #ifndef NDEBUG
144 Int debug_ok, n_row, n_col, *Row_degree ;
145 Row_degree = Numeric->Rperm ; /* for NON_PIVOTAL_ROW macro only */
146 #endif
147
148 /* ---------------------------------------------------------------------- */
149 /* get parameters */
150 /* ---------------------------------------------------------------------- */
151
152 Memory = Numeric->Memory ;
153 E = Work->E ;
154 Col_degree = Numeric->Cperm ;
155
156 Col_tuples = Numeric->Lip ;
157 Col_tlen = Numeric->Lilen ;
158
159 Wx = Work->Wx ;
160 Wy = Work->Wy ;
161 Wp = Work->Wp ;
162 Wm = Work->Wm ;
163 Woi = Work->Woi ;
164 Wio = Work->Wio ;
165 Woo = Work->Woo ;
166 Fcpos = Work->Fcpos ;
167 Frpos = Work->Frpos ;
168 Frows = Work->Frows ;
169 Fcols = Work->Fcols ;
170 fnrows = Work->fnrows ;
171 fncols = Work->fncols ;
172 nb = Work->nb ;
173 fnr_curr = Work->fnr_curr ;
174 fnc_curr = Work->fnc_curr ;
175 fnpiv = Work->fnpiv ;
176 nothing [0] = EMPTY ;
177 nothing [1] = EMPTY ;
178 relax1 = (Symbolic->prefer_diagonal) ? SYM_RELAX1 : RELAX1 ;
179 fnzeros = Work->fnzeros ;
180 new_fnzeros = fnzeros ;
181 jj = EMPTY ;
182
183 Fcblock = Work->Fcblock ; /* current contribution block */
184 Flblock = Work->Flblock ; /* current L block */
185 Fublock = Work->Fublock ; /* current U block */
186 Flublock = Work->Flublock ; /* current LU block */
187
188 /* The pivot column degree cannot exceed max_cdeg */
189 max_cdeg = Work->fnrows_max ;
190 ASSERT (Work->fnrows_max <= Symbolic->maxnrows) ;
191 ASSERT (Work->fncols_max <= Symbolic->maxncols) ;
192
193 if (fnrows == 0 && fncols == 0)
194 {
195 /* frontal matrix is empty */
196 Work->firstsuper = Work->ksuper ;
197 }
198
199 #ifndef NDEBUG
200 n_row = Work->n_row ;
201 n_col = Work->n_col ;
202 DEBUG2 (("\n========LOCAL SEARCH: current frontal matrix: ========= \n")) ;
203 UMF_dump_current_front (Numeric, Work, TRUE) ;
204 if (UMF_debug > 0 || MAX (n_row, n_col) < 1000)
205 {
206 for (i = 0 ; i < MAX (n_row, n_col) ; i++)
207 {
208 ASSERT (Wp [i] < 0) ;
209 }
210 }
211
212 DEBUGm2 ((ID" pivot column Candidates: lo "ID" hi "ID"\n",
213 Work->nCandidates, Work->lo, Work->hi)) ;
214 for (j = 0 ; j < Work->nCandidates ; j++)
215 {
216 col = Work->Candidates [j] ;
217 DEBUGm2 ((ID" ", col));
218 ASSERT (col >= 0 && col < n_col) ;
219 ASSERT (NON_PIVOTAL_COL (col)) ;
220 ASSERT (col >= Work->lo && col <= Work->hi) ;
221 }
222
223 DEBUGm2 (("\n")) ;
224 /* there are no 0-by-c or r-by-0 fronts, where c and r are > 0 */
225 /* a front is either 0-by-0, or r-by-c */
226 DEBUG2 (("\n\n::: "ID" : Npiv: "ID" + fnpiv "ID" = "ID". "
227 "size "ID"-by-"ID"\n", Work->frontid,
228 Work->npiv, Work->fnpiv, Work->npiv + Work->fnpiv, fnrows, fncols)) ;
229 ASSERT ((fnrows == 0 && fncols == 0) ||(fnrows != 0 && fncols != 0)) ;
230 #endif
231
232 /* ====================================================================== */
233 /* === PIVOT SEARCH ===================================================== */
234 /* ====================================================================== */
235
236 /* initialize */
237
238 pivcol [IN] = EMPTY ;
239 pivcol [OUT] = EMPTY ;
240
241 cdeg_in = Int_MAX ;
242 cdeg_out = Int_MAX ;
243
244 pivrow [IN][IN] = EMPTY ;
245 pivrow [IN][OUT] = EMPTY ;
246 pivrow [OUT][IN] = EMPTY ;
247 pivrow [OUT][OUT] = EMPTY ;
248
249 rdeg [IN][IN] = Int_MAX ;
250 rdeg [IN][OUT] = Int_MAX ;
251 rdeg [OUT][IN] = Int_MAX ;
252 rdeg [OUT][OUT] = Int_MAX ;
253
254 freebie [IN] = FALSE ;
255 freebie [OUT] = FALSE ;
256
257 Work->pivot_case = EMPTY ;
258 bestcost = EMPTY ;
259
260 nr_out = EMPTY ;
261 nr_in = EMPTY ;
262
263 jcand [IN] = EMPTY ;
264 jcand [OUT] = EMPTY ;
265
266 fnrows_new [IN][IN] = EMPTY ;
267 fnrows_new [IN][OUT] = EMPTY ;
268 fnrows_new [OUT][IN] = EMPTY ;
269 fnrows_new [OUT][OUT] = EMPTY ;
270
271 fncols_new [IN][IN] = EMPTY ;
272 fncols_new [IN][OUT] = EMPTY ;
273 fncols_new [OUT][IN] = EMPTY ;
274 fncols_new [OUT][OUT] = EMPTY ;
275
276 #ifndef NDEBUG
277 /* check Frpos */
278 DEBUG4 (("Check Frpos : fnrows "ID" col "ID" maxcdeg "ID"\n",
279 fnrows, pivcol [IN], max_cdeg)) ;
280 for (i = 0 ; i < fnrows ; i++)
281 {
282 row = Frows [i] ;
283 DEBUG4 ((" row: "ID"\n", row)) ;
284 ASSERT (row >= 0 && row < n_row) ;
285 ASSERT (Frpos [row] == i) ;
286 }
287 DEBUG4 (("All:\n")) ;
288 if (UMF_debug > 0 || n_row < 1000)
289 {
290 Int cnt = fnrows ;
291 for (row = 0 ; row < n_row ; row++)
292 {
293 if (Frpos [row] == EMPTY)
294 {
295 cnt++ ;
296 }
297 else
298 {
299 DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
300 }
301 }
302 ASSERT (cnt == n_row) ;
303 }
304 #endif
305
306 /* ---------------------------------------------------------------------- */
307 /* find shortest column in the front, and shortest column not in the */
308 /* front, from the candidate pivot column set */
309 /* ---------------------------------------------------------------------- */
310
311 /* If there are too many candidates, then only look at the first */
312 /* MAX_CANDIDATES of them. Otherwise, if there are O(n) candidates, */
313 /* this code could take O(n^2) time. */
314
315 /* ------------------------------------------------------------------ */
316 /* look in the candidate set for the best column */
317 /* ------------------------------------------------------------------ */
318
319 DEBUG2 (("Max candidates %d, Work->ncand "ID" jmax "ID"\n",
320 MAX_CANDIDATES, Work->ncand, Work->nCandidates)) ;
321 col = Work->Candidates [0] ;
322 ASSERT (Work->nCandidates > 0) ;
323 DEBUG3 (("Pivot column candidate: "ID" j = "ID"\n", col, j)) ;
324 ASSERT (col >= 0 && col < n_col) ;
325
326 /* there is no Col_degree if fixQ is true */
327 deg = Symbolic->fixQ ? EMPTY : Col_degree [col] ;
328
329 #ifndef NDEBUG
330 DEBUG3 (("Pivot column candidate: "ID" cost: "ID" Fcpos[col] "ID"\n",
331 col, deg, Fcpos [col])) ;
332 UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
333 if (Symbolic->fixQ)
334 {
335 DEBUG1 (("FIXQ: Candidates "ID" pivcol "ID" npiv "ID" fnpiv "ID
336 " ndiscard "ID "\n", Work->nCandidates, col, Work->npiv,
337 Work->fnpiv, Work->ndiscard)) ;
338 ASSERT (Work->nCandidates == 1) ;
339 ASSERT (col == Work->npiv + Work->fnpiv + Work->ndiscard) ;
340 }
341 #endif
342
343 if (Fcpos [col] >= 0)
344 {
345 /* best column in front, so far */
346 pivcol [IN] = col ;
347 cdeg_in = deg ; /* ignored, if fixQ is true */
348 jcand [IN] = 0 ;
349 }
350 else
351 {
352 /* best column not in front, so far */
353 pivcol [OUT] = col ;
354 cdeg_out = deg ; /* ignored, if fixQ is true */
355 jcand [OUT] = 0 ;
356 }
357
358 /* look at the rest of the candidates */
359 for (j = 1 ; j < Work->nCandidates ; j++)
360 {
361 col = Work->Candidates [j] ;
362
363 DEBUG3 (("Pivot col candidate: "ID" j = "ID"\n", col, j)) ;
364 ASSERT (col >= 0 && col < n_col) ;
365 ASSERT (!Symbolic->fixQ) ;
366 deg = Col_degree [col] ;
367 #ifndef NDEBUG
368 DEBUG3 (("Pivot col candidate: "ID" cost: "ID" Fcpos[col] "ID"\n",
369 col, deg, Fcpos [col])) ;
370 UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
371 #endif
372 if (Fcpos [col] >= 0)
373 {
374 #ifndef NDEBUG
375 Int fs ;
376 fs = Fcpos [col] / fnr_curr ;
377 ASSERT (fs >= 0 && fs < fncols) ;
378 #endif
379 if (deg < cdeg_in || (deg == cdeg_in && col < pivcol [IN]))
380 {
381 /* best column in front, so far */
382 pivcol [IN] = col ;
383 cdeg_in = deg ;
384 jcand [IN] = j ;
385 }
386 }
387 else
388 {
389 if (deg < cdeg_out || (deg == cdeg_out && col < pivcol [OUT]))
390 {
391 /* best column not in front, so far */
392 pivcol [OUT] = col ;
393 cdeg_out = deg ;
394 jcand [OUT] = j ;
395 }
396 }
397 }
398
399 DEBUG2 (("Pivcol in "ID" out "ID"\n", pivcol [IN], pivcol [OUT])) ;
400 ASSERT ((pivcol [IN] >= 0 && pivcol [IN] < n_col)
401 || (pivcol [OUT] >= 0 && pivcol [OUT] < n_col)) ;
402
403 cdeg_in = EMPTY ;
404 cdeg_out = EMPTY ;
405
406 /* ---------------------------------------------------------------------- */
407 /* construct candidate column in front, and search for pivot rows */
408 /* ---------------------------------------------------------------------- */
409
410 #ifndef NDEBUG
411 /* check Frpos */
412 DEBUG4 (("Prior to col update: fnrows "ID" col "ID" maxcdeg "ID"\n",
413 fnrows, pivcol [IN], max_cdeg)) ;
414 for (i = 0 ; i < fnrows ; i++)
415 {
416 row = Frows [i] ;
417 DEBUG4 ((" row: "ID"\n", row)) ;
418 ASSERT (row >= 0 && row < n_row) ;
419 ASSERT (Frpos [row] == i) ;
420 }
421 DEBUG4 (("All:\n")) ;
422 if (UMF_debug > 0 || n_row < 1000)
423 {
424 Int cnt = fnrows ;
425 for (row = 0 ; row < n_row ; row++)
426 {
427 if (Frpos [row] == EMPTY)
428 {
429 cnt++ ;
430 }
431 else
432 {
433 DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
434 }
435 }
436 ASSERT (cnt == n_row) ;
437 }
438 #endif
439
440 if (pivcol [IN] != EMPTY)
441 {
442
443 #ifndef NDEBUG
444 DEBUG2 (("col[IN] column "ID" in front at position = "ID"\n",
445 pivcol [IN], Fcpos [pivcol [IN]])) ;
446 UMF_dump_rowcol (1, Numeric, Work, pivcol [IN], !Symbolic->fixQ) ;
447 #endif
448
449 /* the only way we can have a pivcol[IN] is if the front is not empty */
450 ASSERT (fnrows > 0 && fncols > 0) ;
451
452 DEBUG4 (("Update pivot column:\n")) ;
453 Fs = Fcblock + Fcpos [pivcol [IN]] ;
454 Fu = Fublock + (Fcpos [pivcol [IN]] / fnr_curr) ;
455 Flu = Flublock + fnpiv * nb ;
456
457 /* ------------------------------------------------------------------ */
458 /* copy the pivot column from the U block into the LU block */
459 /* ------------------------------------------------------------------ */
460
461 /* This copy is permanent if the pivcol [IN] is chosen. */
462 for (i = 0 ; i < fnpiv ; i++)
463 {
464 Flu [i] = Fu [i*fnc_curr] ;
465 }
466
467 /* ------------------------------------------------------------------ */
468 /* update the pivot column in the LU block using a triangular solve */
469 /* ------------------------------------------------------------------ */
470
471 /* This work will be discarded if the pivcol [OUT] is chosen instead.
472 * It is permanent if the pivcol [IN] is chosen. */
473
474 if (fnpiv > 1)
475 {
476 /* solve Lx=b, where b = U (:,k), stored in the LU block */
477
478 #ifndef NBLAS
479 BLAS_TRSV (fnpiv, Flublock, Flu, nb) ;
480 #endif
481 if (!blas_ok)
482 {
483 /* use plain C code if no BLAS, or on integer overflow */
484 Entry *Flub = Flublock ;
485 for (j = 0 ; j < fnpiv ; j++)
486 {
487 Entry Fuj = Flu [j] ;
488 #pragma ivdep
489 for (i = j+1 ; i < fnpiv ; i++)
490 {
491 /* Flu [i] -= Flublock [i + j*nb] * Flu [j] ; */
492 MULT_SUB (Flu [i], Flub [i], Fuj) ;
493 }
494 Flub += nb ;
495 }
496 }
497 }
498
499 /* ------------------------------------------------------------------ */
500 /* copy the pivot column from the C block into Wy */
501 /* ------------------------------------------------------------------ */
502
503 for (i = 0 ; i < fnrows ; i++)
504 {
505 Wy [i] = Fs [i] ;
506 }
507
508 /* ------------------------------------------------------------------ */
509 /* update the pivot column of L using a matrix-vector multiply */
510 /* ------------------------------------------------------------------ */
511
512 /* this work will be discarded if the pivcol [OUT] is chosen instead */
513
514 #ifdef NBLAS
515 /* no BLAS available - use plain C code instead */
516 for (j = 0 ; j < fnpiv ; j++)
517 {
518 Entry Fuj, *Flub = Flblock + j * fnr_curr ;
519 Fuj = Flu [j] ;
520 if (IS_NONZERO (Fuj))
521 {
522 #pragma ivdep
523 for (i = 0 ; i < fnrows ; i++)
524 {
525 /* Wy [i] -= Flblock [i+j*fnr_curr] * Fuj ; */
526 MULT_SUB (Wy [i], Flub [i], Fuj) ;
527 }
528 }
529 /* Flblock += fnr_curr ; */
530 }
531 #else
532 /* Using 1-based notation:
533 * Wy (1:fnrows) -= Flblock (1:fnrows,1:fnpiv) * Flu (1:fnpiv) */
534 BLAS_GEMV (fnrows, fnpiv, Flblock, Flu, Wy, fnr_curr) ;
535 #endif
536
537 /* ------------------------------------------------------------------ */
538
539 #ifndef NDEBUG
540 DEBUG2 (("Wy after update: fnrows="ID"\n", fnrows)) ;
541 DEBUG4 ((" fnpiv="ID" \n", fnpiv)) ;
542 for (i = 0 ; i < fnrows ; i++)
543 {
544 DEBUG4 ((ID" "ID" "ID, i, Frows [i], Frpos [Frows [i]])) ;
545 EDEBUG4 (Wy [i]) ;
546 DEBUG4 (("\n")) ;
547 }
548 #endif
549
550 /* ------------------------------------------------------------------ */
551 /* construct the candidate column */
552 /* ------------------------------------------------------------------ */
553
554 cdeg_in = fnrows ;
555
556 #ifndef NDEBUG
557 /* check Frpos */
558 DEBUG4 (("After col update: fnrows "ID" col "ID" maxcdeg "ID"\n",
559 fnrows, pivcol [IN], max_cdeg)) ;
560 for (i = 0 ; i < fnrows ; i++)
561 {
562 row = Frows [i] ;
563 DEBUG4 ((" row: "ID"\n", row)) ;
564 ASSERT (row >= 0 && row < n_row) ;
565 ASSERT (Frpos [row] == i) ;
566 }
567 DEBUG4 (("All:\n")) ;
568 if (UMF_debug > 0 || n_row < 1000)
569 {
570 Int cnt = fnrows ;
571 for (row = 0 ; row < n_row ; row++)
572 {
573 if (Frpos [row] == EMPTY)
574 {
575 cnt++ ;
576 }
577 else
578 {
579 DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
580 }
581 }
582 ASSERT (cnt == n_row) ;
583 }
584 #endif
585
586 #ifndef NDEBUG
587 /* check Frpos */
588 DEBUG4 (("COL ASSEMBLE: cdeg "ID"\nREDUCE COL in "ID" max_cdeg "ID"\n",
589 cdeg_in, pivcol [IN], max_cdeg)) ;
590 for (i = 0 ; i < cdeg_in ; i++)
591 {
592 row = Frows [i] ;
593 ASSERT (row >= 0 && row < n_row) ;
594 ASSERT (Frpos [row] == i) ;
595 }
596 if (UMF_debug > 0 || n_row < 1000)
597 {
598 Int cnt = cdeg_in ;
599 for (row = 0 ; row < n_row ; row++)
600 {
601 if (Frpos [row] == EMPTY) cnt++ ;
602 }
603 ASSERT (cnt == n_row) ;
604 }
605 #endif
606
607 /* assemble column into Wy */
608
609 ASSERT (pivcol [IN] >= 0 && pivcol [IN] < n_col) ;
610 ASSERT (NON_PIVOTAL_COL (pivcol [IN])) ;
611
612 tpi = Col_tuples [pivcol [IN]] ;
613 if (tpi)
614 {
615 tp = (Tuple *) (Memory + tpi) ;
616 tp1 = tp ;
617 tp2 = tp ;
618 tpend = tp + Col_tlen [pivcol [IN]] ;
619 for ( ; tp < tpend ; tp++)
620 {
621 e = tp->e ;
622 ASSERT (e > 0 && e <= Work->nel) ;
623 if (!E [e]) continue ; /* element already deallocated */
624 f = tp->f ;
625 p = Memory + E [e] ;
626 ep = (Element *) p ;
627 p += UNITS (Element, 1) ;
628 Cols = (Int *) p ;
629 if (Cols [f] == EMPTY) continue ; /* column already assembled */
630 ASSERT (pivcol [IN] == Cols [f]) ;
631
632 Rows = Cols + ep->ncols ;
633 nrows = ep->nrows ;
634 p += UNITS (Int, ep->ncols + nrows) ;
635 C = ((Entry *) p) + f * nrows ;
636
637 for (i = 0 ; i < nrows ; i++)
638 {
639 row = Rows [i] ;
640 if (row >= 0) /* skip this if already gone from element */
641 {
642 ASSERT (row < n_row) ;
643 pos = Frpos [row] ;
644 if (pos < 0)
645 {
646 /* new entry in the pattern - save Frpos */
647 ASSERT (cdeg_in < n_row) ;
648 if (cdeg_in >= max_cdeg)
649 {
650 /* :: pattern change (cdeg in failure) :: */
651 DEBUGm4 (("cdeg_in failure\n")) ;
652 return (UMFPACK_ERROR_different_pattern) ;
653 }
654 Frpos [row] = cdeg_in ;
655 Frows [cdeg_in] = row ;
656 Wy [cdeg_in++] = C [i] ;
657 }
658 else
659 {
660 /* entry already in pattern - sum values in Wy */
661 /* Wy [pos] += C [i] ; */
662 ASSERT (pos < max_cdeg) ;
663 ASSEMBLE (Wy [pos], C [i]) ;
664 }
665 }
666 }
667 *tp2++ = *tp ; /* leave the tuple in the list */
668 }
669 Col_tlen [pivcol [IN]] = tp2 - tp1 ;
670 }
671
672 /* ------------------------------------------------------------------ */
673
674 #ifndef NDEBUG
675 /* check Frpos again */
676 DEBUG4 (("COL DONE: cdeg "ID"\nREDUCE COL in "ID" max_cdeg "ID"\n",
677 cdeg_in, pivcol [IN], max_cdeg)) ;
678 for (i = 0 ; i < cdeg_in ; i++)
679 {
680 row = Frows [i] ;
681 ASSERT (row >= 0 && row < n_row) ;
682 ASSERT (Frpos [row] == i) ;
683 }
684 if (UMF_debug > 0 || n_row < 1000)
685 {
686 Int cnt = cdeg_in ;
687 for (row = 0 ; row < n_row ; row++)
688 {
689 if (Frpos [row] == EMPTY) cnt++ ;
690 }
691 ASSERT (cnt == n_row) ;
692 }
693 #endif
694
695 #ifndef NDEBUG
696 DEBUG4 (("Reduced column: cdeg in "ID" fnrows_max "ID"\n",
697 cdeg_in, Work->fnrows_max)) ;
698 for (i = 0 ; i < cdeg_in ; i++)
699 {
700 DEBUG4 ((" "ID" "ID" "ID, i, Frows [i], Frpos [Frows [i]])) ;
701 EDEBUG4 (Wy [i]) ;
702 DEBUG4 (("\n")) ;
703 ASSERT (i == Frpos [Frows [i]]) ;
704 }
705 ASSERT (cdeg_in <= Work->fnrows_max) ;
706 #endif
707
708 /* ------------------------------------------------------------------ */
709 /* cdeg_in is now the exact degree of this column */
710 /* ------------------------------------------------------------------ */
711
712 nr_in = cdeg_in - fnrows ;
713
714 /* since there are no 0-by-x fronts, if there is a pivcol [IN] the */
715 /* front must have at least one row. */
716 ASSERT (cdeg_in > 0) ;
717
718 /* new degree of pivcol [IN], excluding current front is nr_in */
719 /* column expands by nr_in rows */
720
721 /* ------------------------------------------------------------------ */
722 /* search for two candidate pivot rows */
723 /* ------------------------------------------------------------------ */
724
725 /* for the IN_IN pivot row (if any), */
726 /* extend the pattern in place, using Fcols */
727 status = UMF_row_search (Numeric, Work, Symbolic,
728 fnrows, cdeg_in, Frows, Frpos, /* pattern of column to search */
729 pivrow [IN], rdeg [IN], Fcols, Wio, nothing, Wy,
730 pivcol [IN], freebie) ;
731 ASSERT (!freebie [IN] && !freebie [OUT]) ;
732
733 /* ------------------------------------------------------------------ */
734 /* fatal error if matrix pattern has changed since symbolic analysis */
735 /* ------------------------------------------------------------------ */
736
737 if (status == UMFPACK_ERROR_different_pattern)
738 {
739 /* :: pattern change (row search IN failure) :: */
740 DEBUGm4 (("row search IN failure\n")) ;
741 return (UMFPACK_ERROR_different_pattern) ;
742 }
743
744 /* ------------------------------------------------------------------ */
745 /* we now must have a structural pivot */
746 /* ------------------------------------------------------------------ */
747
748 /* Since the pivcol[IN] exists, there must be at least one row in the */
749 /* current frontal matrix, and so we must have found a structural */
750 /* pivot. The numerical value might be zero, of course. */
751
752 ASSERT (status != UMFPACK_WARNING_singular_matrix) ;
753
754 /* ------------------------------------------------------------------ */
755 /* evaluate IN_IN option */
756 /* ------------------------------------------------------------------ */
757
758 if (pivrow [IN][IN] != EMPTY)
759 {
760 /* The current front would become an (implicit) LUson.
761 * Both candidate pivot row and column are in the current front.
762 * Cost is how much the current front would expand */
763
764 /* pivrow[IN][IN] candidates are not found via row merge search */
765
766 ASSERT (fnrows >= 0 && fncols >= 0) ;
767
768 ASSERT (cdeg_in > 0) ;
769 nc = rdeg [IN][IN] - fncols ;
770
771 thiscost =
772 /* each column in front (except pivot column) grows by nr_in: */
773 (nr_in * (fncols - 1)) +
774 /* new columns not in old front: */
775 (nc * (cdeg_in - 1)) ;
776
777 /* no extra cost to relaxed amalgamation */
778
779 ASSERT (fnrows + nr_in == cdeg_in) ;
780 ASSERT (fncols + nc == rdeg [IN][IN]) ;
781
782 /* size of relaxed front (after pivot row column removed): */
783 fnrows_new [IN][IN] = (fnrows-1) + nr_in ;
784 fncols_new [IN][IN] = (fncols-1) + nc ;
785 /* relaxed_front = fnrows_new [IN][IN] * fncols_new [IN][IN] ; */
786
787 do_extend = TRUE ;
788
789 DEBUG2 (("Evaluating option IN-IN:\n")) ;
790 DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_in "ID" nc "ID"\n",
791 Work->fnzeros, fnpiv, nr_in, nc)) ;
792 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
793
794 /* determine if BLAS-3 update should be applied before extending. */
795 /* update if too many zero entries accumulate in the LU block */
796 fnzeros = Work->fnzeros + fnpiv * (nr_in + nc) ;
797
798 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
799
800 new_LUsize = (fnpiv+1) * (fnrows + nr_in + fncols + nc) ;
801
802 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
803
804 /* There are fnpiv pivots currently in the front. This one
805 * will be the (fnpiv+1)st pivot, if it is extended. */
806
807 /* RELAX2 parameter uses a double relop, but ignore NaN case: */
808 do_update = fnpiv > 0 &&
809 (((double) fnzeros) / ((double) new_LUsize)) > RELAX2 ;
810
811 DEBUG2 (("do_update "ID"\n", do_update))
812
813 DEBUG2 (("option IN IN : nr "ID" nc "ID" cost "ID"(0) relax "ID
814 "\n", nr_in, nc, thiscost, do_extend)) ;
815
816 /* this is the best option seen so far */
817 Work->pivot_case = IN_IN ;
818 bestcost = thiscost ;
819
820 /* do the amalgamation and extend the front */
821 Work->do_extend = do_extend ;
822 Work->do_update = do_update ;
823 new_fnzeros = fnzeros ;
824
825 }
826
827 /* ------------------------------------------------------------------ */
828 /* evaluate IN_OUT option */
829 /* ------------------------------------------------------------------ */
830
831 if (pivrow [IN][OUT] != EMPTY)
832 {
833 /* The current front would become a Uson of the new front.
834 * Candidate pivot column is in the current front, but the
835 * candidate pivot row is not. */
836
837 ASSERT (fnrows >= 0 && fncols > 0) ;
838 ASSERT (cdeg_in > 0) ;
839
840 /* must be at least one row outside the front */
841 /* (the pivrow [IN][OUT] itself) */
842 ASSERT (nr_in >= 1) ;
843
844 /* count columns not in current front */
845 nc = 0 ;
846 #ifndef NDEBUG
847 debug_ok = FALSE ;
848 #endif
849 for (i = 0 ; i < rdeg [IN][OUT] ; i++)
850 {
851 col = Wio [i] ;
852 DEBUG4 (("counting col "ID" Fcpos[] = "ID"\n", col,
853 Fcpos [col])) ;
854 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
855 if (Fcpos [col] < 0) nc++ ;
856 #ifndef NDEBUG
857 /* we must see the pivot column somewhere */
858 if (col == pivcol [IN])
859 {
860 ASSERT (Fcpos [col] >= 0) ;
861 debug_ok = TRUE ;
862 }
863 #endif
864 }
865 ASSERT (debug_ok) ;
866
867 thiscost =
868 /* each row in front grows by nc: */
869 (nc * fnrows) +
870 /* new rows not affected by front: */
871 ((nr_in-1) * (rdeg [IN][OUT]-1)) ;
872
873 /* check the cost of relaxed IN_OUT amalgamation */
874
875 extra_cols = ((fncols-1) + nc ) - (rdeg [IN][OUT] - 1) ;
876 ASSERT (extra_cols >= 0) ;
877 ASSERT (fncols + nc == extra_cols + rdeg [IN][OUT]) ;
878 extra_zeros = (nr_in-1) * extra_cols ; /* symbolic fill-in */
879
880 ASSERT (fnrows + nr_in == cdeg_in) ;
881 ASSERT (fncols + nc == rdeg [IN][OUT] + extra_cols) ;
882
883 /* size of relaxed front (after pivot column removed): */
884 fnrows_new [IN][OUT] = fnrows + (nr_in-1) ;
885 fncols_new [IN][OUT] = (fncols-1) + nc ;
886 relaxed_front = fnrows_new [IN][OUT] * fncols_new [IN][OUT] ;
887
888 /* do relaxed amalgamation if the extra zeros are no more */
889 /* than a fraction (default 0.25) of the relaxed front */
890 /* if relax = 0: no extra zeros allowed */
891 /* if relax = +inf: always amalgamate */
892
893 /* relax parameter uses a double relop, but ignore NaN case: */
894 if (extra_zeros == 0)
895 {
896 do_extend = TRUE ;
897 }
898 else
899 {
900 do_extend = ((double) extra_zeros) <
901 (relax1 * (double) relaxed_front) ;
902 }
903
904 if (do_extend)
905 {
906 /* count the cost of relaxed amalgamation */
907 thiscost += extra_zeros ;
908
909 DEBUG2 (("Evaluating option IN-OUT:\n")) ;
910 DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_in "ID" nc "ID"\n",
911 Work->fnzeros, fnpiv, nr_in, nc)) ;
912 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
913
914 /* determine if BLAS-3 update to be applied before extending. */
915 /* update if too many zero entries accumulate in the LU block */
916 fnzeros = Work->fnzeros + fnpiv * (nr_in + nc) ;
917
918 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
919
920 new_LUsize = (fnpiv+1) * (fnrows + nr_in + fncols + nc) ;
921
922 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
923
924 /* RELAX3 parameter uses a double relop, ignore NaN case: */
925 do_update = fnpiv > 0 &&
926 (((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
927 DEBUG2 (("do_update "ID"\n", do_update))
928
929 }
930 else
931 {
932 /* the current front would not be extended */
933 do_update = fnpiv > 0 ;
934 fnzeros = 0 ;
935 DEBUG2 (("IN-OUT do_update forced true: "ID"\n", do_update)) ;
936
937 /* The new front would be just big enough to hold the new
938 * pivot row and column. */
939 fnrows_new [IN][OUT] = cdeg_in - 1 ;
940 fncols_new [IN][OUT] = rdeg [IN][OUT] - 1 ;
941
942 }
943
944 DEBUG2 (("option IN OUT: nr "ID" nc "ID" cost "ID"("ID") relax "ID
945 "\n", nr_in, nc, thiscost, extra_zeros, do_extend)) ;
946
947 if (bestcost == EMPTY || thiscost < bestcost)
948 {
949 /* this is the best option seen so far */
950 Work->pivot_case = IN_OUT ;
951 bestcost = thiscost ;
952 Work->do_extend = do_extend ;
953 Work->do_update = do_update ;
954 new_fnzeros = fnzeros ;
955 }
956 }
957 }
958
959 /* ---------------------------------------------------------------------- */
960 /* construct candidate column not in front, and search for pivot rows */
961 /* ---------------------------------------------------------------------- */
962
963 search_pivcol_out = (bestcost != 0 && pivcol [OUT] != EMPTY) ;
964 if (Symbolic->prefer_diagonal)
965 {
966 search_pivcol_out = search_pivcol_out && (pivrow [IN][IN] == EMPTY) ;
967 }
968
969 if (search_pivcol_out)
970 {
971
972 #ifndef NDEBUG
973 DEBUG2 (("out_col column "ID" NOT in front at position = "ID"\n",
974 pivcol [OUT], Fcpos [pivcol [OUT]])) ;
975 UMF_dump_rowcol (1, Numeric, Work, pivcol [OUT], !Symbolic->fixQ) ;
976 DEBUG2 (("fncols "ID" fncols_max "ID"\n", fncols, Work->fncols_max)) ;
977 ASSERT (fncols < Work->fncols_max) ;
978 #endif
979
980 /* Use Wx as temporary workspace to construct the pivcol [OUT] */
981
982
983 /* ------------------------------------------------------------------ */
984 /* construct the candidate column (currently not in the front) */
985 /* ------------------------------------------------------------------ */
986
987 /* Construct the column in Wx, Wm, using Wp for the positions: */
988 /* Wm [0..cdeg_out-1] list of row indices in the column */
989 /* Wx [0..cdeg_out-1] list of corresponding numerical values */
990 /* Wp [0..n-1] starts as all negative, and ends that way too. */
991
992 cdeg_out = 0 ;
993
994 #ifndef NDEBUG
995 /* check Wp */
996 DEBUG4 (("COL ASSEMBLE: cdeg 0\nREDUCE COL out "ID"\n", pivcol [OUT])) ;
997 if (UMF_debug > 0 || MAX (n_row, n_col) < 1000)
998 {
999 for (i = 0 ; i < MAX (n_row, n_col) ; i++)
1000 {
1001 ASSERT (Wp [i] < 0) ;
1002 }
1003 }
1004 DEBUG4 (("max_cdeg: "ID"\n", max_cdeg)) ;
1005 #endif
1006
1007 ASSERT (pivcol [OUT] >= 0 && pivcol [OUT] < n_col) ;
1008 ASSERT (NON_PIVOTAL_COL (pivcol [OUT])) ;
1009
1010 tpi = Col_tuples [pivcol [OUT]] ;
1011 if (tpi)
1012 {
1013 tp = (Tuple *) (Memory + tpi) ;
1014 tp1 = tp ;
1015 tp2 = tp ;
1016 tpend = tp + Col_tlen [pivcol [OUT]] ;
1017 for ( ; tp < tpend ; tp++)
1018 {
1019 e = tp->e ;
1020 ASSERT (e > 0 && e <= Work->nel) ;
1021 if (!E [e]) continue ; /* element already deallocated */
1022 f = tp->f ;
1023 p = Memory + E [e] ;
1024 ep = (Element *) p ;
1025 p += UNITS (Element, 1) ;
1026 Cols = (Int *) p ;
1027 if (Cols [f] == EMPTY) continue ; /* column already assembled */
1028 ASSERT (pivcol [OUT] == Cols [f]) ;
1029
1030 Rows = Cols + ep->ncols ;
1031 nrows = ep->nrows ;
1032 p += UNITS (Int, ep->ncols + nrows) ;
1033 C = ((Entry *) p) + f * nrows ;
1034
1035 for (i = 0 ; i < nrows ; i++)
1036 {
1037 row = Rows [i] ;
1038 if (row >= 0) /* skip this if already gone from element */
1039 {
1040 ASSERT (row < n_row) ;
1041 pos = Wp [row] ;
1042 if (pos < 0)
1043 {
1044 /* new entry in the pattern - save Wp */
1045 ASSERT (cdeg_out < n_row) ;
1046 if (cdeg_out >= max_cdeg)
1047 {
1048 /* :: pattern change (cdeg out failure) :: */
1049 DEBUGm4 (("cdeg out failure\n")) ;
1050 return (UMFPACK_ERROR_different_pattern) ;
1051 }
1052 Wp [row] = cdeg_out ;
1053 Wm [cdeg_out] = row ;
1054 Wx [cdeg_out++] = C [i] ;
1055 }
1056 else
1057 {
1058 /* entry already in pattern - sum the values */
1059 /* Wx [pos] += C [i] ; */
1060 ASSEMBLE (Wx [pos], C [i]) ;
1061 }
1062 }
1063 }
1064 *tp2++ = *tp ; /* leave the tuple in the list */
1065 }
1066 Col_tlen [pivcol [OUT]] = tp2 - tp1 ;
1067 }
1068
1069 /* ------------------------------------------------------------------ */
1070
1071 #ifndef NDEBUG
1072 DEBUG4 (("Reduced column: cdeg out "ID"\n", cdeg_out)) ;
1073 for (i = 0 ; i < cdeg_out ; i++)
1074 {
1075 DEBUG4 ((" "ID" "ID" "ID, i, Wm [i], Wp [Wm [i]])) ;
1076 EDEBUG4 (Wx [i]) ;
1077 DEBUG4 (("\n")) ;
1078 ASSERT (i == Wp [Wm [i]]) ;
1079 }
1080 #endif
1081
1082 /* ------------------------------------------------------------------ */
1083 /* new degree of pivcol [OUT] is cdeg_out */
1084 /* ------------------------------------------------------------------ */
1085
1086 /* search for two candidate pivot rows */
1087 status = UMF_row_search (Numeric, Work, Symbolic,
1088 0, cdeg_out, Wm, Wp, /* pattern of column to search */
1089 pivrow [OUT], rdeg [OUT], Woi, Woo, pivrow [IN], Wx,
1090 pivcol [OUT], freebie) ;
1091
1092 /* ------------------------------------------------------------------ */
1093 /* fatal error if matrix pattern has changed since symbolic analysis */
1094 /* ------------------------------------------------------------------ */
1095
1096 if (status == UMFPACK_ERROR_different_pattern)
1097 {
1098 /* :: pattern change detected in umf_local_search :: */
1099 return (UMFPACK_ERROR_different_pattern) ;
1100 }
1101
1102 /* ------------------------------------------------------------------ */
1103 /* Clear Wp */
1104 /* ------------------------------------------------------------------ */
1105
1106 for (i = 0 ; i < cdeg_out ; i++)
1107 {
1108 Wp [Wm [i]] = EMPTY ; /* clear Wp */
1109 }
1110
1111 /* ------------------------------------------------------------------ */
1112 /* check for rectangular, singular matrix */
1113 /* ------------------------------------------------------------------ */
1114
1115 if (status == UMFPACK_WARNING_singular_matrix)
1116 {
1117 /* Pivot column is empty, and row-merge set is empty too. The
1118 * matrix is structurally singular. The current frontal matrix must
1119 * be empty, too. It it weren't, and pivcol [OUT] exists, then
1120 * there would be at least one row that could be selected. Since
1121 * the current front is empty, pivcol [IN] must also be EMPTY.
1122 */
1123
1124 DEBUGm4 (("Note: pivcol [OUT]: "ID" discard\n", pivcol [OUT])) ;
1125 ASSERT ((Work->fnrows == 0 && Work->fncols == 0)) ;
1126 ASSERT (pivcol [IN] == EMPTY) ;
1127
1128 /* remove the failed pivcol [OUT] from candidate set */
1129 ASSERT (pivcol [OUT] == Work->Candidates [jcand [OUT]]) ;
1130 remove_candidate (jcand [OUT], Work, Symbolic) ;
1131 Work->ndiscard++ ;
1132
1133 /* delete all of the tuples, and all contributions to this column */
1134 DEBUG1 (("Prune tuples of dead outcol: "ID"\n", pivcol [OUT])) ;
1135 Col_tlen [pivcol [OUT]] = 0 ;
1136 UMF_mem_free_tail_block (Numeric, Col_tuples [pivcol [OUT]]) ;
1137 Col_tuples [pivcol [OUT]] = 0 ;
1138
1139 /* no pivot found at all */
1140 return (UMFPACK_WARNING_singular_matrix) ;
1141 }
1142
1143 /* ------------------------------------------------------------------ */
1144
1145 if (freebie [IN])
1146 {
1147 /* the "in" row is the same as the "in" row for the "in" column */
1148 Woi = Fcols ;
1149 rdeg [OUT][IN] = rdeg [IN][IN] ;
1150 DEBUG4 (("Freebie in, row "ID"\n", pivrow [IN][IN])) ;
1151 }
1152
1153 if (freebie [OUT])
1154 {
1155 /* the "out" row is the same as the "out" row for the "in" column */
1156 Woo = Wio ;
1157 rdeg [OUT][OUT] = rdeg [IN][OUT] ;
1158 DEBUG4 (("Freebie out, row "ID"\n", pivrow [IN][OUT])) ;
1159 }
1160
1161 /* ------------------------------------------------------------------ */
1162 /* evaluate OUT_IN option */
1163 /* ------------------------------------------------------------------ */
1164
1165 if (pivrow [OUT][IN] != EMPTY)
1166 {
1167 /* The current front would become an Lson of the new front.
1168 * The candidate pivot row is in the current front, but the
1169 * candidate pivot column is not. */
1170
1171 ASSERT (fnrows > 0 && fncols >= 0) ;
1172
1173 did_rowmerge = (cdeg_out == 0) ;
1174 if (did_rowmerge)
1175 {
1176 /* pivrow [OUT][IN] was found via row merge search */
1177 /* it is not (yet) in the pivot column pattern (add it now) */
1178 DEBUGm4 (("did row merge OUT col, IN row\n")) ;
1179 Wm [0] = pivrow [OUT][IN] ;
1180 CLEAR (Wx [0]) ;
1181 cdeg_out = 1 ;
1182 ASSERT (nr_out == EMPTY) ;
1183 }
1184
1185 nc = rdeg [OUT][IN] - fncols ;
1186 ASSERT (nc >= 1) ;
1187
1188 /* count rows not in current front */
1189 nr_out = 0 ;
1190 #ifndef NDEBUG
1191 debug_ok = FALSE ;
1192 #endif
1193 for (i = 0 ; i < cdeg_out ; i++)
1194 {
1195 row = Wm [i] ;
1196 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1197 if (Frpos [row] < 0 || Frpos [row] >= fnrows) nr_out++ ;
1198 #ifndef NDEBUG
1199 /* we must see the pivot row somewhere */
1200 if (row == pivrow [OUT][IN])
1201 {
1202 ASSERT (Frpos [row] >= 0) ;
1203 debug_ok = TRUE ;
1204 }
1205 #endif
1206 }
1207 ASSERT (debug_ok) ;
1208
1209 thiscost =
1210 /* each column in front grows by nr_out: */
1211 (nr_out * fncols) +
1212 /* new cols not affected by front: */
1213 ((nc-1) * (cdeg_out-1)) ;
1214
1215 /* check the cost of relaxed OUT_IN amalgamation */
1216
1217 extra_rows = ((fnrows-1) + nr_out) - (cdeg_out - 1) ;
1218 ASSERT (extra_rows >= 0) ;
1219 ASSERT (fnrows + nr_out == extra_rows + cdeg_out) ;
1220 extra_zeros = (nc-1) * extra_rows ; /* symbolic fill-in */
1221
1222 ASSERT (fnrows + nr_out == cdeg_out + extra_rows) ;
1223 ASSERT (fncols + nc == rdeg [OUT][IN]) ;
1224
1225 /* size of relaxed front (after pivot row removed): */
1226 fnrows_new [OUT][IN] = (fnrows-1) + nr_out ;
1227 fncols_new [OUT][IN] = fncols + (nc-1) ;
1228 relaxed_front = fnrows_new [OUT][IN] * fncols_new [OUT][IN] ;
1229
1230 /* do relaxed amalgamation if the extra zeros are no more */
1231 /* than a fraction (default 0.25) of the relaxed front */
1232 /* if relax = 0: no extra zeros allowed */
1233 /* if relax = +inf: always amalgamate */
1234 if (did_rowmerge)
1235 {
1236 do_extend = FALSE ;
1237 }
1238 else
1239 {
1240 /* relax parameter uses a double relop, but ignore NaN case: */
1241 if (extra_zeros == 0)
1242 {
1243 do_extend = TRUE ;
1244 }
1245 else
1246 {
1247 do_extend = ((double) extra_zeros) <
1248 (relax1 * (double) relaxed_front) ;
1249 }
1250 }
1251
1252 if (do_extend)
1253 {
1254 /* count the cost of relaxed amalgamation */
1255 thiscost += extra_zeros ;
1256
1257 DEBUG2 (("Evaluating option OUT-IN:\n")) ;
1258 DEBUG2 ((" Work->fnzeros "ID" fnpiv "ID" nr_out "ID" nc "ID"\n",
1259 Work->fnzeros, fnpiv, nr_out, nc)) ;
1260 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
1261
1262 /* determine if BLAS-3 update to be applied before extending. */
1263 /* update if too many zero entries accumulate in the LU block */
1264 fnzeros = Work->fnzeros + fnpiv * (nr_out + nc) ;
1265
1266 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
1267
1268 new_LUsize = (fnpiv+1) * (fnrows + nr_out + fncols + nc) ;
1269
1270 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
1271
1272 /* RELAX3 parameter uses a double relop, ignore NaN case: */
1273 do_update = fnpiv > 0 &&
1274 (((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
1275 DEBUG2 (("do_update "ID"\n", do_update))
1276 }
1277 else
1278 {
1279 /* the current front would not be extended */
1280 do_update = fnpiv > 0 ;
1281 fnzeros = 0 ;
1282 DEBUG2 (("OUT-IN do_update forced true: "ID"\n", do_update)) ;
1283
1284 /* The new front would be just big enough to hold the new
1285 * pivot row and column. */
1286 fnrows_new [OUT][IN] = cdeg_out - 1 ;
1287 fncols_new [OUT][IN] = rdeg [OUT][IN] - 1 ;
1288 }
1289
1290 DEBUG2 (("option OUT IN : nr "ID" nc "ID" cost "ID"("ID") relax "ID
1291 "\n", nr_out, nc, thiscost, extra_zeros, do_extend)) ;
1292
1293 if (bestcost == EMPTY || thiscost < bestcost)
1294 {
1295 /* this is the best option seen so far */
1296 Work->pivot_case = OUT_IN ;
1297 bestcost = thiscost ;
1298 Work->do_extend = do_extend ;
1299 Work->do_update = do_update ;
1300 new_fnzeros = fnzeros ;
1301 }
1302 }
1303
1304 /* ------------------------------------------------------------------ */
1305 /* evaluate OUT_OUT option */
1306 /* ------------------------------------------------------------------ */
1307
1308 if (pivrow [OUT][OUT] != EMPTY)
1309 {
1310 /* Neither the candidate pivot row nor the candidate pivot column
1311 * are in the current front. */
1312
1313 ASSERT (fnrows >= 0 && fncols >= 0) ;
1314
1315 did_rowmerge = (cdeg_out == 0) ;
1316 if (did_rowmerge)
1317 {
1318 /* pivrow [OUT][OUT] was found via row merge search */
1319 /* it is not (yet) in the pivot column pattern (add it now) */
1320 DEBUGm4 (("did row merge OUT col, OUT row\n")) ;
1321 Wm [0] = pivrow [OUT][OUT] ;
1322 CLEAR (Wx [0]) ;
1323 cdeg_out = 1 ;
1324 ASSERT (nr_out == EMPTY) ;
1325 nr_out = 1 ;
1326 }
1327
1328 if (fnrows == 0 && fncols == 0)
1329 {
1330 /* the current front is completely empty */
1331 ASSERT (fnpiv == 0) ;
1332 nc = rdeg [OUT][OUT] ;
1333 extra_cols = 0 ;
1334 nr_out = cdeg_out ;
1335 extra_rows = 0 ;
1336 extra_zeros = 0 ;
1337
1338 thiscost = (nc-1) * (cdeg_out-1) ; /* new columns only */
1339
1340 /* size of new front: */
1341 fnrows_new [OUT][OUT] = nr_out-1 ;
1342 fncols_new [OUT][OUT] = nc-1 ;
1343 relaxed_front = fnrows_new [OUT][OUT] * fncols_new [OUT][OUT] ;
1344 }
1345 else
1346 {
1347
1348 /* count rows not in current front */
1349 if (nr_out == EMPTY)
1350 {
1351 nr_out = 0 ;
1352 #ifndef NDEBUG
1353 debug_ok = FALSE ;
1354 #endif
1355 for (i = 0 ; i < cdeg_out ; i++)
1356 {
1357 row = Wm [i] ;
1358 ASSERT (row >= 0 && row < n_row) ;
1359 ASSERT (NON_PIVOTAL_ROW (row)) ;
1360 if (Frpos [row] < 0 || Frpos [row] >= fnrows) nr_out++ ;
1361 #ifndef NDEBUG
1362 /* we must see the pivot row somewhere */
1363 if (row == pivrow [OUT][OUT])
1364 {
1365 ASSERT (Frpos [row] < 0 || Frpos [row] >= fnrows) ;
1366 debug_ok = TRUE ;
1367 }
1368 #endif
1369 }
1370 ASSERT (debug_ok) ;
1371 }
1372
1373 /* count columns not in current front */
1374 nc = 0 ;
1375 #ifndef NDEBUG
1376 debug_ok = FALSE ;
1377 #endif
1378 for (i = 0 ; i < rdeg [OUT][OUT] ; i++)
1379 {
1380 col = Woo [i] ;
1381 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1382 if (Fcpos [col] < 0) nc++ ;
1383 #ifndef NDEBUG
1384 /* we must see the pivot column somewhere */
1385 if (col == pivcol [OUT])
1386 {
1387 ASSERT (Fcpos [col] < 0) ;
1388 debug_ok = TRUE ;
1389 }
1390 #endif
1391 }
1392 ASSERT (debug_ok) ;
1393
1394 extra_cols = (fncols + (nc-1)) - (rdeg [OUT][OUT] - 1) ;
1395 extra_rows = (fnrows + (nr_out-1)) - (cdeg_out - 1) ;
1396 ASSERT (extra_rows >= 0) ;
1397 ASSERT (extra_cols >= 0) ;
1398 extra_zeros = ((nc-1) * extra_rows) + ((nr_out-1) * extra_cols);
1399
1400 ASSERT (fnrows + nr_out == cdeg_out + extra_rows) ;
1401 ASSERT (fncols + nc == rdeg [OUT][OUT] + extra_cols) ;
1402
1403 thiscost =
1404 /* new columns: */
1405 ((nc-1) * (cdeg_out-1)) +
1406 /* old columns in front grow by nr_out-1: */
1407 ((nr_out-1) * (fncols - extra_cols)) ;
1408
1409 /* size of relaxed front: */
1410 fnrows_new [OUT][OUT] = fnrows + (nr_out-1) ;
1411 fncols_new [OUT][OUT] = fncols + (nc-1) ;
1412 relaxed_front = fnrows_new [OUT][OUT] * fncols_new [OUT][OUT] ;
1413
1414 }
1415
1416 /* do relaxed amalgamation if the extra zeros are no more */
1417 /* than a fraction (default 0.25) of the relaxed front */
1418 /* if relax = 0: no extra zeros allowed */
1419 /* if relax = +inf: always amalgamate */
1420 if (did_rowmerge)
1421 {
1422 do_extend = FALSE ;
1423 }
1424 else
1425 {
1426 /* relax parameter uses a double relop, but ignore NaN case: */
1427 if (extra_zeros == 0)
1428 {
1429 do_extend = TRUE ;
1430 }
1431 else
1432 {
1433 do_extend = ((double) extra_zeros) <
1434 (relax1 * (double) relaxed_front) ;
1435 }
1436 }
1437
1438 if (do_extend)
1439 {
1440 /* count the cost of relaxed amalgamation */
1441 thiscost += extra_zeros ;
1442
1443 DEBUG2 (("Evaluating option OUT-OUT:\n")) ;
1444 DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_out "ID" nc "ID"\n",
1445 Work->fnzeros, fnpiv, nr_out, nc)) ;
1446 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
1447
1448 /* determine if BLAS-3 update to be applied before extending. */
1449 /* update if too many zero entries accumulate in the LU block */
1450 fnzeros = Work->fnzeros + fnpiv * (nr_out + nc) ;
1451
1452 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
1453
1454 new_LUsize = (fnpiv+1) * (fnrows + nr_out + fncols + nc) ;
1455
1456 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
1457
1458 /* RELAX3 parameter uses a double relop, ignore NaN case: */
1459 do_update = fnpiv > 0 &&
1460 (((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
1461 DEBUG2 (("do_update "ID"\n", do_update))
1462 }
1463 else
1464 {
1465 /* the current front would not be extended */
1466 do_update = fnpiv > 0 ;
1467 fnzeros = 0 ;
1468 DEBUG2 (("OUT-OUT do_update forced true: "ID"\n", do_update)) ;
1469
1470 /* The new front would be just big enough to hold the new
1471 * pivot row and column. */
1472 fnrows_new [OUT][OUT] = cdeg_out - 1 ;
1473 fncols_new [OUT][OUT] = rdeg [OUT][OUT] - 1 ;
1474 }
1475
1476 DEBUG2 (("option OUT OUT: nr "ID" nc "ID" cost "ID"\n",
1477 rdeg [OUT][OUT], cdeg_out, thiscost)) ;
1478
1479 if (bestcost == EMPTY || thiscost < bestcost)
1480 {
1481 /* this is the best option seen so far */
1482 Work->pivot_case = OUT_OUT ;
1483 bestcost = thiscost ;
1484 Work->do_extend = do_extend ;
1485 Work->do_update = do_update ;
1486 new_fnzeros = fnzeros ;
1487 }
1488 }
1489 }
1490
1491 /* At this point, a structural pivot has been found. */
1492 /* It may be numerically zero, however. */
1493 ASSERT (Work->pivot_case != EMPTY) ;
1494 DEBUG2 (("local search, best option "ID", best cost "ID"\n",
1495 Work->pivot_case, bestcost)) ;
1496
1497 /* ====================================================================== */
1498 /* Pivot row and column, and extension, now determined */
1499 /* ====================================================================== */
1500
1501 Work->fnzeros = new_fnzeros ;
1502
1503 /* ---------------------------------------------------------------------- */
1504 /* finalize the pivot row and column */
1505 /* ---------------------------------------------------------------------- */
1506
1507 switch (Work->pivot_case)
1508 {
1509 case IN_IN:
1510 DEBUG2 (("IN-IN option selected\n")) ;
1511 ASSERT (fnrows > 0 && fncols > 0) ;
1512 Work->pivcol_in_front = TRUE ;
1513 Work->pivrow_in_front = TRUE ;
1514 Work->pivcol = pivcol [IN] ;
1515 Work->pivrow = pivrow [IN][IN] ;
1516 Work->ccdeg = nr_in ;
1517 Work->Wrow = Fcols ;
1518 Work->rrdeg = rdeg [IN][IN] ;
1519 jj = jcand [IN] ;
1520 Work->fnrows_new = fnrows_new [IN][IN] ;
1521 Work->fncols_new = fncols_new [IN][IN] ;
1522 break ;
1523
1524 case IN_OUT:
1525 DEBUG2 (("IN-OUT option selected\n")) ;
1526 ASSERT (fnrows >= 0 && fncols > 0) ;
1527 Work->pivcol_in_front = TRUE ;
1528 Work->pivrow_in_front = FALSE ;
1529 Work->pivcol = pivcol [IN] ;
1530 Work->pivrow = pivrow [IN][OUT] ;
1531 Work->ccdeg = nr_in ;
1532 Work->Wrow = Wio ;
1533 Work->rrdeg = rdeg [IN][OUT] ;
1534 jj = jcand [IN] ;
1535 Work->fnrows_new = fnrows_new [IN][OUT] ;
1536 Work->fncols_new = fncols_new [IN][OUT] ;
1537 break ;
1538
1539 case OUT_IN:
1540 DEBUG2 (("OUT-IN option selected\n")) ;
1541 ASSERT (fnrows > 0 && fncols >= 0) ;
1542 Work->pivcol_in_front = FALSE ;
1543 Work->pivrow_in_front = TRUE ;
1544 Work->pivcol = pivcol [OUT] ;
1545 Work->pivrow = pivrow [OUT][IN] ;
1546 Work->ccdeg = cdeg_out ;
1547 /* Wrow might be equivalenced to Fcols (Freebie in): */
1548 Work->Wrow = Woi ;
1549 Work->rrdeg = rdeg [OUT][IN] ;
1550 /* Work->Wrow[0..fncols-1] is not there. See Fcols instead */
1551 jj = jcand [OUT] ;
1552 Work->fnrows_new = fnrows_new [OUT][IN] ;
1553 Work->fncols_new = fncols_new [OUT][IN] ;
1554 break ;
1555
1556 case OUT_OUT:
1557 DEBUG2 (("OUT-OUT option selected\n")) ;
1558 ASSERT (fnrows >= 0 && fncols >= 0) ;
1559 Work->pivcol_in_front = FALSE ;
1560 Work->pivrow_in_front = FALSE ;
1561 Work->pivcol = pivcol [OUT] ;
1562 Work->pivrow = pivrow [OUT][OUT] ;
1563 Work->ccdeg = cdeg_out ;
1564 /* Wrow might be equivalenced to Wio (Freebie out): */
1565 Work->Wrow = Woo ;
1566 Work->rrdeg = rdeg [OUT][OUT] ;
1567 jj = jcand [OUT] ;
1568 Work->fnrows_new = fnrows_new [OUT][OUT] ;
1569 Work->fncols_new = fncols_new [OUT][OUT] ;
1570 break ;
1571
1572 }
1573
1574 ASSERT (IMPLIES (fnrows == 0 && fncols == 0, Work->pivot_case == OUT_OUT)) ;
1575
1576 if (!Work->pivcol_in_front && pivcol [IN] != EMPTY)
1577 {
1578 /* clear Frpos if pivcol [IN] was searched, but not selected */
1579 for (i = fnrows ; i < cdeg_in ; i++)
1580 {
1581 Frpos [Frows [i]] = EMPTY;
1582 }
1583 }
1584
1585 /* ---------------------------------------------------------------------- */
1586 /* Pivot row and column have been found */
1587 /* ---------------------------------------------------------------------- */
1588
1589 /* ---------------------------------------------------------------------- */
1590 /* remove pivot column from candidate pivot column set */
1591 /* ---------------------------------------------------------------------- */
1592
1593 ASSERT (jj >= 0 && jj < Work->nCandidates) ;
1594 ASSERT (Work->pivcol == Work->Candidates [jj]) ;
1595 remove_candidate (jj, Work, Symbolic) ;
1596
1597 /* ---------------------------------------------------------------------- */
1598 /* check for frontal matrix growth */
1599 /* ---------------------------------------------------------------------- */
1600
1601 DEBUG1 (("Check frontal growth:\n")) ;
1602 DEBUG1 (("fnrows_new "ID" + 1 = "ID", fnr_curr "ID"\n",
1603 Work->fnrows_new, Work->fnrows_new + 1, fnr_curr)) ;
1604 DEBUG1 (("fncols_new "ID" + 1 = "ID", fnc_curr "ID"\n",
1605 Work->fncols_new, Work->fncols_new + 1, fnc_curr)) ;
1606
1607 Work->do_grow = (Work->fnrows_new + 1 > fnr_curr
1608 || Work->fncols_new + 1 > fnc_curr) ;
1609 if (Work->do_grow)
1610 {
1611 DEBUG0 (("\nNeed to grow frontal matrix, force do_update true\n")) ;
1612 /* If the front must grow, then apply the pending updates and remove
1613 * the current pivot rows/columns from the front prior to growing the
1614 * front. This frees up as much space as possible for the new front. */
1615 if (!Work->do_update && fnpiv > 0)
1616 {
1617 /* This update would not have to be done if the current front
1618 * was big enough. */
1619 Work->nforced++ ;
1620 Work->do_update = TRUE ;
1621 }
1622 }
1623
1624 /* ---------------------------------------------------------------------- */
1625 /* current pivot column */
1626 /* ---------------------------------------------------------------------- */
1627
1628 /*
1629 c1) If pivot column index is in the current front:
1630
1631 The pivot column pattern is in Frows [0 .. fnrows-1] and
1632 the extension is in Frows [fnrows ... fnrows+ccdeg-1].
1633
1634 Frpos [Frows [0 .. fnrows+ccdeg-1]] is
1635 equal to 0 .. fnrows+ccdeg-1. Wm is not needed.
1636
1637 The values are in Wy [0 .. fnrows+ccdeg-1].
1638
1639 c2) Otherwise, if the pivot column index is not in the current front:
1640
1641 c2a) If the front is being extended, old row indices in the the
1642 pivot column pattern are in Frows [0 .. fnrows-1].
1643
1644 All entries are in Wm [0 ... ccdeg-1], with values in
1645 Wx [0 .. ccdeg-1]. These may include entries already in
1646 Frows [0 .. fnrows-1].
1647
1648 Frpos [Frows [0 .. fnrows-1]] is equal to 0 .. fnrows-1.
1649 Frpos [Wm [0 .. ccdeg-1]] for new entries is < 0.
1650
1651 c2b) If the front is not being extended, then the entire pivot
1652 column pattern is in Wm [0 .. ccdeg-1]. It includes
1653 the pivot row index. It is does not contain the pattern
1654 Frows [0..fnrows-1]. The intersection of these two
1655 sets may or may not be empty. The values are in Wx [0..ccdeg-1]
1656
1657 In both cases c1 and c2, Frpos [Frows [0 .. fnrows-1]] is equal
1658 to 0 .. fnrows-1, which is the pattern of the current front.
1659 Any entry of Frpos that is not specified above is < 0.
1660 */
1661
1662
1663 #ifndef NDEBUG
1664 DEBUG2 (("\n\nSEARCH DONE: Pivot col "ID" in: ("ID") pivot row "ID" in: ("ID
1665 ") extend: "ID"\n\n", Work->pivcol, Work->pivcol_in_front,
1666 Work->pivrow, Work->pivrow_in_front, Work->do_extend)) ;
1667 UMF_dump_rowcol (1, Numeric, Work, Work->pivcol, !Symbolic->fixQ) ;
1668 DEBUG2 (("Pivot col "ID": fnrows "ID" ccdeg "ID"\n", Work->pivcol, fnrows,
1669 Work->ccdeg)) ;
1670 if (Work->pivcol_in_front) /* case c1 */
1671 {
1672 Int found = FALSE ;
1673 DEBUG3 (("Pivcol in front\n")) ;
1674 for (i = 0 ; i < fnrows ; i++)
1675 {
1676 row = Frows [i] ;
1677 DEBUG3 ((ID": row:: "ID" in front ", i, row)) ;
1678 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1679 ASSERT (Frpos [row] == i) ;
1680 EDEBUG3 (Wy [i]) ;
1681 if (row == Work->pivrow)
1682 {
1683 DEBUG3 ((" <- pivrow")) ;
1684 found = TRUE ;
1685 }
1686 DEBUG3 (("\n")) ;
1687 }
1688 ASSERT (found == Work->pivrow_in_front) ;
1689 found = FALSE ;
1690 for (i = fnrows ; i < fnrows + Work->ccdeg ; i++)
1691 {
1692 row = Frows [i] ;
1693 DEBUG3 ((ID": row:: "ID" (new)", i, row)) ;
1694 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1695 ASSERT (Frpos [row] == i) ;
1696 EDEBUG3 (Wy [i]) ;
1697 if (row == Work->pivrow)
1698 {
1699 DEBUG3 ((" <- pivrow")) ;
1700 found = TRUE ;
1701 }
1702 DEBUG3 (("\n")) ;
1703 }
1704 ASSERT (found == !Work->pivrow_in_front) ;
1705 }
1706 else
1707 {
1708 if (Work->do_extend)
1709 {
1710 Int found = FALSE ;
1711 DEBUG3 (("Pivcol not in front (extend)\n")) ;
1712 for (i = 0 ; i < fnrows ; i++)
1713 {
1714 row = Frows [i] ;
1715 DEBUG3 ((ID": row:: "ID" in front ", i, row)) ;
1716 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1717 ASSERT (Frpos [row] == i) ;
1718 if (row == Work->pivrow)
1719 {
1720 DEBUG3 ((" <- pivrow")) ;
1721 found = TRUE ;
1722 }
1723 DEBUG3 (("\n")) ;
1724 }
1725 ASSERT (found == Work->pivrow_in_front) ;
1726 found = FALSE ;
1727 DEBUG3 (("----\n")) ;
1728 for (i = 0 ; i < Work->ccdeg ; i++)
1729 {
1730 row = Wm [i] ;
1731 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1732 DEBUG3 ((ID": row:: "ID" ", i, row)) ;
1733 EDEBUG3 (Wx [i]) ;
1734 if (Frpos [row] < 0)
1735 {
1736 DEBUG3 ((" (new) ")) ;
1737 }
1738 if (row == Work->pivrow)
1739 {
1740 DEBUG3 ((" <- pivrow")) ;
1741 found = TRUE ;
1742 /* ... */
1743 if (Work->pivrow_in_front) ASSERT (Frpos [row] >= 0) ;
1744 else ASSERT (Frpos [row] < 0) ;
1745 }
1746 DEBUG3 (("\n")) ;
1747 }
1748 ASSERT (found) ;
1749 }
1750 else
1751 {
1752 Int found = FALSE ;
1753 DEBUG3 (("Pivcol not in front (no extend)\n")) ;
1754 for (i = 0 ; i < Work->ccdeg ; i++)
1755 {
1756 row = Wm [i] ;
1757 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1758 DEBUG3 ((ID": row:: "ID" ", i, row)) ;
1759 EDEBUG3 (Wx [i]) ;
1760 DEBUG3 ((" (new) ")) ;
1761 if (row == Work->pivrow)
1762 {
1763 DEBUG3 ((" <- pivrow")) ;
1764 found = TRUE ;
1765 }
1766 DEBUG3 (("\n")) ;
1767 }
1768 ASSERT (found) ;
1769 }
1770 }
1771 #endif
1772
1773 /* ---------------------------------------------------------------------- */
1774 /* current pivot row */
1775 /* ---------------------------------------------------------------------- */
1776
1777 /*
1778 r1) If the pivot row index is in the current front:
1779
1780 The pivot row pattern is in Fcols [0..fncols-1] and the extenson is
1781 in Wrow [fncols .. rrdeg-1]. If the pivot column is in the current
1782 front, then Fcols and Wrow are equivalenced.
1783
1784 r2) If the pivot row index is not in the current front:
1785
1786 r2a) If the front is being extended, the pivot row pattern is in
1787 Fcols [0 .. fncols-1]. New entries are in Wrow [0 .. rrdeg-1],
1788 but these may include entries already in Fcols [0 .. fncols-1].
1789
1790 r2b) Otherwise, the pivot row pattern is Wrow [0 .. rrdeg-1].
1791
1792 Fcpos [Fcols [0..fncols-1]] is (0..fncols-1) * fnr_curr.
1793 All other entries in Fcpos are < 0.
1794
1795 These conditions are asserted below.
1796
1797 ------------------------------------------------------------------------
1798 Other items in Work structure that are relevant:
1799
1800 pivcol: the pivot column index
1801 pivrow: the pivot column index
1802
1803 rrdeg:
1804 ccdeg:
1805
1806 fnrows: the number of rows in the currnt contribution block
1807 fncols: the number of columns in the current contribution block
1808
1809 fnrows_new: the number of rows in the new contribution block
1810 fncols_new: the number of rows in the new contribution block
1811
1812 ------------------------------------------------------------------------
1813 */
1814
1815
1816 #ifndef NDEBUG
1817 UMF_dump_rowcol (0, Numeric, Work, Work->pivrow, TRUE) ;
1818 DEBUG2 (("Pivot row "ID":\n", Work->pivrow)) ;
1819 if (Work->pivrow_in_front)
1820 {
1821 Int found = FALSE ;
1822 for (i = 0 ; i < fncols ; i++)
1823 {
1824 col = Fcols [i] ;
1825 DEBUG3 ((" col:: "ID" in front\n", col)) ;
1826 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1827 ASSERT (Fcpos [col] == i * fnr_curr) ;
1828 if (col == Work->pivcol) found = TRUE ;
1829 }
1830 ASSERT (found == Work->pivcol_in_front) ;
1831 found = FALSE ;
1832 ASSERT (IMPLIES (Work->pivcol_in_front, Fcols == Work->Wrow)) ;
1833 for (i = fncols ; i < Work->rrdeg ; i++)
1834 {
1835 col = Work->Wrow [i] ;
1836 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1837 ASSERT (Fcpos [col] < 0) ;
1838 if (col == Work->pivcol) found = TRUE ;
1839 else DEBUG3 ((" col:: "ID" (new)\n", col)) ;
1840 }
1841 ASSERT (found == !Work->pivcol_in_front) ;
1842 }
1843 else
1844 {
1845 if (Work->do_extend)
1846 {
1847 Int found = FALSE ;
1848 for (i = 0 ; i < fncols ; i++)
1849 {
1850 col = Fcols [i] ;
1851 DEBUG3 ((" col:: "ID" in front\n", col)) ;
1852 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1853 ASSERT (Fcpos [col] == i * fnr_curr) ;
1854 if (col == Work->pivcol) found = TRUE ;
1855 }
1856 ASSERT (found == Work->pivcol_in_front) ;
1857 found = FALSE ;
1858 for (i = 0 ; i < Work->rrdeg ; i++)
1859 {
1860 col = Work->Wrow [i] ;
1861 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1862 if (Fcpos [col] >= 0) continue ;
1863 if (col == Work->pivcol) found = TRUE ;
1864 else DEBUG3 ((" col:: "ID" (new, extend)\n", col)) ;
1865 }
1866 ASSERT (found == !Work->pivcol_in_front) ;
1867 }
1868 else
1869 {
1870 Int found = FALSE ;
1871 for (i = 0 ; i < Work->rrdeg ; i++)
1872 {
1873 col = Work->Wrow [i] ;
1874 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1875 if (col == Work->pivcol) found = TRUE ;
1876 else DEBUG3 ((" col:: "ID" (all new)\n", col)) ;
1877 }
1878 ASSERT (found) ;
1879 }
1880 }
1881 #endif
1882
1883 /* ---------------------------------------------------------------------- */
1884 /* determine whether to do scan2-row and scan2-col */
1885 /* ---------------------------------------------------------------------- */
1886
1887 if (Work->do_extend)
1888 {
1889 Work->do_scan2row = (fncols > 0) ;
1890 Work->do_scan2col = (fnrows > 0) ;
1891 }
1892 else
1893 {
1894 Work->do_scan2row = (fncols > 0) && Work->pivrow_in_front ;
1895 Work->do_scan2col = (fnrows > 0) && Work->pivcol_in_front ;
1896 }
1897
1898 /* ---------------------------------------------------------------------- */
1899
1900 DEBUG2 (("LOCAL SEARCH DONE: pivot column "ID" pivot row: "ID"\n",
1901 Work->pivcol, Work->pivrow)) ;
1902 DEBUG2 (("do_extend: "ID"\n", Work->do_extend)) ;
1903 DEBUG2 (("do_update: "ID"\n", Work->do_update)) ;
1904 DEBUG2 (("do_grow: "ID"\n", Work->do_grow)) ;
1905
1906 /* ---------------------------------------------------------------------- */
1907 /* keep track of the diagonal */
1908 /* ---------------------------------------------------------------------- */
1909
1910 if (Symbolic->prefer_diagonal
1911 && Work->pivcol < Work->n_col - Symbolic->nempty_col)
1912 {
1913 Diagonal_map = Work->Diagonal_map ;
1914 Diagonal_imap = Work->Diagonal_imap ;
1915 ASSERT (Diagonal_map != (Int *) NULL) ;
1916 ASSERT (Diagonal_imap != (Int *) NULL) ;
1917
1918 row2 = Diagonal_map [Work->pivcol] ;
1919 col2 = Diagonal_imap [Work->pivrow] ;
1920
1921 if (row2 < 0)
1922 {
1923 /* this was an off-diagonal pivot row */
1924 Work->noff_diagonal++ ;
1925 row2 = UNFLIP (row2) ;
1926 }
1927
1928 ASSERT (Diagonal_imap [row2] == Work->pivcol) ;
1929 ASSERT (UNFLIP (Diagonal_map [col2]) == Work->pivrow) ;
1930
1931 if (row2 != Work->pivrow)
1932 {
1933 /* swap the diagonal map to attempt to maintain symmetry later on.
1934 * Also mark the map for col2 (via FLIP) to denote that the entry
1935 * now on the diagonal is not the original entry on the diagonal. */
1936
1937 DEBUG0 (("Unsymmetric pivot\n")) ;
1938 Diagonal_map [Work->pivcol] = FLIP (Work->pivrow) ;
1939 Diagonal_imap [Work->pivrow] = Work->pivcol ;
1940
1941 Diagonal_map [col2] = FLIP (row2) ;
1942 Diagonal_imap [row2] = col2 ;
1943
1944 }
1945 ASSERT (n_row == n_col) ;
1946 #ifndef NDEBUG
1947 UMF_dump_diagonal_map (Diagonal_map, Diagonal_imap, Symbolic->n1,
1948 Symbolic->n_col, Symbolic->nempty_col) ;
1949 #endif
1950 }
1951
1952 return (UMFPACK_OK) ;
1953 }
1954