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