1
/* ========================================================================== */
2
/* === UMF_local_search ===================================================== */
3
/* ========================================================================== */
5
/* -------------------------------------------------------------------------- */
6
/* UMFPACK Version 4.1 (Apr. 30, 2003), Copyright (c) 2003 by Timothy A. */
7
/* Davis. All Rights Reserved. See ../README for License. */
8
/* email: davis@cise.ufl.edu CISE Department, Univ. of Florida. */
9
/* web: http://www.cise.ufl.edu/research/sparse/umfpack */
10
/* -------------------------------------------------------------------------- */
13
Perform pivot search to find pivot row and pivot column.
14
The pivot column is selected from the candidate set. The candidate set
15
corresponds to a supercolumn from colamd or UMF_analyze. The pivot column
16
is then removed from that set. Constructs the pivot column pattern and
17
values. Called by umf_kernel. Returns UMFPACK_OK if successful, or
18
UMFPACK_WARNING_singular_matrix or UMFPACK_ERROR_different_pattern if not.
21
#include "umf_internal.h"
22
#include "umf_row_search.h"
23
#include "umf_mem_free_tail_block.h"
25
/* Version 4.1: relaxed amalgamation control parameters are now fixed, and
26
* cannot be changed via Control [..] settings, as they could in Version 4.0. */
27
#define RELAX1 0.25 /* this was UMFPACK_DEFAULT_RELAXED_AMALGAMATION */
28
#define SYM_RELAX1 0.0 /* this is new to Version 4.1 */
29
#define RELAX2 0.1 /* this was UMFPACK_DEFAULT_RELAXED2_AMALGAMATION */
30
#define RELAX3 0.125 /* this was UMFPACK_DEFAULT_RELAXED3_AMALGAMATION */
32
/* ========================================================================== */
33
/* === remove_candidate ===================================================== */
34
/* ========================================================================== */
36
/* Remove a column from the set of candidate pivot columns. */
38
PRIVATE void remove_candidate (Int jj, WorkType *Work, SymbolicType *Symbolic)
43
DEBUGm2 (("pivot column Candidates before remove: nCand "ID" ncand "ID
44
" lo "ID" hi "ID" jj "ID"\n", Work->nCandidates, Work->ncand,
45
Work->lo, Work->hi, jj)) ;
46
for (j = 0 ; j < Work->nCandidates ; j++)
48
Int col = Work->Candidates [j] ;
49
DEBUGm2 ((ID" ", col));
50
ASSERT (col >= 0 && col < Work->n_col) ;
51
/* ASSERT (NON_PIVOTAL_COL (col)) ; */
52
ASSERT (col >= Work->lo && col <= Work->hi) ;
59
DEBUGm2 (("FixQ\n")) ;
60
/* do not modify the column ordering */
61
ASSERT (Work->nCandidates == 1) ;
65
Work->Candidates [0] = Work->nextcand++ ;
69
Work->nCandidates = 0 ;
74
/* place the next candidate in the set */
75
if (Work->ncand > MAX_CANDIDATES)
77
Work->Candidates [jj] = Work->nextcand++ ;
81
ASSERT (Work->nCandidates == Work->ncand) ;
82
Work->Candidates [jj] = Work->Candidates [Work->ncand - 1] ;
83
Work->Candidates [Work->ncand - 1] = EMPTY ;
90
DEBUGm2 (("pivot column Candidates after remove: nCand "ID" ncand "ID
91
" lo "ID" hi "ID" jj "ID"\n", Work->nCandidates, Work->ncand, Work->lo,
93
for (j = 0 ; j < Work->nCandidates ; j++)
95
Int col = Work->Candidates [j] ;
96
DEBUGm2 ((ID" ", col));
97
ASSERT (col >= 0 && col < Work->n_col) ;
98
/* ASSERT (NON_PIVOTAL_COL (col)) ; */
99
ASSERT (col >= Work->lo && col <= Work->hi) ;
102
ASSERT (Work->ncand >= 0) ;
103
ASSERT (Work->nCandidates <= Work->ncand) ;
107
/* ========================================================================== */
108
/* === UMF_local_search ===================================================== */
109
/* ========================================================================== */
111
GLOBAL Int UMF_local_search
113
NumericType *Numeric,
115
SymbolicType *Symbolic
118
/* ---------------------------------------------------------------------- */
119
/* local variables */
120
/* ---------------------------------------------------------------------- */
122
Entry *Flblock, *Fublock, *Fs, *Fcblock, *C, *Wx, *Wy, *Fu, *Flublock,
125
Int pos, nrows, *Cols, *Rows, e, f, status, max_cdeg, fnzeros, nb, j, col,
126
i, row, cdeg_in, rdeg [2][2], fnpiv, nothing [2], new_LUsize,
127
pivrow [2][2], pivcol [2], *Wp, *Fcpos, *Frpos, new_fnzeros, cdeg_out,
128
*Wm, *Wio, *Woi, *Woo, *Frows, *Fcols, fnrows, fncols, *E, deg, nr_in,
129
nc, thiscost, bestcost, nr_out, do_update, extra_cols, extra_rows,
130
extra_zeros, relaxed_front, do_extend, fnr_curr, fnc_curr, tpi,
131
*Col_tuples, *Col_degree, *Col_tlen, jj, jcand [2], freebie [2],
132
did_rowmerge, fnrows_new [2][2], fncols_new [2][2], search_pivcol_out,
133
*Diagonal_map, *Diagonal_imap, row2, col2 ;
135
Tuple *tp, *tpend, *tp1, *tp2 ;
139
Int debug_ok, n_row, n_col, *Row_degree ;
140
Row_degree = Numeric->Rperm ; /* for NON_PIVOTAL_ROW macro only */
143
/* ---------------------------------------------------------------------- */
145
/* ---------------------------------------------------------------------- */
147
Memory = Numeric->Memory ;
149
Col_degree = Numeric->Cperm ;
151
Col_tuples = Numeric->Lip ;
152
Col_tlen = Numeric->Lilen ;
161
Fcpos = Work->Fcpos ;
162
Frpos = Work->Frpos ;
163
Frows = Work->Frows ;
164
Fcols = Work->Fcols ;
165
fnrows = Work->fnrows ;
166
fncols = Work->fncols ;
168
fnr_curr = Work->fnr_curr ;
169
fnc_curr = Work->fnc_curr ;
170
fnpiv = Work->fnpiv ;
171
nothing [0] = EMPTY ;
172
nothing [1] = EMPTY ;
173
relax1 = (Symbolic->prefer_diagonal) ? SYM_RELAX1 : RELAX1 ;
174
fnzeros = Work->fnzeros ;
175
new_fnzeros = fnzeros ;
178
Fcblock = Work->Fcblock ; /* current contribution block */
179
Flblock = Work->Flblock ; /* current L block */
180
Fublock = Work->Fublock ; /* current U block */
181
Flublock = Work->Flublock ; /* current LU block */
183
/* The pivot column degree cannot exceed max_cdeg */
184
max_cdeg = Work->fnrows_max ;
185
ASSERT (Work->fnrows_max <= Symbolic->maxnrows) ;
186
ASSERT (Work->fncols_max <= Symbolic->maxncols) ;
188
if (fnrows == 0 && fncols == 0)
190
/* frontal matrix is empty */
191
Work->firstsuper = Work->ksuper ;
195
n_row = Work->n_row ;
196
n_col = Work->n_col ;
197
DEBUG2 (("\n========LOCAL SEARCH: current frontal matrix: ========= \n")) ;
198
UMF_dump_current_front (Numeric, Work, TRUE) ;
199
if (UMF_debug > 0 || MAX (n_row, n_col) < 1000)
201
for (i = 0 ; i < MAX (n_row, n_col) ; i++)
203
ASSERT (Wp [i] < 0) ;
207
DEBUGm2 ((ID" pivot column Candidates: lo "ID" hi "ID"\n",
208
Work->nCandidates, Work->lo, Work->hi)) ;
209
for (j = 0 ; j < Work->nCandidates ; j++)
211
col = Work->Candidates [j] ;
212
DEBUGm2 ((ID" ", col));
213
ASSERT (col >= 0 && col < n_col) ;
214
ASSERT (NON_PIVOTAL_COL (col)) ;
215
ASSERT (col >= Work->lo && col <= Work->hi) ;
219
/* there are no 0-by-c or r-by-0 fronts, where c and r are > 0 */
220
/* a front is either 0-by-0, or r-by-c */
221
DEBUG2 (("\n\n::: "ID" : Npiv: "ID" + fnpiv "ID" = "ID". "
222
"size "ID"-by-"ID"\n", Work->frontid,
223
Work->npiv, Work->fnpiv, Work->npiv + Work->fnpiv, fnrows, fncols)) ;
224
ASSERT ((fnrows == 0 && fncols == 0) ||(fnrows != 0 && fncols != 0)) ;
227
/* ====================================================================== */
228
/* === PIVOT SEARCH ===================================================== */
229
/* ====================================================================== */
233
pivcol [IN] = EMPTY ;
234
pivcol [OUT] = EMPTY ;
239
pivrow [IN][IN] = EMPTY ;
240
pivrow [IN][OUT] = EMPTY ;
241
pivrow [OUT][IN] = EMPTY ;
242
pivrow [OUT][OUT] = EMPTY ;
244
rdeg [IN][IN] = Int_MAX ;
245
rdeg [IN][OUT] = Int_MAX ;
246
rdeg [OUT][IN] = Int_MAX ;
247
rdeg [OUT][OUT] = Int_MAX ;
249
freebie [IN] = FALSE ;
250
freebie [OUT] = FALSE ;
252
Work->pivot_case = EMPTY ;
259
jcand [OUT] = EMPTY ;
261
fnrows_new [IN][IN] = EMPTY ;
262
fnrows_new [IN][OUT] = EMPTY ;
263
fnrows_new [OUT][IN] = EMPTY ;
264
fnrows_new [OUT][OUT] = EMPTY ;
266
fncols_new [IN][IN] = EMPTY ;
267
fncols_new [IN][OUT] = EMPTY ;
268
fncols_new [OUT][IN] = EMPTY ;
269
fncols_new [OUT][OUT] = EMPTY ;
273
DEBUG4 (("Check Frpos : fnrows "ID" col "ID" maxcdeg "ID"\n",
274
fnrows, pivcol [IN], max_cdeg)) ;
275
for (i = 0 ; i < fnrows ; i++)
278
DEBUG4 ((" row: "ID"\n", row)) ;
279
ASSERT (row >= 0 && row < n_row) ;
280
ASSERT (Frpos [row] == i) ;
282
DEBUG4 (("All:\n")) ;
283
if (UMF_debug > 0 || n_row < 1000)
286
for (row = 0 ; row < n_row ; row++)
288
if (Frpos [row] == EMPTY)
294
DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
297
ASSERT (cnt == n_row) ;
301
/* ---------------------------------------------------------------------- */
302
/* find shortest column in the front, and shortest column not in the */
303
/* front, from the candidate pivot column set */
304
/* ---------------------------------------------------------------------- */
306
/* If there are too many candidates, then only look at the first */
307
/* MAX_CANDIDATES of them. Otherwise, if there are O(n) candidates, */
308
/* this code could take O(n^2) time. */
310
/* ------------------------------------------------------------------ */
311
/* look in the candidate set for the best column */
312
/* ------------------------------------------------------------------ */
314
DEBUG2 (("Max candidates %d, Work->ncand "ID" jmax "ID"\n",
315
MAX_CANDIDATES, Work->ncand, Work->nCandidates)) ;
316
col = Work->Candidates [0] ;
317
ASSERT (Work->nCandidates > 0) ;
318
DEBUG3 (("Pivot column candidate: "ID" j = "ID"\n", col, j)) ;
319
ASSERT (col >= 0 && col < n_col) ;
321
/* there is no Col_degree if fixQ is true */
322
deg = Symbolic->fixQ ? EMPTY : Col_degree [col] ;
325
DEBUG3 (("Pivot column candidate: "ID" cost: "ID" Fcpos[col] "ID"\n",
326
col, deg, Fcpos [col])) ;
327
UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
330
DEBUG1 (("FIXQ: Candidates "ID" pivcol "ID" npiv "ID" fnpiv "ID
331
" ndiscard "ID "\n", Work->nCandidates, col, Work->npiv,
332
Work->fnpiv, Work->ndiscard)) ;
333
ASSERT (Work->nCandidates == 1) ;
334
ASSERT (col == Work->npiv + Work->fnpiv + Work->ndiscard) ;
338
if (Fcpos [col] >= 0)
340
/* best column in front, so far */
342
cdeg_in = deg ; /* ignored, if fixQ is true */
347
/* best column not in front, so far */
349
cdeg_out = deg ; /* ignored, if fixQ is true */
353
/* look at the rest of the candidates */
354
for (j = 1 ; j < Work->nCandidates ; j++)
356
col = Work->Candidates [j] ;
358
DEBUG3 (("Pivot col candidate: "ID" j = "ID"\n", col, j)) ;
359
ASSERT (col >= 0 && col < n_col) ;
360
ASSERT (!Symbolic->fixQ) ;
361
deg = Col_degree [col] ;
363
DEBUG3 (("Pivot col candidate: "ID" cost: "ID" Fcpos[col] "ID"\n",
364
col, deg, Fcpos [col])) ;
365
UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
367
if (Fcpos [col] >= 0)
371
fs = Fcpos [col] / fnr_curr ;
372
ASSERT (fs >= 0 && fs < fncols) ;
374
if (deg < cdeg_in || (deg == cdeg_in && col < pivcol [IN]))
376
/* best column in front, so far */
384
if (deg < cdeg_out || (deg == cdeg_out && col < pivcol [OUT]))
386
/* best column not in front, so far */
394
DEBUG2 (("Pivcol in "ID" out "ID"\n", pivcol [IN], pivcol [OUT])) ;
395
ASSERT ((pivcol [IN] >= 0 && pivcol [IN] < n_col)
396
|| (pivcol [OUT] >= 0 && pivcol [OUT] < n_col)) ;
401
/* ---------------------------------------------------------------------- */
402
/* construct candidate column in front, and search for pivot rows */
403
/* ---------------------------------------------------------------------- */
407
DEBUG4 (("Prior to col update: fnrows "ID" col "ID" maxcdeg "ID"\n",
408
fnrows, pivcol [IN], max_cdeg)) ;
409
for (i = 0 ; i < fnrows ; i++)
412
DEBUG4 ((" row: "ID"\n", row)) ;
413
ASSERT (row >= 0 && row < n_row) ;
414
ASSERT (Frpos [row] == i) ;
416
DEBUG4 (("All:\n")) ;
417
if (UMF_debug > 0 || n_row < 1000)
420
for (row = 0 ; row < n_row ; row++)
422
if (Frpos [row] == EMPTY)
428
DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
431
ASSERT (cnt == n_row) ;
435
if (pivcol [IN] != EMPTY)
439
DEBUG2 (("col[IN] column "ID" in front at position = "ID"\n",
440
pivcol [IN], Fcpos [pivcol [IN]])) ;
441
UMF_dump_rowcol (1, Numeric, Work, pivcol [IN], !Symbolic->fixQ) ;
444
/* the only way we can have a pivcol[IN] is if the front is not empty */
445
ASSERT (fnrows > 0 && fncols > 0) ;
447
DEBUG4 (("Update pivot column:\n")) ;
448
Fs = Fcblock + Fcpos [pivcol [IN]] ;
449
Fu = Fublock + (Fcpos [pivcol [IN]] / fnr_curr) ;
450
Flu = Flublock + fnpiv * nb ;
452
/* ------------------------------------------------------------------ */
453
/* copy the pivot column from the U block into the LU block */
454
/* ------------------------------------------------------------------ */
456
/* This copy is permanent if the pivcol [IN] is chosen. */
457
for (i = 0 ; i < fnpiv ; i++)
459
Flu [i] = Fu [i*fnc_curr] ;
462
/* ------------------------------------------------------------------ */
463
/* update the pivot column in the LU block using a triangular solve */
464
/* ------------------------------------------------------------------ */
466
/* This work will be discarded if the pivcol [OUT] is chosen instead.
467
* It is permanent if the pivcol [IN] is chosen. */
471
/* solve Lx=b, where b = U (:,k), stored in the LU block */
475
/* no BLAS available - use plain C code instead */
476
Entry *Flub = Flublock ;
477
for (j = 0 ; j < fnpiv ; j++)
479
Entry Fuj = Flu [j] ;
481
for (i = j+1 ; i < fnpiv ; i++)
483
/* Flu [i] -= Flublock [i + j*nb] * Flu [j] ; */
484
MULT_SUB (Flu [i], Flub [i], Fuj) ;
490
BLAS_TRSV (fnpiv, Flublock, Flu, nb) ;
495
/* ------------------------------------------------------------------ */
496
/* copy the pivot column from the C block into Wy */
497
/* ------------------------------------------------------------------ */
499
for (i = 0 ; i < fnrows ; i++)
504
/* ------------------------------------------------------------------ */
505
/* update the pivot column of L using a matrix-vector multiply */
506
/* ------------------------------------------------------------------ */
508
/* this work will be discarded if the pivcol [OUT] is chosen instead */
511
/* no BLAS available - use plain C code instead */
512
for (j = 0 ; j < fnpiv ; j++)
514
Entry Fuj, *Flub = Flblock + j * fnr_curr ;
516
if (IS_NONZERO (Fuj))
519
for (i = 0 ; i < fnrows ; i++)
521
/* Wy [i] -= Flblock [i+j*fnr_curr] * Fuj ; */
522
MULT_SUB (Wy [i], Flub [i], Fuj) ;
525
/* Flblock += fnr_curr ; */
528
/* Using 1-based notation:
529
* Wy (1:fnrows) -= Flblock (1:fnrows,1:fnpiv) * Flu (1:fnpiv) */
530
BLAS_GEMV (fnrows, fnpiv, Flblock, Flu, Wy, fnr_curr) ;
533
/* ------------------------------------------------------------------ */
536
DEBUG2 (("Wy after update: fnrows="ID"\n", fnrows)) ;
537
DEBUG4 ((" fnpiv="ID" \n", fnpiv)) ;
538
for (i = 0 ; i < fnrows ; i++)
540
DEBUG4 ((ID" "ID" "ID, i, Frows [i], Frpos [Frows [i]])) ;
546
/* ------------------------------------------------------------------ */
547
/* construct the candidate column */
548
/* ------------------------------------------------------------------ */
554
DEBUG4 (("After col update: fnrows "ID" col "ID" maxcdeg "ID"\n",
555
fnrows, pivcol [IN], max_cdeg)) ;
556
for (i = 0 ; i < fnrows ; i++)
559
DEBUG4 ((" row: "ID"\n", row)) ;
560
ASSERT (row >= 0 && row < n_row) ;
561
ASSERT (Frpos [row] == i) ;
563
DEBUG4 (("All:\n")) ;
564
if (UMF_debug > 0 || n_row < 1000)
567
for (row = 0 ; row < n_row ; row++)
569
if (Frpos [row] == EMPTY)
575
DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
578
ASSERT (cnt == n_row) ;
584
DEBUG4 (("COL ASSEMBLE: cdeg "ID"\nREDUCE COL in "ID" max_cdeg "ID"\n",
585
cdeg_in, pivcol [IN], max_cdeg)) ;
586
for (i = 0 ; i < cdeg_in ; i++)
589
ASSERT (row >= 0 && row < n_row) ;
590
ASSERT (Frpos [row] == i) ;
592
if (UMF_debug > 0 || n_row < 1000)
595
for (row = 0 ; row < n_row ; row++)
597
if (Frpos [row] == EMPTY) cnt++ ;
599
ASSERT (cnt == n_row) ;
603
/* assemble column into Wy */
605
ASSERT (pivcol [IN] >= 0 && pivcol [IN] < n_col) ;
606
ASSERT (NON_PIVOTAL_COL (pivcol [IN])) ;
608
tpi = Col_tuples [pivcol [IN]] ;
611
tp = (Tuple *) (Memory + tpi) ;
614
tpend = tp + Col_tlen [pivcol [IN]] ;
615
for ( ; tp < tpend ; tp++)
618
ASSERT (e > 0 && e <= Work->nel) ;
619
if (!E [e]) continue ; /* element already deallocated */
623
p += UNITS (Element, 1) ;
625
if (Cols [f] == EMPTY) continue ; /* column already assembled */
626
ASSERT (pivcol [IN] == Cols [f]) ;
628
Rows = Cols + ep->ncols ;
630
p += UNITS (Int, ep->ncols + nrows) ;
631
C = ((Entry *) p) + f * nrows ;
633
for (i = 0 ; i < nrows ; i++)
636
if (row >= 0) /* skip this if already gone from element */
638
ASSERT (row < n_row) ;
642
/* new entry in the pattern - save Frpos */
643
ASSERT (cdeg_in < n_row) ;
644
if (cdeg_in >= max_cdeg)
646
/* :: pattern change (cdeg in failure) :: */
647
DEBUGm4 (("cdeg_in failure\n")) ;
648
return (UMFPACK_ERROR_different_pattern) ;
650
Frpos [row] = cdeg_in ;
651
Frows [cdeg_in] = row ;
652
Wy [cdeg_in++] = C [i] ;
656
/* entry already in pattern - sum values in Wy */
657
/* Wy [pos] += C [i] ; */
658
ASSERT (pos < max_cdeg) ;
659
ASSEMBLE (Wy [pos], C [i]) ;
663
*tp2++ = *tp ; /* leave the tuple in the list */
665
Col_tlen [pivcol [IN]] = tp2 - tp1 ;
668
/* ------------------------------------------------------------------ */
671
/* check Frpos again */
672
DEBUG4 (("COL DONE: cdeg "ID"\nREDUCE COL in "ID" max_cdeg "ID"\n",
673
cdeg_in, pivcol [IN], max_cdeg)) ;
674
for (i = 0 ; i < cdeg_in ; i++)
677
ASSERT (row >= 0 && row < n_row) ;
678
ASSERT (Frpos [row] == i) ;
680
if (UMF_debug > 0 || n_row < 1000)
683
for (row = 0 ; row < n_row ; row++)
685
if (Frpos [row] == EMPTY) cnt++ ;
687
ASSERT (cnt == n_row) ;
692
DEBUG4 (("Reduced column: cdeg in "ID" fnrows_max "ID"\n",
693
cdeg_in, Work->fnrows_max)) ;
694
for (i = 0 ; i < cdeg_in ; i++)
696
DEBUG4 ((" "ID" "ID" "ID, i, Frows [i], Frpos [Frows [i]])) ;
699
ASSERT (i == Frpos [Frows [i]]) ;
701
ASSERT (cdeg_in <= Work->fnrows_max) ;
704
/* ------------------------------------------------------------------ */
705
/* cdeg_in is now the exact degree of this column */
706
/* ------------------------------------------------------------------ */
708
nr_in = cdeg_in - fnrows ;
710
/* since there are no 0-by-x fronts, if there is a pivcol [IN] the */
711
/* front must have at least one row. */
712
ASSERT (cdeg_in > 0) ;
714
/* new degree of pivcol [IN], excluding current front is nr_in */
715
/* column expands by nr_in rows */
717
/* ------------------------------------------------------------------ */
718
/* search for two candidate pivot rows */
719
/* ------------------------------------------------------------------ */
721
/* for the IN_IN pivot row (if any), */
722
/* extend the pattern in place, using Fcols */
723
status = UMF_row_search (Numeric, Work, Symbolic,
724
fnrows, cdeg_in, Frows, Frpos, /* pattern of column to search */
725
pivrow [IN], rdeg [IN], Fcols, Wio, nothing, Wy,
726
pivcol [IN], freebie) ;
727
ASSERT (!freebie [IN] && !freebie [OUT]) ;
729
/* ------------------------------------------------------------------ */
730
/* fatal error if matrix pattern has changed since symbolic analysis */
731
/* ------------------------------------------------------------------ */
733
if (status == UMFPACK_ERROR_different_pattern)
735
/* :: pattern change (row search IN failure) :: */
736
DEBUGm4 (("row search IN failure\n")) ;
737
return (UMFPACK_ERROR_different_pattern) ;
740
/* ------------------------------------------------------------------ */
741
/* we now must have a structural pivot */
742
/* ------------------------------------------------------------------ */
744
/* Since the pivcol[IN] exists, there must be at least one row in the */
745
/* current frontal matrix, and so we must have found a structural */
746
/* pivot. The numerical value might be zero, of course. */
748
ASSERT (status != UMFPACK_WARNING_singular_matrix) ;
750
/* ------------------------------------------------------------------ */
751
/* evaluate IN_IN option */
752
/* ------------------------------------------------------------------ */
754
if (pivrow [IN][IN] != EMPTY)
756
/* The current front would become an (implicit) LUson.
757
* Both candidate pivot row and column are in the current front.
758
* Cost is how much the current front would expand */
760
/* pivrow[IN][IN] candidates are not found via row merge search */
762
ASSERT (fnrows >= 0 && fncols >= 0) ;
764
ASSERT (cdeg_in > 0) ;
765
nc = rdeg [IN][IN] - fncols ;
768
/* each column in front (except pivot column) grows by nr_in: */
769
(nr_in * (fncols - 1)) +
770
/* new columns not in old front: */
771
(nc * (cdeg_in - 1)) ;
773
/* no extra cost to relaxed amalgamation */
775
ASSERT (fnrows + nr_in == cdeg_in) ;
776
ASSERT (fncols + nc == rdeg [IN][IN]) ;
778
/* size of relaxed front (after pivot row column removed): */
779
fnrows_new [IN][IN] = (fnrows-1) + nr_in ;
780
fncols_new [IN][IN] = (fncols-1) + nc ;
781
/* relaxed_front = fnrows_new [IN][IN] * fncols_new [IN][IN] ; */
785
DEBUG2 (("Evaluating option IN-IN:\n")) ;
786
DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_in "ID" nc "ID"\n",
787
Work->fnzeros, fnpiv, nr_in, nc)) ;
788
DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
790
/* determine if BLAS-3 update should be applied before extending. */
791
/* update if too many zero entries accumulate in the LU block */
792
fnzeros = Work->fnzeros + fnpiv * (nr_in + nc) ;
794
DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
796
new_LUsize = (fnpiv+1) * (fnrows + nr_in + fncols + nc) ;
798
DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
800
/* There are fnpiv pivots currently in the front. This one
801
* will be the (fnpiv+1)st pivot, if it is extended. */
803
/* RELAX2 parameter uses a double relop, but ignore NaN case: */
804
do_update = fnpiv > 0 &&
805
(((double) fnzeros) / ((double) new_LUsize)) > RELAX2 ;
807
DEBUG2 (("do_update "ID"\n", do_update))
809
DEBUG2 (("option IN IN : nr "ID" nc "ID" cost "ID"(0) relax "ID
810
"\n", nr_in, nc, thiscost, do_extend)) ;
812
/* this is the best option seen so far */
813
Work->pivot_case = IN_IN ;
814
bestcost = thiscost ;
816
/* do the amalgamation and extend the front */
817
Work->do_extend = do_extend ;
818
Work->do_update = do_update ;
819
new_fnzeros = fnzeros ;
823
/* ------------------------------------------------------------------ */
824
/* evaluate IN_OUT option */
825
/* ------------------------------------------------------------------ */
827
if (pivrow [IN][OUT] != EMPTY)
829
/* The current front would become a Uson of the new front.
830
* Candidate pivot column is in the current front, but the
831
* candidate pivot row is not. */
833
ASSERT (fnrows >= 0 && fncols > 0) ;
834
ASSERT (cdeg_in > 0) ;
836
/* must be at least one row outside the front */
837
/* (the pivrow [IN][OUT] itself) */
838
ASSERT (nr_in >= 1) ;
840
/* count columns not in current front */
845
for (i = 0 ; i < rdeg [IN][OUT] ; i++)
848
DEBUG4 (("counting col "ID" Fcpos[] = "ID"\n", col,
850
ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
851
if (Fcpos [col] < 0) nc++ ;
853
/* we must see the pivot column somewhere */
854
if (col == pivcol [IN])
856
ASSERT (Fcpos [col] >= 0) ;
864
/* each row in front grows by nc: */
866
/* new rows not affected by front: */
867
((nr_in-1) * (rdeg [IN][OUT]-1)) ;
869
/* check the cost of relaxed IN_OUT amalgamation */
871
extra_cols = ((fncols-1) + nc ) - (rdeg [IN][OUT] - 1) ;
872
ASSERT (extra_cols >= 0) ;
873
ASSERT (fncols + nc == extra_cols + rdeg [IN][OUT]) ;
874
extra_zeros = (nr_in-1) * extra_cols ; /* symbolic fill-in */
876
ASSERT (fnrows + nr_in == cdeg_in) ;
877
ASSERT (fncols + nc == rdeg [IN][OUT] + extra_cols) ;
879
/* size of relaxed front (after pivot column removed): */
880
fnrows_new [IN][OUT] = fnrows + (nr_in-1) ;
881
fncols_new [IN][OUT] = (fncols-1) + nc ;
882
relaxed_front = fnrows_new [IN][OUT] * fncols_new [IN][OUT] ;
884
/* do relaxed amalgamation if the extra zeros are no more */
885
/* than a fraction (default 0.25) of the relaxed front */
886
/* if relax = 0: no extra zeros allowed */
887
/* if relax = +inf: always amalgamate */
889
/* relax parameter uses a double relop, but ignore NaN case: */
890
if (extra_zeros == 0)
896
do_extend = ((double) extra_zeros) <
897
(relax1 * (double) relaxed_front) ;
902
/* count the cost of relaxed amalgamation */
903
thiscost += extra_zeros ;
905
DEBUG2 (("Evaluating option IN-OUT:\n")) ;
906
DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_in "ID" nc "ID"\n",
907
Work->fnzeros, fnpiv, nr_in, nc)) ;
908
DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
910
/* determine if BLAS-3 update to be applied before extending. */
911
/* update if too many zero entries accumulate in the LU block */
912
fnzeros = Work->fnzeros + fnpiv * (nr_in + nc) ;
914
DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
916
new_LUsize = (fnpiv+1) * (fnrows + nr_in + fncols + nc) ;
918
DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
920
/* RELAX3 parameter uses a double relop, ignore NaN case: */
921
do_update = fnpiv > 0 &&
922
(((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
923
DEBUG2 (("do_update "ID"\n", do_update))
928
/* the current front would not be extended */
929
do_update = fnpiv > 0 ;
931
DEBUG2 (("IN-OUT do_update forced true: "ID"\n", do_update)) ;
933
/* The new front would be just big enough to hold the new
934
* pivot row and column. */
935
fnrows_new [IN][OUT] = cdeg_in - 1 ;
936
fncols_new [IN][OUT] = rdeg [IN][OUT] - 1 ;
940
DEBUG2 (("option IN OUT: nr "ID" nc "ID" cost "ID"("ID") relax "ID
941
"\n", nr_in, nc, thiscost, extra_zeros, do_extend)) ;
943
if (bestcost == EMPTY || thiscost < bestcost)
945
/* this is the best option seen so far */
946
Work->pivot_case = IN_OUT ;
947
bestcost = thiscost ;
948
Work->do_extend = do_extend ;
949
Work->do_update = do_update ;
950
new_fnzeros = fnzeros ;
955
/* ---------------------------------------------------------------------- */
956
/* construct candidate column not in front, and search for pivot rows */
957
/* ---------------------------------------------------------------------- */
959
search_pivcol_out = (bestcost != 0 && pivcol [OUT] != EMPTY) ;
960
if (Symbolic->prefer_diagonal)
962
search_pivcol_out = search_pivcol_out && (pivrow [IN][IN] == EMPTY) ;
965
if (search_pivcol_out)
969
DEBUG2 (("out_col column "ID" NOT in front at position = "ID"\n",
970
pivcol [OUT], Fcpos [pivcol [OUT]])) ;
971
UMF_dump_rowcol (1, Numeric, Work, pivcol [OUT], !Symbolic->fixQ) ;
972
DEBUG2 (("fncols "ID" fncols_max "ID"\n", fncols, Work->fncols_max)) ;
973
ASSERT (fncols < Work->fncols_max) ;
976
/* Use Wx as temporary workspace to construct the pivcol [OUT] */
979
/* ------------------------------------------------------------------ */
980
/* construct the candidate column (currently not in the front) */
981
/* ------------------------------------------------------------------ */
983
/* Construct the column in Wx, Wm, using Wp for the positions: */
984
/* Wm [0..cdeg_out-1] list of row indices in the column */
985
/* Wx [0..cdeg_out-1] list of corresponding numerical values */
986
/* Wp [0..n-1] starts as all negative, and ends that way too. */
992
DEBUG4 (("COL ASSEMBLE: cdeg 0\nREDUCE COL out "ID"\n", pivcol [OUT])) ;
993
if (UMF_debug > 0 || MAX (n_row, n_col) < 1000)
995
for (i = 0 ; i < MAX (n_row, n_col) ; i++)
997
ASSERT (Wp [i] < 0) ;
1000
DEBUG4 (("max_cdeg: "ID"\n", max_cdeg)) ;
1003
ASSERT (pivcol [OUT] >= 0 && pivcol [OUT] < n_col) ;
1004
ASSERT (NON_PIVOTAL_COL (pivcol [OUT])) ;
1006
tpi = Col_tuples [pivcol [OUT]] ;
1009
tp = (Tuple *) (Memory + tpi) ;
1012
tpend = tp + Col_tlen [pivcol [OUT]] ;
1013
for ( ; tp < tpend ; tp++)
1016
ASSERT (e > 0 && e <= Work->nel) ;
1017
if (!E [e]) continue ; /* element already deallocated */
1019
p = Memory + E [e] ;
1020
ep = (Element *) p ;
1021
p += UNITS (Element, 1) ;
1023
if (Cols [f] == EMPTY) continue ; /* column already assembled */
1024
ASSERT (pivcol [OUT] == Cols [f]) ;
1026
Rows = Cols + ep->ncols ;
1028
p += UNITS (Int, ep->ncols + nrows) ;
1029
C = ((Entry *) p) + f * nrows ;
1031
for (i = 0 ; i < nrows ; i++)
1034
if (row >= 0) /* skip this if already gone from element */
1036
ASSERT (row < n_row) ;
1040
/* new entry in the pattern - save Wp */
1041
ASSERT (cdeg_out < n_row) ;
1042
if (cdeg_out >= max_cdeg)
1044
/* :: pattern change (cdeg out failure) :: */
1045
DEBUGm4 (("cdeg out failure\n")) ;
1046
return (UMFPACK_ERROR_different_pattern) ;
1048
Wp [row] = cdeg_out ;
1049
Wm [cdeg_out] = row ;
1050
Wx [cdeg_out++] = C [i] ;
1054
/* entry already in pattern - sum the values */
1055
/* Wx [pos] += C [i] ; */
1056
ASSEMBLE (Wx [pos], C [i]) ;
1060
*tp2++ = *tp ; /* leave the tuple in the list */
1062
Col_tlen [pivcol [OUT]] = tp2 - tp1 ;
1065
/* ------------------------------------------------------------------ */
1068
DEBUG4 (("Reduced column: cdeg out "ID"\n", cdeg_out)) ;
1069
for (i = 0 ; i < cdeg_out ; i++)
1071
DEBUG4 ((" "ID" "ID" "ID, i, Wm [i], Wp [Wm [i]])) ;
1074
ASSERT (i == Wp [Wm [i]]) ;
1078
/* ------------------------------------------------------------------ */
1079
/* new degree of pivcol [OUT] is cdeg_out */
1080
/* ------------------------------------------------------------------ */
1082
/* search for two candidate pivot rows */
1083
status = UMF_row_search (Numeric, Work, Symbolic,
1084
0, cdeg_out, Wm, Wp, /* pattern of column to search */
1085
pivrow [OUT], rdeg [OUT], Woi, Woo, pivrow [IN], Wx,
1086
pivcol [OUT], freebie) ;
1088
/* ------------------------------------------------------------------ */
1089
/* fatal error if matrix pattern has changed since symbolic analysis */
1090
/* ------------------------------------------------------------------ */
1092
if (status == UMFPACK_ERROR_different_pattern)
1094
/* :: pattern change detected in umf_local_search :: */
1095
return (UMFPACK_ERROR_different_pattern) ;
1098
/* ------------------------------------------------------------------ */
1100
/* ------------------------------------------------------------------ */
1102
for (i = 0 ; i < cdeg_out ; i++)
1104
Wp [Wm [i]] = EMPTY ; /* clear Wp */
1107
/* ------------------------------------------------------------------ */
1108
/* check for rectangular, singular matrix */
1109
/* ------------------------------------------------------------------ */
1111
if (status == UMFPACK_WARNING_singular_matrix)
1113
/* Pivot column is empty, and row-merge set is empty too. The
1114
* matrix is structurally singular. The current frontal matrix must
1115
* be empty, too. It it weren't, and pivcol [OUT] exists, then
1116
* there would be at least one row that could be selected. Since
1117
* the current front is empty, pivcol [IN] must also be EMPTY.
1120
DEBUGm4 (("Note: pivcol [OUT]: "ID" discard\n", pivcol [OUT])) ;
1121
ASSERT ((Work->fnrows == 0 && Work->fncols == 0)) ;
1122
ASSERT (pivcol [IN] == EMPTY) ;
1124
/* remove the failed pivcol [OUT] from candidate set */
1125
ASSERT (pivcol [OUT] == Work->Candidates [jcand [OUT]]) ;
1126
remove_candidate (jcand [OUT], Work, Symbolic) ;
1129
/* delete all of the tuples, and all contributions to this column */
1130
DEBUG1 (("Prune tuples of dead outcol: "ID"\n", pivcol [OUT])) ;
1131
Col_tlen [pivcol [OUT]] = 0 ;
1132
UMF_mem_free_tail_block (Numeric, Col_tuples [pivcol [OUT]]) ;
1133
Col_tuples [pivcol [OUT]] = 0 ;
1135
/* no pivot found at all */
1136
return (UMFPACK_WARNING_singular_matrix) ;
1139
/* ------------------------------------------------------------------ */
1143
/* the "in" row is the same as the "in" row for the "in" column */
1145
rdeg [OUT][IN] = rdeg [IN][IN] ;
1146
DEBUG4 (("Freebie in, row "ID"\n", pivrow [IN][IN])) ;
1151
/* the "out" row is the same as the "out" row for the "in" column */
1153
rdeg [OUT][OUT] = rdeg [IN][OUT] ;
1154
DEBUG4 (("Freebie out, row "ID"\n", pivrow [IN][OUT])) ;
1157
/* ------------------------------------------------------------------ */
1158
/* evaluate OUT_IN option */
1159
/* ------------------------------------------------------------------ */
1161
if (pivrow [OUT][IN] != EMPTY)
1163
/* The current front would become an Lson of the new front.
1164
* The candidate pivot row is in the current front, but the
1165
* candidate pivot column is not. */
1167
ASSERT (fnrows > 0 && fncols >= 0) ;
1169
did_rowmerge = (cdeg_out == 0) ;
1172
/* pivrow [OUT][IN] was found via row merge search */
1173
/* it is not (yet) in the pivot column pattern (add it now) */
1174
DEBUGm4 (("did row merge OUT col, IN row\n")) ;
1175
Wm [0] = pivrow [OUT][IN] ;
1178
ASSERT (nr_out == EMPTY) ;
1181
nc = rdeg [OUT][IN] - fncols ;
1184
/* count rows not in current front */
1189
for (i = 0 ; i < cdeg_out ; i++)
1192
ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1193
if (Frpos [row] < 0 || Frpos [row] >= fnrows) nr_out++ ;
1195
/* we must see the pivot row somewhere */
1196
if (row == pivrow [OUT][IN])
1198
ASSERT (Frpos [row] >= 0) ;
1206
/* each column in front grows by nr_out: */
1208
/* new cols not affected by front: */
1209
((nc-1) * (cdeg_out-1)) ;
1211
/* check the cost of relaxed OUT_IN amalgamation */
1213
extra_rows = ((fnrows-1) + nr_out) - (cdeg_out - 1) ;
1214
ASSERT (extra_rows >= 0) ;
1215
ASSERT (fnrows + nr_out == extra_rows + cdeg_out) ;
1216
extra_zeros = (nc-1) * extra_rows ; /* symbolic fill-in */
1218
ASSERT (fnrows + nr_out == cdeg_out + extra_rows) ;
1219
ASSERT (fncols + nc == rdeg [OUT][IN]) ;
1221
/* size of relaxed front (after pivot row removed): */
1222
fnrows_new [OUT][IN] = (fnrows-1) + nr_out ;
1223
fncols_new [OUT][IN] = fncols + (nc-1) ;
1224
relaxed_front = fnrows_new [OUT][IN] * fncols_new [OUT][IN] ;
1226
/* do relaxed amalgamation if the extra zeros are no more */
1227
/* than a fraction (default 0.25) of the relaxed front */
1228
/* if relax = 0: no extra zeros allowed */
1229
/* if relax = +inf: always amalgamate */
1236
/* relax parameter uses a double relop, but ignore NaN case: */
1237
if (extra_zeros == 0)
1243
do_extend = ((double) extra_zeros) <
1244
(relax1 * (double) relaxed_front) ;
1250
/* count the cost of relaxed amalgamation */
1251
thiscost += extra_zeros ;
1253
DEBUG2 (("Evaluating option OUT-IN:\n")) ;
1254
DEBUG2 ((" Work->fnzeros "ID" fnpiv "ID" nr_out "ID" nc "ID"\n",
1255
Work->fnzeros, fnpiv, nr_out, nc)) ;
1256
DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
1258
/* determine if BLAS-3 update to be applied before extending. */
1259
/* update if too many zero entries accumulate in the LU block */
1260
fnzeros = Work->fnzeros + fnpiv * (nr_out + nc) ;
1262
DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
1264
new_LUsize = (fnpiv+1) * (fnrows + nr_out + fncols + nc) ;
1266
DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
1268
/* RELAX3 parameter uses a double relop, ignore NaN case: */
1269
do_update = fnpiv > 0 &&
1270
(((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
1271
DEBUG2 (("do_update "ID"\n", do_update))
1275
/* the current front would not be extended */
1276
do_update = fnpiv > 0 ;
1278
DEBUG2 (("OUT-IN do_update forced true: "ID"\n", do_update)) ;
1280
/* The new front would be just big enough to hold the new
1281
* pivot row and column. */
1282
fnrows_new [OUT][IN] = cdeg_out - 1 ;
1283
fncols_new [OUT][IN] = rdeg [OUT][IN] - 1 ;
1286
DEBUG2 (("option OUT IN : nr "ID" nc "ID" cost "ID"("ID") relax "ID
1287
"\n", nr_out, nc, thiscost, extra_zeros, do_extend)) ;
1289
if (bestcost == EMPTY || thiscost < bestcost)
1291
/* this is the best option seen so far */
1292
Work->pivot_case = OUT_IN ;
1293
bestcost = thiscost ;
1294
Work->do_extend = do_extend ;
1295
Work->do_update = do_update ;
1296
new_fnzeros = fnzeros ;
1300
/* ------------------------------------------------------------------ */
1301
/* evaluate OUT_OUT option */
1302
/* ------------------------------------------------------------------ */
1304
if (pivrow [OUT][OUT] != EMPTY)
1306
/* Neither the candidate pivot row nor the candidate pivot column
1307
* are in the current front. */
1309
ASSERT (fnrows >= 0 && fncols >= 0) ;
1311
did_rowmerge = (cdeg_out == 0) ;
1314
/* pivrow [OUT][OUT] was found via row merge search */
1315
/* it is not (yet) in the pivot column pattern (add it now) */
1316
DEBUGm4 (("did row merge OUT col, OUT row\n")) ;
1317
Wm [0] = pivrow [OUT][OUT] ;
1320
ASSERT (nr_out == EMPTY) ;
1324
if (fnrows == 0 && fncols == 0)
1326
/* the current front is completely empty */
1327
ASSERT (fnpiv == 0) ;
1328
nc = rdeg [OUT][OUT] ;
1334
thiscost = (nc-1) * (cdeg_out-1) ; /* new columns only */
1336
/* size of new front: */
1337
fnrows_new [OUT][OUT] = nr_out-1 ;
1338
fncols_new [OUT][OUT] = nc-1 ;
1339
relaxed_front = fnrows_new [OUT][OUT] * fncols_new [OUT][OUT] ;
1344
/* count rows not in current front */
1345
if (nr_out == EMPTY)
1351
for (i = 0 ; i < cdeg_out ; i++)
1354
ASSERT (row >= 0 && row < n_row) ;
1355
ASSERT (NON_PIVOTAL_ROW (row)) ;
1356
if (Frpos [row] < 0 || Frpos [row] >= fnrows) nr_out++ ;
1358
/* we must see the pivot row somewhere */
1359
if (row == pivrow [OUT][OUT])
1361
ASSERT (Frpos [row] < 0 || Frpos [row] >= fnrows) ;
1369
/* count columns not in current front */
1374
for (i = 0 ; i < rdeg [OUT][OUT] ; i++)
1377
ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1378
if (Fcpos [col] < 0) nc++ ;
1380
/* we must see the pivot column somewhere */
1381
if (col == pivcol [OUT])
1383
ASSERT (Fcpos [col] < 0) ;
1390
extra_cols = (fncols + (nc-1)) - (rdeg [OUT][OUT] - 1) ;
1391
extra_rows = (fnrows + (nr_out-1)) - (cdeg_out - 1) ;
1392
ASSERT (extra_rows >= 0) ;
1393
ASSERT (extra_cols >= 0) ;
1394
extra_zeros = ((nc-1) * extra_rows) + ((nr_out-1) * extra_cols);
1396
ASSERT (fnrows + nr_out == cdeg_out + extra_rows) ;
1397
ASSERT (fncols + nc == rdeg [OUT][OUT] + extra_cols) ;
1401
((nc-1) * (cdeg_out-1)) +
1402
/* old columns in front grow by nr_out-1: */
1403
((nr_out-1) * (fncols - extra_cols)) ;
1405
/* size of relaxed front: */
1406
fnrows_new [OUT][OUT] = fnrows + (nr_out-1) ;
1407
fncols_new [OUT][OUT] = fncols + (nc-1) ;
1408
relaxed_front = fnrows_new [OUT][OUT] * fncols_new [OUT][OUT] ;
1412
/* do relaxed amalgamation if the extra zeros are no more */
1413
/* than a fraction (default 0.25) of the relaxed front */
1414
/* if relax = 0: no extra zeros allowed */
1415
/* if relax = +inf: always amalgamate */
1422
/* relax parameter uses a double relop, but ignore NaN case: */
1423
if (extra_zeros == 0)
1429
do_extend = ((double) extra_zeros) <
1430
(relax1 * (double) relaxed_front) ;
1436
/* count the cost of relaxed amalgamation */
1437
thiscost += extra_zeros ;
1439
DEBUG2 (("Evaluating option OUT-OUT:\n")) ;
1440
DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_out "ID" nc "ID"\n",
1441
Work->fnzeros, fnpiv, nr_out, nc)) ;
1442
DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
1444
/* determine if BLAS-3 update to be applied before extending. */
1445
/* update if too many zero entries accumulate in the LU block */
1446
fnzeros = Work->fnzeros + fnpiv * (nr_out + nc) ;
1448
DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
1450
new_LUsize = (fnpiv+1) * (fnrows + nr_out + fncols + nc) ;
1452
DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
1454
/* RELAX3 parameter uses a double relop, ignore NaN case: */
1455
do_update = fnpiv > 0 &&
1456
(((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
1457
DEBUG2 (("do_update "ID"\n", do_update))
1461
/* the current front would not be extended */
1462
do_update = fnpiv > 0 ;
1464
DEBUG2 (("OUT-OUT do_update forced true: "ID"\n", do_update)) ;
1466
/* The new front would be just big enough to hold the new
1467
* pivot row and column. */
1468
fnrows_new [OUT][OUT] = cdeg_out - 1 ;
1469
fncols_new [OUT][OUT] = rdeg [OUT][OUT] - 1 ;
1472
DEBUG2 (("option OUT OUT: nr "ID" nc "ID" cost "ID"\n",
1473
rdeg [OUT][OUT], cdeg_out, thiscost)) ;
1475
if (bestcost == EMPTY || thiscost < bestcost)
1477
/* this is the best option seen so far */
1478
Work->pivot_case = OUT_OUT ;
1479
bestcost = thiscost ;
1480
Work->do_extend = do_extend ;
1481
Work->do_update = do_update ;
1482
new_fnzeros = fnzeros ;
1487
/* At this point, a structural pivot has been found. */
1488
/* It may be numerically zero, however. */
1489
ASSERT (Work->pivot_case != EMPTY) ;
1490
DEBUG2 (("local search, best option "ID", best cost "ID"\n",
1491
Work->pivot_case, bestcost)) ;
1493
/* ====================================================================== */
1494
/* Pivot row and column, and extension, now determined */
1495
/* ====================================================================== */
1497
Work->fnzeros = new_fnzeros ;
1499
/* ---------------------------------------------------------------------- */
1500
/* finalize the pivot row and column */
1501
/* ---------------------------------------------------------------------- */
1503
switch (Work->pivot_case)
1506
DEBUG2 (("IN-IN option selected\n")) ;
1507
ASSERT (fnrows > 0 && fncols > 0) ;
1508
Work->pivcol_in_front = TRUE ;
1509
Work->pivrow_in_front = TRUE ;
1510
Work->pivcol = pivcol [IN] ;
1511
Work->pivrow = pivrow [IN][IN] ;
1512
Work->ccdeg = nr_in ;
1513
Work->Wrow = Fcols ;
1514
Work->rrdeg = rdeg [IN][IN] ;
1516
Work->fnrows_new = fnrows_new [IN][IN] ;
1517
Work->fncols_new = fncols_new [IN][IN] ;
1521
DEBUG2 (("IN-OUT option selected\n")) ;
1522
ASSERT (fnrows >= 0 && fncols > 0) ;
1523
Work->pivcol_in_front = TRUE ;
1524
Work->pivrow_in_front = FALSE ;
1525
Work->pivcol = pivcol [IN] ;
1526
Work->pivrow = pivrow [IN][OUT] ;
1527
Work->ccdeg = nr_in ;
1529
Work->rrdeg = rdeg [IN][OUT] ;
1531
Work->fnrows_new = fnrows_new [IN][OUT] ;
1532
Work->fncols_new = fncols_new [IN][OUT] ;
1536
DEBUG2 (("OUT-IN option selected\n")) ;
1537
ASSERT (fnrows > 0 && fncols >= 0) ;
1538
Work->pivcol_in_front = FALSE ;
1539
Work->pivrow_in_front = TRUE ;
1540
Work->pivcol = pivcol [OUT] ;
1541
Work->pivrow = pivrow [OUT][IN] ;
1542
Work->ccdeg = cdeg_out ;
1543
/* Wrow might be equivalenced to Fcols (Freebie in): */
1545
Work->rrdeg = rdeg [OUT][IN] ;
1546
/* Work->Wrow[0..fncols-1] is not there. See Fcols instead */
1548
Work->fnrows_new = fnrows_new [OUT][IN] ;
1549
Work->fncols_new = fncols_new [OUT][IN] ;
1553
DEBUG2 (("OUT-OUT option selected\n")) ;
1554
ASSERT (fnrows >= 0 && fncols >= 0) ;
1555
Work->pivcol_in_front = FALSE ;
1556
Work->pivrow_in_front = FALSE ;
1557
Work->pivcol = pivcol [OUT] ;
1558
Work->pivrow = pivrow [OUT][OUT] ;
1559
Work->ccdeg = cdeg_out ;
1560
/* Wrow might be equivalenced to Wio (Freebie out): */
1562
Work->rrdeg = rdeg [OUT][OUT] ;
1564
Work->fnrows_new = fnrows_new [OUT][OUT] ;
1565
Work->fncols_new = fncols_new [OUT][OUT] ;
1570
ASSERT (IMPLIES (fnrows == 0 && fncols == 0, Work->pivot_case == OUT_OUT)) ;
1572
if (!Work->pivcol_in_front && pivcol [IN] != EMPTY)
1574
/* clear Frpos if pivcol [IN] was searched, but not selected */
1575
for (i = fnrows ; i < cdeg_in ; i++)
1577
Frpos [Frows [i]] = EMPTY;
1581
/* ---------------------------------------------------------------------- */
1582
/* Pivot row and column have been found */
1583
/* ---------------------------------------------------------------------- */
1585
/* ---------------------------------------------------------------------- */
1586
/* remove pivot column from candidate pivot column set */
1587
/* ---------------------------------------------------------------------- */
1589
ASSERT (jj >= 0 && jj < Work->nCandidates) ;
1590
ASSERT (Work->pivcol == Work->Candidates [jj]) ;
1591
remove_candidate (jj, Work, Symbolic) ;
1593
/* ---------------------------------------------------------------------- */
1594
/* check for frontal matrix growth */
1595
/* ---------------------------------------------------------------------- */
1597
DEBUG1 (("Check frontal growth:\n")) ;
1598
DEBUG1 (("fnrows_new "ID" + 1 = "ID", fnr_curr "ID"\n",
1599
Work->fnrows_new, Work->fnrows_new + 1, fnr_curr)) ;
1600
DEBUG1 (("fncols_new "ID" + 1 = "ID", fnc_curr "ID"\n",
1601
Work->fncols_new, Work->fncols_new + 1, fnc_curr)) ;
1603
Work->do_grow = (Work->fnrows_new + 1 > fnr_curr
1604
|| Work->fncols_new + 1 > fnc_curr) ;
1607
DEBUG0 (("\nNeed to grow frontal matrix, force do_update true\n")) ;
1608
/* If the front must grow, then apply the pending updates and remove
1609
* the current pivot rows/columns from the front prior to growing the
1610
* front. This frees up as much space as possible for the new front. */
1611
if (!Work->do_update && fnpiv > 0)
1613
/* This update would not have to be done if the current front
1614
* was big enough. */
1616
Work->do_update = TRUE ;
1620
/* ---------------------------------------------------------------------- */
1621
/* current pivot column */
1622
/* ---------------------------------------------------------------------- */
1625
c1) If pivot column index is in the current front:
1627
The pivot column pattern is in Frows [0 .. fnrows-1] and
1628
the extension is in Frows [fnrows ... fnrows+ccdeg-1].
1630
Frpos [Frows [0 .. fnrows+ccdeg-1]] is
1631
equal to 0 .. fnrows+ccdeg-1. Wm is not needed.
1633
The values are in Wy [0 .. fnrows+ccdeg-1].
1635
c2) Otherwise, if the pivot column index is not in the current front:
1637
c2a) If the front is being extended, old row indices in the the
1638
pivot column pattern are in Frows [0 .. fnrows-1].
1640
All entries are in Wm [0 ... ccdeg-1], with values in
1641
Wx [0 .. ccdeg-1]. These may include entries already in
1642
Frows [0 .. fnrows-1].
1644
Frpos [Frows [0 .. fnrows-1]] is equal to 0 .. fnrows-1.
1645
Frpos [Wm [0 .. ccdeg-1]] for new entries is < 0.
1647
c2b) If the front is not being extended, then the entire pivot
1648
column pattern is in Wm [0 .. ccdeg-1]. It includes
1649
the pivot row index. It is does not contain the pattern
1650
Frows [0..fnrows-1]. The intersection of these two
1651
sets may or may not be empty. The values are in Wx [0..ccdeg-1]
1653
In both cases c1 and c2, Frpos [Frows [0 .. fnrows-1]] is equal
1654
to 0 .. fnrows-1, which is the pattern of the current front.
1655
Any entry of Frpos that is not specified above is < 0.
1660
DEBUG2 (("\n\nSEARCH DONE: Pivot col "ID" in: ("ID") pivot row "ID" in: ("ID
1661
") extend: "ID"\n\n", Work->pivcol, Work->pivcol_in_front,
1662
Work->pivrow, Work->pivrow_in_front, Work->do_extend)) ;
1663
UMF_dump_rowcol (1, Numeric, Work, Work->pivcol, !Symbolic->fixQ) ;
1664
DEBUG2 (("Pivot col "ID": fnrows "ID" ccdeg "ID"\n", Work->pivcol, fnrows,
1666
if (Work->pivcol_in_front) /* case c1 */
1669
DEBUG3 (("Pivcol in front\n")) ;
1670
for (i = 0 ; i < fnrows ; i++)
1673
DEBUG3 ((ID": row:: "ID" in front ", i, row)) ;
1674
ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1675
ASSERT (Frpos [row] == i) ;
1677
if (row == Work->pivrow)
1679
DEBUG3 ((" <- pivrow")) ;
1684
ASSERT (found == Work->pivrow_in_front) ;
1686
for (i = fnrows ; i < fnrows + Work->ccdeg ; i++)
1689
DEBUG3 ((ID": row:: "ID" (new)", i, row)) ;
1690
ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1691
ASSERT (Frpos [row] == i) ;
1693
if (row == Work->pivrow)
1695
DEBUG3 ((" <- pivrow")) ;
1700
ASSERT (found == !Work->pivrow_in_front) ;
1704
if (Work->do_extend)
1707
DEBUG3 (("Pivcol not in front (extend)\n")) ;
1708
for (i = 0 ; i < fnrows ; i++)
1711
DEBUG3 ((ID": row:: "ID" in front ", i, row)) ;
1712
ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1713
ASSERT (Frpos [row] == i) ;
1714
if (row == Work->pivrow)
1716
DEBUG3 ((" <- pivrow")) ;
1721
ASSERT (found == Work->pivrow_in_front) ;
1723
DEBUG3 (("----\n")) ;
1724
for (i = 0 ; i < Work->ccdeg ; i++)
1727
ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1728
DEBUG3 ((ID": row:: "ID" ", i, row)) ;
1730
if (Frpos [row] < 0)
1732
DEBUG3 ((" (new) ")) ;
1734
if (row == Work->pivrow)
1736
DEBUG3 ((" <- pivrow")) ;
1739
if (Work->pivrow_in_front) ASSERT (Frpos [row] >= 0) ;
1740
else ASSERT (Frpos [row] < 0) ;
1749
DEBUG3 (("Pivcol not in front (no extend)\n")) ;
1750
for (i = 0 ; i < Work->ccdeg ; i++)
1753
ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
1754
DEBUG3 ((ID": row:: "ID" ", i, row)) ;
1756
DEBUG3 ((" (new) ")) ;
1757
if (row == Work->pivrow)
1759
DEBUG3 ((" <- pivrow")) ;
1769
/* ---------------------------------------------------------------------- */
1770
/* current pivot row */
1771
/* ---------------------------------------------------------------------- */
1774
r1) If the pivot row index is in the current front:
1776
The pivot row pattern is in Fcols [0..fncols-1] and the extenson is
1777
in Wrow [fncols .. rrdeg-1]. If the pivot column is in the current
1778
front, then Fcols and Wrow are equivalenced.
1780
r2) If the pivot row index is not in the current front:
1782
r2a) If the front is being extended, the pivot row pattern is in
1783
Fcols [0 .. fncols-1]. New entries are in Wrow [0 .. rrdeg-1],
1784
but these may include entries already in Fcols [0 .. fncols-1].
1786
r2b) Otherwise, the pivot row pattern is Wrow [0 .. rrdeg-1].
1788
Fcpos [Fcols [0..fncols-1]] is (0..fncols-1) * fnr_curr.
1789
All other entries in Fcpos are < 0.
1791
These conditions are asserted below.
1793
------------------------------------------------------------------------
1794
Other items in Work structure that are relevant:
1796
pivcol: the pivot column index
1797
pivrow: the pivot column index
1802
fnrows: the number of rows in the currnt contribution block
1803
fncols: the number of columns in the current contribution block
1805
fnrows_new: the number of rows in the new contribution block
1806
fncols_new: the number of rows in the new contribution block
1808
------------------------------------------------------------------------
1813
UMF_dump_rowcol (0, Numeric, Work, Work->pivrow, TRUE) ;
1814
DEBUG2 (("Pivot row "ID":\n", Work->pivrow)) ;
1815
if (Work->pivrow_in_front)
1818
for (i = 0 ; i < fncols ; i++)
1821
DEBUG3 ((" col:: "ID" in front\n", col)) ;
1822
ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1823
ASSERT (Fcpos [col] == i * fnr_curr) ;
1824
if (col == Work->pivcol) found = TRUE ;
1826
ASSERT (found == Work->pivcol_in_front) ;
1828
ASSERT (IMPLIES (Work->pivcol_in_front, Fcols == Work->Wrow)) ;
1829
for (i = fncols ; i < Work->rrdeg ; i++)
1831
col = Work->Wrow [i] ;
1832
ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1833
ASSERT (Fcpos [col] < 0) ;
1834
if (col == Work->pivcol) found = TRUE ;
1835
else DEBUG3 ((" col:: "ID" (new)\n", col)) ;
1837
ASSERT (found == !Work->pivcol_in_front) ;
1841
if (Work->do_extend)
1844
for (i = 0 ; i < fncols ; i++)
1847
DEBUG3 ((" col:: "ID" in front\n", col)) ;
1848
ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1849
ASSERT (Fcpos [col] == i * fnr_curr) ;
1850
if (col == Work->pivcol) found = TRUE ;
1852
ASSERT (found == Work->pivcol_in_front) ;
1854
for (i = 0 ; i < Work->rrdeg ; i++)
1856
col = Work->Wrow [i] ;
1857
ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1858
if (Fcpos [col] >= 0) continue ;
1859
if (col == Work->pivcol) found = TRUE ;
1860
else DEBUG3 ((" col:: "ID" (new, extend)\n", col)) ;
1862
ASSERT (found == !Work->pivcol_in_front) ;
1867
for (i = 0 ; i < Work->rrdeg ; i++)
1869
col = Work->Wrow [i] ;
1870
ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
1871
if (col == Work->pivcol) found = TRUE ;
1872
else DEBUG3 ((" col:: "ID" (all new)\n", col)) ;
1879
/* ---------------------------------------------------------------------- */
1880
/* determine whether to do scan2-row and scan2-col */
1881
/* ---------------------------------------------------------------------- */
1883
if (Work->do_extend)
1885
Work->do_scan2row = (fncols > 0) ;
1886
Work->do_scan2col = (fnrows > 0) ;
1890
Work->do_scan2row = (fncols > 0) && Work->pivrow_in_front ;
1891
Work->do_scan2col = (fnrows > 0) && Work->pivcol_in_front ;
1894
/* ---------------------------------------------------------------------- */
1896
DEBUG2 (("LOCAL SEARCH DONE: pivot column "ID" pivot row: "ID"\n",
1897
Work->pivcol, Work->pivrow)) ;
1898
DEBUG2 (("do_extend: "ID"\n", Work->do_extend)) ;
1899
DEBUG2 (("do_update: "ID"\n", Work->do_update)) ;
1900
DEBUG2 (("do_grow: "ID"\n", Work->do_grow)) ;
1902
/* ---------------------------------------------------------------------- */
1903
/* keep track of the diagonal */
1904
/* ---------------------------------------------------------------------- */
1906
if (Symbolic->prefer_diagonal
1907
&& Work->pivcol < Work->n_col - Symbolic->nempty_col)
1909
Diagonal_map = Work->Diagonal_map ;
1910
Diagonal_imap = Work->Diagonal_imap ;
1911
ASSERT (Diagonal_map != (Int *) NULL) ;
1912
ASSERT (Diagonal_imap != (Int *) NULL) ;
1914
row2 = Diagonal_map [Work->pivcol] ;
1915
col2 = Diagonal_imap [Work->pivrow] ;
1919
/* this was an off-diagonal pivot row */
1920
Work->noff_diagonal++ ;
1921
row2 = UNFLIP (row2) ;
1924
ASSERT (Diagonal_imap [row2] == Work->pivcol) ;
1925
ASSERT (UNFLIP (Diagonal_map [col2]) == Work->pivrow) ;
1927
if (row2 != Work->pivrow)
1929
/* swap the diagonal map to attempt to maintain symmetry later on.
1930
* Also mark the map for col2 (via FLIP) to denote that the entry
1931
* now on the diagonal is not the original entry on the diagonal. */
1933
DEBUG0 (("Unsymmetric pivot\n")) ;
1934
Diagonal_map [Work->pivcol] = FLIP (Work->pivrow) ;
1935
Diagonal_imap [Work->pivrow] = Work->pivcol ;
1937
Diagonal_map [col2] = FLIP (row2) ;
1938
Diagonal_imap [row2] = col2 ;
1941
ASSERT (n_row == n_col) ;
1943
UMF_dump_diagonal_map (Diagonal_map, Diagonal_imap, Symbolic->n1,
1944
Symbolic->n_col, Symbolic->nempty_col) ;
1948
return (UMFPACK_OK) ;
added
removed
Lib/sandbox/pysparse/umfpack/umf_local_search.c
