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