00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #include "umf_internal.h"
00021 #include "umf_row_search.h"
00022 #include "umf_mem_free_tail_block.h"
00023
00024
00025 #define RELAX1 0.25
00026 #define SYM_RELAX1 0.0
00027 #define RELAX2 0.1
00028 #define RELAX3 0.125
00029
00030
00031
00032
00033
00034
00035
00036 PRIVATE void remove_candidate (Int jj, WorkType *Work, SymbolicType *Symbolic)
00037 {
00038
00039 #ifndef NDEBUG
00040 Int j ;
00041 DEBUGm2 (("pivot column Candidates before remove: nCand "ID" ncand "ID
00042 " lo "ID" hi "ID" jj "ID"\n", Work->nCandidates, Work->ncand,
00043 Work->lo, Work->hi, jj)) ;
00044 for (j = 0 ; j < Work->nCandidates ; j++)
00045 {
00046 Int col = Work->Candidates [j] ;
00047 DEBUGm2 ((ID" ", col));
00048 ASSERT (col >= 0 && col < Work->n_col) ;
00049
00050 ASSERT (col >= Work->lo && col <= Work->hi) ;
00051 }
00052 DEBUGm2 (("\n")) ;
00053 #endif
00054
00055 if (Symbolic->fixQ)
00056 {
00057 DEBUGm2 (("FixQ\n")) ;
00058
00059 ASSERT (Work->nCandidates == 1) ;
00060 ASSERT (jj == 0) ;
00061 if (Work->ncand > 1)
00062 {
00063 Work->Candidates [0] = Work->nextcand++ ;
00064 }
00065 else
00066 {
00067 Work->nCandidates = 0 ;
00068 }
00069 }
00070 else
00071 {
00072
00073 if (Work->ncand > MAX_CANDIDATES)
00074 {
00075 Work->Candidates [jj] = Work->nextcand++ ;
00076 }
00077 else
00078 {
00079 ASSERT (Work->nCandidates == Work->ncand) ;
00080 Work->Candidates [jj] = Work->Candidates [Work->ncand - 1] ;
00081 Work->Candidates [Work->ncand - 1] = EMPTY ;
00082 Work->nCandidates-- ;
00083 }
00084 }
00085 Work->ncand-- ;
00086
00087 #ifndef NDEBUG
00088 DEBUGm2 (("pivot column Candidates after remove: nCand "ID" ncand "ID
00089 " lo "ID" hi "ID" jj "ID"\n", Work->nCandidates, Work->ncand, Work->lo,
00090 Work->hi, jj)) ;
00091 for (j = 0 ; j < Work->nCandidates ; j++)
00092 {
00093 Int col = Work->Candidates [j] ;
00094 DEBUGm2 ((ID" ", col));
00095 ASSERT (col >= 0 && col < Work->n_col) ;
00096
00097 ASSERT (col >= Work->lo && col <= Work->hi) ;
00098 }
00099 DEBUGm2 (("\n")) ;
00100 ASSERT (Work->ncand >= 0) ;
00101 ASSERT (Work->nCandidates <= Work->ncand) ;
00102 #endif
00103 }
00104
00105
00106
00107
00108
00109 GLOBAL Int UMF_local_search
00110 (
00111 NumericType *Numeric,
00112 WorkType *Work,
00113 SymbolicType *Symbolic
00114 )
00115 {
00116
00117
00118
00119
00120 double relax1 ;
00121 Entry *Flblock, *Fublock, *Fs, *Fcblock, *C, *Wx, *Wy, *Fu, *Flublock,
00122 *Flu ;
00123 Int pos, nrows, *Cols, *Rows, e, f, status, max_cdeg, fnzeros, nb, j, col,
00124 i, row, cdeg_in, rdeg [2][2], fnpiv, nothing [2], new_LUsize,
00125 pivrow [2][2], pivcol [2], *Wp, *Fcpos, *Frpos, new_fnzeros, cdeg_out,
00126 *Wm, *Wio, *Woi, *Woo, *Frows, *Fcols, fnrows, fncols, *E, deg, nr_in,
00127 nc, thiscost, bestcost, nr_out, do_update, extra_cols, extra_rows,
00128 extra_zeros, relaxed_front, do_extend, fnr_curr, fnc_curr, tpi,
00129 *Col_tuples, *Col_degree, *Col_tlen, jj, jcand [2], freebie [2],
00130 did_rowmerge, fnrows_new [2][2], fncols_new [2][2], search_pivcol_out,
00131 *Diagonal_map, *Diagonal_imap, row2, col2 ;
00132 Unit *Memory, *p ;
00133 Tuple *tp, *tpend, *tp1, *tp2 ;
00134 Element *ep ;
00135
00136 #ifndef NBLAS
00137 Int blas_ok = TRUE ;
00138 #else
00139 #define blas_ok FALSE
00140 #endif
00141
00142 #ifndef NDEBUG
00143 Int debug_ok, n_row, n_col, *Row_degree ;
00144 Row_degree = Numeric->Rperm ;
00145 #endif
00146
00147
00148
00149
00150
00151 Memory = Numeric->Memory ;
00152 E = Work->E ;
00153 Col_degree = Numeric->Cperm ;
00154
00155 Col_tuples = Numeric->Lip ;
00156 Col_tlen = Numeric->Lilen ;
00157
00158 Wx = Work->Wx ;
00159 Wy = Work->Wy ;
00160 Wp = Work->Wp ;
00161 Wm = Work->Wm ;
00162 Woi = Work->Woi ;
00163 Wio = Work->Wio ;
00164 Woo = Work->Woo ;
00165 Fcpos = Work->Fcpos ;
00166 Frpos = Work->Frpos ;
00167 Frows = Work->Frows ;
00168 Fcols = Work->Fcols ;
00169 fnrows = Work->fnrows ;
00170 fncols = Work->fncols ;
00171 nb = Work->nb ;
00172 fnr_curr = Work->fnr_curr ;
00173 fnc_curr = Work->fnc_curr ;
00174 fnpiv = Work->fnpiv ;
00175 nothing [0] = EMPTY ;
00176 nothing [1] = EMPTY ;
00177 relax1 = (Symbolic->prefer_diagonal) ? SYM_RELAX1 : RELAX1 ;
00178 fnzeros = Work->fnzeros ;
00179 new_fnzeros = fnzeros ;
00180 jj = EMPTY ;
00181
00182 Fcblock = Work->Fcblock ;
00183 Flblock = Work->Flblock ;
00184 Fublock = Work->Fublock ;
00185 Flublock = Work->Flublock ;
00186
00187
00188 max_cdeg = Work->fnrows_max ;
00189 ASSERT (Work->fnrows_max <= Symbolic->maxnrows) ;
00190 ASSERT (Work->fncols_max <= Symbolic->maxncols) ;
00191
00192 if (fnrows == 0 && fncols == 0)
00193 {
00194
00195 Work->firstsuper = Work->ksuper ;
00196 }
00197
00198 #ifndef NDEBUG
00199 n_row = Work->n_row ;
00200 n_col = Work->n_col ;
00201 DEBUG2 (("\n========LOCAL SEARCH: current frontal matrix: ========= \n")) ;
00202 UMF_dump_current_front (Numeric, Work, TRUE) ;
00203 if (UMF_debug > 0 || MAX (n_row, n_col) < 1000)
00204 {
00205 for (i = 0 ; i < MAX (n_row, n_col) ; i++)
00206 {
00207 ASSERT (Wp [i] < 0) ;
00208 }
00209 }
00210
00211 DEBUGm2 ((ID" pivot column Candidates: lo "ID" hi "ID"\n",
00212 Work->nCandidates, Work->lo, Work->hi)) ;
00213 for (j = 0 ; j < Work->nCandidates ; j++)
00214 {
00215 col = Work->Candidates [j] ;
00216 DEBUGm2 ((ID" ", col));
00217 ASSERT (col >= 0 && col < n_col) ;
00218 ASSERT (NON_PIVOTAL_COL (col)) ;
00219 ASSERT (col >= Work->lo && col <= Work->hi) ;
00220 }
00221
00222 DEBUGm2 (("\n")) ;
00223
00224
00225 DEBUG2 (("\n\n::: "ID" : Npiv: "ID" + fnpiv "ID" = "ID". "
00226 "size "ID"-by-"ID"\n", Work->frontid,
00227 Work->npiv, Work->fnpiv, Work->npiv + Work->fnpiv, fnrows, fncols)) ;
00228 ASSERT ((fnrows == 0 && fncols == 0) ||(fnrows != 0 && fncols != 0)) ;
00229 #endif
00230
00231
00232
00233
00234
00235
00236
00237 pivcol [IN] = EMPTY ;
00238 pivcol [OUT] = EMPTY ;
00239
00240 cdeg_in = Int_MAX ;
00241 cdeg_out = Int_MAX ;
00242
00243 pivrow [IN][IN] = EMPTY ;
00244 pivrow [IN][OUT] = EMPTY ;
00245 pivrow [OUT][IN] = EMPTY ;
00246 pivrow [OUT][OUT] = EMPTY ;
00247
00248 rdeg [IN][IN] = Int_MAX ;
00249 rdeg [IN][OUT] = Int_MAX ;
00250 rdeg [OUT][IN] = Int_MAX ;
00251 rdeg [OUT][OUT] = Int_MAX ;
00252
00253 freebie [IN] = FALSE ;
00254 freebie [OUT] = FALSE ;
00255
00256 Work->pivot_case = EMPTY ;
00257 bestcost = EMPTY ;
00258
00259 nr_out = EMPTY ;
00260 nr_in = EMPTY ;
00261
00262 jcand [IN] = EMPTY ;
00263 jcand [OUT] = EMPTY ;
00264
00265 fnrows_new [IN][IN] = EMPTY ;
00266 fnrows_new [IN][OUT] = EMPTY ;
00267 fnrows_new [OUT][IN] = EMPTY ;
00268 fnrows_new [OUT][OUT] = EMPTY ;
00269
00270 fncols_new [IN][IN] = EMPTY ;
00271 fncols_new [IN][OUT] = EMPTY ;
00272 fncols_new [OUT][IN] = EMPTY ;
00273 fncols_new [OUT][OUT] = EMPTY ;
00274
00275 #ifndef NDEBUG
00276
00277 DEBUG4 (("Check Frpos : fnrows "ID" col "ID" maxcdeg "ID"\n",
00278 fnrows, pivcol [IN], max_cdeg)) ;
00279 for (i = 0 ; i < fnrows ; i++)
00280 {
00281 row = Frows [i] ;
00282 DEBUG4 ((" row: "ID"\n", row)) ;
00283 ASSERT (row >= 0 && row < n_row) ;
00284 ASSERT (Frpos [row] == i) ;
00285 }
00286 DEBUG4 (("All:\n")) ;
00287 if (UMF_debug > 0 || n_row < 1000)
00288 {
00289 Int cnt = fnrows ;
00290 for (row = 0 ; row < n_row ; row++)
00291 {
00292 if (Frpos [row] == EMPTY)
00293 {
00294 cnt++ ;
00295 }
00296 else
00297 {
00298 DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
00299 }
00300 }
00301 ASSERT (cnt == n_row) ;
00302 }
00303 #endif
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318 DEBUG2 (("Max candidates %d, Work->ncand "ID" jmax "ID"\n",
00319 MAX_CANDIDATES, Work->ncand, Work->nCandidates)) ;
00320 col = Work->Candidates [0] ;
00321 ASSERT (Work->nCandidates > 0) ;
00322 DEBUG3 (("Pivot column candidate: "ID" j = "ID"\n", col, j)) ;
00323 ASSERT (col >= 0 && col < n_col) ;
00324
00325
00326 deg = Symbolic->fixQ ? EMPTY : Col_degree [col] ;
00327
00328 #ifndef NDEBUG
00329 DEBUG3 (("Pivot column candidate: "ID" cost: "ID" Fcpos[col] "ID"\n",
00330 col, deg, Fcpos [col])) ;
00331 UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
00332 if (Symbolic->fixQ)
00333 {
00334 DEBUG1 (("FIXQ: Candidates "ID" pivcol "ID" npiv "ID" fnpiv "ID
00335 " ndiscard "ID "\n", Work->nCandidates, col, Work->npiv,
00336 Work->fnpiv, Work->ndiscard)) ;
00337 ASSERT (Work->nCandidates == 1) ;
00338 ASSERT (col == Work->npiv + Work->fnpiv + Work->ndiscard) ;
00339 }
00340 #endif
00341
00342 if (Fcpos [col] >= 0)
00343 {
00344
00345 pivcol [IN] = col ;
00346 cdeg_in = deg ;
00347 jcand [IN] = 0 ;
00348 }
00349 else
00350 {
00351
00352 pivcol [OUT] = col ;
00353 cdeg_out = deg ;
00354 jcand [OUT] = 0 ;
00355 }
00356
00357
00358 for (j = 1 ; j < Work->nCandidates ; j++)
00359 {
00360 col = Work->Candidates [j] ;
00361
00362 DEBUG3 (("Pivot col candidate: "ID" j = "ID"\n", col, j)) ;
00363 ASSERT (col >= 0 && col < n_col) ;
00364 ASSERT (!Symbolic->fixQ) ;
00365 deg = Col_degree [col] ;
00366 #ifndef NDEBUG
00367 DEBUG3 (("Pivot col candidate: "ID" cost: "ID" Fcpos[col] "ID"\n",
00368 col, deg, Fcpos [col])) ;
00369 UMF_dump_rowcol (1, Numeric, Work, col, !Symbolic->fixQ) ;
00370 #endif
00371 if (Fcpos [col] >= 0)
00372 {
00373 #ifndef NDEBUG
00374 Int fs ;
00375 fs = Fcpos [col] / fnr_curr ;
00376 ASSERT (fs >= 0 && fs < fncols) ;
00377 #endif
00378 if (deg < cdeg_in || (deg == cdeg_in && col < pivcol [IN]))
00379 {
00380
00381 pivcol [IN] = col ;
00382 cdeg_in = deg ;
00383 jcand [IN] = j ;
00384 }
00385 }
00386 else
00387 {
00388 if (deg < cdeg_out || (deg == cdeg_out && col < pivcol [OUT]))
00389 {
00390
00391 pivcol [OUT] = col ;
00392 cdeg_out = deg ;
00393 jcand [OUT] = j ;
00394 }
00395 }
00396 }
00397
00398 DEBUG2 (("Pivcol in "ID" out "ID"\n", pivcol [IN], pivcol [OUT])) ;
00399 ASSERT ((pivcol [IN] >= 0 && pivcol [IN] < n_col)
00400 || (pivcol [OUT] >= 0 && pivcol [OUT] < n_col)) ;
00401
00402 cdeg_in = EMPTY ;
00403 cdeg_out = EMPTY ;
00404
00405
00406
00407
00408
00409 #ifndef NDEBUG
00410
00411 DEBUG4 (("Prior to col update: fnrows "ID" col "ID" maxcdeg "ID"\n",
00412 fnrows, pivcol [IN], max_cdeg)) ;
00413 for (i = 0 ; i < fnrows ; i++)
00414 {
00415 row = Frows [i] ;
00416 DEBUG4 ((" row: "ID"\n", row)) ;
00417 ASSERT (row >= 0 && row < n_row) ;
00418 ASSERT (Frpos [row] == i) ;
00419 }
00420 DEBUG4 (("All:\n")) ;
00421 if (UMF_debug > 0 || n_row < 1000)
00422 {
00423 Int cnt = fnrows ;
00424 for (row = 0 ; row < n_row ; row++)
00425 {
00426 if (Frpos [row] == EMPTY)
00427 {
00428 cnt++ ;
00429 }
00430 else
00431 {
00432 DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
00433 }
00434 }
00435 ASSERT (cnt == n_row) ;
00436 }
00437 #endif
00438
00439 if (pivcol [IN] != EMPTY)
00440 {
00441
00442 #ifndef NDEBUG
00443 DEBUG2 (("col[IN] column "ID" in front at position = "ID"\n",
00444 pivcol [IN], Fcpos [pivcol [IN]])) ;
00445 UMF_dump_rowcol (1, Numeric, Work, pivcol [IN], !Symbolic->fixQ) ;
00446 #endif
00447
00448
00449 ASSERT (fnrows > 0 && fncols > 0) ;
00450
00451 DEBUG4 (("Update pivot column:\n")) ;
00452 Fs = Fcblock + Fcpos [pivcol [IN]] ;
00453 Fu = Fublock + (Fcpos [pivcol [IN]] / fnr_curr) ;
00454 Flu = Flublock + fnpiv * nb ;
00455
00456
00457
00458
00459
00460
00461 for (i = 0 ; i < fnpiv ; i++)
00462 {
00463 Flu [i] = Fu [i*fnc_curr] ;
00464 }
00465
00466
00467
00468
00469
00470
00471
00472
00473 if (fnpiv > 1)
00474 {
00475
00476
00477 #ifndef NBLAS
00478 BLAS_TRSV (fnpiv, Flublock, Flu, nb) ;
00479 #endif
00480 if (!blas_ok)
00481 {
00482
00483 Entry *Flub = Flublock ;
00484 for (j = 0 ; j < fnpiv ; j++)
00485 {
00486 Entry Fuj = Flu [j] ;
00487 #pragma ivdep
00488 for (i = j+1 ; i < fnpiv ; i++)
00489 {
00490
00491 MULT_SUB (Flu [i], Flub [i], Fuj) ;
00492 }
00493 Flub += nb ;
00494 }
00495 }
00496 }
00497
00498
00499
00500
00501
00502 for (i = 0 ; i < fnrows ; i++)
00503 {
00504 Wy [i] = Fs [i] ;
00505 }
00506
00507
00508
00509
00510
00511
00512
00513 #ifdef NBLAS
00514
00515 for (j = 0 ; j < fnpiv ; j++)
00516 {
00517 Entry Fuj, *Flub = Flblock + j * fnr_curr ;
00518 Fuj = Flu [j] ;
00519 if (IS_NONZERO (Fuj))
00520 {
00521 #pragma ivdep
00522 for (i = 0 ; i < fnrows ; i++)
00523 {
00524
00525 MULT_SUB (Wy [i], Flub [i], Fuj) ;
00526 }
00527 }
00528
00529 }
00530 #else
00531
00532
00533 BLAS_GEMV (fnrows, fnpiv, Flblock, Flu, Wy, fnr_curr) ;
00534 #endif
00535
00536
00537
00538 #ifndef NDEBUG
00539 DEBUG2 (("Wy after update: fnrows="ID"\n", fnrows)) ;
00540 DEBUG4 ((" fnpiv="ID" \n", fnpiv)) ;
00541 for (i = 0 ; i < fnrows ; i++)
00542 {
00543 DEBUG4 ((ID" "ID" "ID, i, Frows [i], Frpos [Frows [i]])) ;
00544 EDEBUG4 (Wy [i]) ;
00545 DEBUG4 (("\n")) ;
00546 }
00547 #endif
00548
00549
00550
00551
00552
00553 cdeg_in = fnrows ;
00554
00555 #ifndef NDEBUG
00556
00557 DEBUG4 (("After col update: fnrows "ID" col "ID" maxcdeg "ID"\n",
00558 fnrows, pivcol [IN], max_cdeg)) ;
00559 for (i = 0 ; i < fnrows ; i++)
00560 {
00561 row = Frows [i] ;
00562 DEBUG4 ((" row: "ID"\n", row)) ;
00563 ASSERT (row >= 0 && row < n_row) ;
00564 ASSERT (Frpos [row] == i) ;
00565 }
00566 DEBUG4 (("All:\n")) ;
00567 if (UMF_debug > 0 || n_row < 1000)
00568 {
00569 Int cnt = fnrows ;
00570 for (row = 0 ; row < n_row ; row++)
00571 {
00572 if (Frpos [row] == EMPTY)
00573 {
00574 cnt++ ;
00575 }
00576 else
00577 {
00578 DEBUG4 ((" row: "ID" pos "ID"\n", row, Frpos [row])) ;
00579 }
00580 }
00581 ASSERT (cnt == n_row) ;
00582 }
00583 #endif
00584
00585 #ifndef NDEBUG
00586
00587 DEBUG4 (("COL ASSEMBLE: cdeg "ID"\nREDUCE COL in "ID" max_cdeg "ID"\n",
00588 cdeg_in, pivcol [IN], max_cdeg)) ;
00589 for (i = 0 ; i < cdeg_in ; i++)
00590 {
00591 row = Frows [i] ;
00592 ASSERT (row >= 0 && row < n_row) ;
00593 ASSERT (Frpos [row] == i) ;
00594 }
00595 if (UMF_debug > 0 || n_row < 1000)
00596 {
00597 Int cnt = cdeg_in ;
00598 for (row = 0 ; row < n_row ; row++)
00599 {
00600 if (Frpos [row] == EMPTY) cnt++ ;
00601 }
00602 ASSERT (cnt == n_row) ;
00603 }
00604 #endif
00605
00606
00607
00608 ASSERT (pivcol [IN] >= 0 && pivcol [IN] < n_col) ;
00609 ASSERT (NON_PIVOTAL_COL (pivcol [IN])) ;
00610
00611 tpi = Col_tuples [pivcol [IN]] ;
00612 if (tpi)
00613 {
00614 tp = (Tuple *) (Memory + tpi) ;
00615 tp1 = tp ;
00616 tp2 = tp ;
00617 tpend = tp + Col_tlen [pivcol [IN]] ;
00618 for ( ; tp < tpend ; tp++)
00619 {
00620 e = tp->e ;
00621 ASSERT (e > 0 && e <= Work->nel) ;
00622 if (!E [e]) continue ;
00623 f = tp->f ;
00624 p = Memory + E [e] ;
00625 ep = (Element *) p ;
00626 p += UNITS (Element, 1) ;
00627 Cols = (Int *) p ;
00628 if (Cols [f] == EMPTY) continue ;
00629 ASSERT (pivcol [IN] == Cols [f]) ;
00630
00631 Rows = Cols + ep->ncols ;
00632 nrows = ep->nrows ;
00633 p += UNITS (Int, ep->ncols + nrows) ;
00634 C = ((Entry *) p) + f * nrows ;
00635
00636 for (i = 0 ; i < nrows ; i++)
00637 {
00638 row = Rows [i] ;
00639 if (row >= 0)
00640 {
00641 ASSERT (row < n_row) ;
00642 pos = Frpos [row] ;
00643 if (pos < 0)
00644 {
00645
00646 ASSERT (cdeg_in < n_row) ;
00647 if (cdeg_in >= max_cdeg)
00648 {
00649
00650 DEBUGm4 (("cdeg_in failure\n")) ;
00651 return (UMFPACK_ERROR_different_pattern) ;
00652 }
00653 Frpos [row] = cdeg_in ;
00654 Frows [cdeg_in] = row ;
00655 Wy [cdeg_in++] = C [i] ;
00656 }
00657 else
00658 {
00659
00660
00661 ASSERT (pos < max_cdeg) ;
00662 ASSEMBLE (Wy [pos], C [i]) ;
00663 }
00664 }
00665 }
00666 *tp2++ = *tp ;
00667 }
00668 Col_tlen [pivcol [IN]] = tp2 - tp1 ;
00669 }
00670
00671
00672
00673 #ifndef NDEBUG
00674
00675 DEBUG4 (("COL DONE: cdeg "ID"\nREDUCE COL in "ID" max_cdeg "ID"\n",
00676 cdeg_in, pivcol [IN], max_cdeg)) ;
00677 for (i = 0 ; i < cdeg_in ; i++)
00678 {
00679 row = Frows [i] ;
00680 ASSERT (row >= 0 && row < n_row) ;
00681 ASSERT (Frpos [row] == i) ;
00682 }
00683 if (UMF_debug > 0 || n_row < 1000)
00684 {
00685 Int cnt = cdeg_in ;
00686 for (row = 0 ; row < n_row ; row++)
00687 {
00688 if (Frpos [row] == EMPTY) cnt++ ;
00689 }
00690 ASSERT (cnt == n_row) ;
00691 }
00692 #endif
00693
00694 #ifndef NDEBUG
00695 DEBUG4 (("Reduced column: cdeg in "ID" fnrows_max "ID"\n",
00696 cdeg_in, Work->fnrows_max)) ;
00697 for (i = 0 ; i < cdeg_in ; i++)
00698 {
00699 DEBUG4 ((" "ID" "ID" "ID, i, Frows [i], Frpos [Frows [i]])) ;
00700 EDEBUG4 (Wy [i]) ;
00701 DEBUG4 (("\n")) ;
00702 ASSERT (i == Frpos [Frows [i]]) ;
00703 }
00704 ASSERT (cdeg_in <= Work->fnrows_max) ;
00705 #endif
00706
00707
00708
00709
00710
00711 nr_in = cdeg_in - fnrows ;
00712
00713
00714
00715 ASSERT (cdeg_in > 0) ;
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726 status = UMF_row_search (Numeric, Work, Symbolic,
00727 fnrows, cdeg_in, Frows, Frpos,
00728 pivrow [IN], rdeg [IN], Fcols, Wio, nothing, Wy,
00729 pivcol [IN], freebie) ;
00730 ASSERT (!freebie [IN] && !freebie [OUT]) ;
00731
00732
00733
00734
00735
00736 if (status == UMFPACK_ERROR_different_pattern)
00737 {
00738
00739 DEBUGm4 (("row search IN failure\n")) ;
00740 return (UMFPACK_ERROR_different_pattern) ;
00741 }
00742
00743
00744
00745
00746
00747
00748
00749
00750
00751 ASSERT (status != UMFPACK_WARNING_singular_matrix) ;
00752
00753
00754
00755
00756
00757 if (pivrow [IN][IN] != EMPTY)
00758 {
00759
00760
00761
00762
00763
00764
00765 ASSERT (fnrows >= 0 && fncols >= 0) ;
00766
00767 ASSERT (cdeg_in > 0) ;
00768 nc = rdeg [IN][IN] - fncols ;
00769
00770 thiscost =
00771
00772 (nr_in * (fncols - 1)) +
00773
00774 (nc * (cdeg_in - 1)) ;
00775
00776
00777
00778 ASSERT (fnrows + nr_in == cdeg_in) ;
00779 ASSERT (fncols + nc == rdeg [IN][IN]) ;
00780
00781
00782 fnrows_new [IN][IN] = (fnrows-1) + nr_in ;
00783 fncols_new [IN][IN] = (fncols-1) + nc ;
00784
00785
00786 do_extend = TRUE ;
00787
00788 DEBUG2 (("Evaluating option IN-IN:\n")) ;
00789 DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_in "ID" nc "ID"\n",
00790 Work->fnzeros, fnpiv, nr_in, nc)) ;
00791 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
00792
00793
00794
00795 fnzeros = Work->fnzeros + fnpiv * (nr_in + nc) ;
00796
00797 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
00798
00799 new_LUsize = (fnpiv+1) * (fnrows + nr_in + fncols + nc) ;
00800
00801 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
00802
00803
00804
00805
00806
00807 do_update = fnpiv > 0 &&
00808 (((double) fnzeros) / ((double) new_LUsize)) > RELAX2 ;
00809
00810 DEBUG2 (("do_update "ID"\n", do_update))
00811
00812 DEBUG2 (("option IN IN : nr "ID" nc "ID" cost "ID"(0) relax "ID
00813 "\n", nr_in, nc, thiscost, do_extend)) ;
00814
00815
00816 Work->pivot_case = IN_IN ;
00817 bestcost = thiscost ;
00818
00819
00820 Work->do_extend = do_extend ;
00821 Work->do_update = do_update ;
00822 new_fnzeros = fnzeros ;
00823
00824 }
00825
00826
00827
00828
00829
00830 if (pivrow [IN][OUT] != EMPTY)
00831 {
00832
00833
00834
00835
00836 ASSERT (fnrows >= 0 && fncols > 0) ;
00837 ASSERT (cdeg_in > 0) ;
00838
00839
00840
00841 ASSERT (nr_in >= 1) ;
00842
00843
00844 nc = 0 ;
00845 #ifndef NDEBUG
00846 debug_ok = FALSE ;
00847 #endif
00848 for (i = 0 ; i < rdeg [IN][OUT] ; i++)
00849 {
00850 col = Wio [i] ;
00851 DEBUG4 (("counting col "ID" Fcpos[] = "ID"\n", col,
00852 Fcpos [col])) ;
00853 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
00854 if (Fcpos [col] < 0) nc++ ;
00855 #ifndef NDEBUG
00856
00857 if (col == pivcol [IN])
00858 {
00859 ASSERT (Fcpos [col] >= 0) ;
00860 debug_ok = TRUE ;
00861 }
00862 #endif
00863 }
00864 ASSERT (debug_ok) ;
00865
00866 thiscost =
00867
00868 (nc * fnrows) +
00869
00870 ((nr_in-1) * (rdeg [IN][OUT]-1)) ;
00871
00872
00873
00874 extra_cols = ((fncols-1) + nc ) - (rdeg [IN][OUT] - 1) ;
00875 ASSERT (extra_cols >= 0) ;
00876 ASSERT (fncols + nc == extra_cols + rdeg [IN][OUT]) ;
00877 extra_zeros = (nr_in-1) * extra_cols ;
00878
00879 ASSERT (fnrows + nr_in == cdeg_in) ;
00880 ASSERT (fncols + nc == rdeg [IN][OUT] + extra_cols) ;
00881
00882
00883 fnrows_new [IN][OUT] = fnrows + (nr_in-1) ;
00884 fncols_new [IN][OUT] = (fncols-1) + nc ;
00885 relaxed_front = fnrows_new [IN][OUT] * fncols_new [IN][OUT] ;
00886
00887
00888
00889
00890
00891
00892
00893 if (extra_zeros == 0)
00894 {
00895 do_extend = TRUE ;
00896 }
00897 else
00898 {
00899 do_extend = ((double) extra_zeros) <
00900 (relax1 * (double) relaxed_front) ;
00901 }
00902
00903 if (do_extend)
00904 {
00905
00906 thiscost += extra_zeros ;
00907
00908 DEBUG2 (("Evaluating option IN-OUT:\n")) ;
00909 DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_in "ID" nc "ID"\n",
00910 Work->fnzeros, fnpiv, nr_in, nc)) ;
00911 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
00912
00913
00914
00915 fnzeros = Work->fnzeros + fnpiv * (nr_in + nc) ;
00916
00917 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
00918
00919 new_LUsize = (fnpiv+1) * (fnrows + nr_in + fncols + nc) ;
00920
00921 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
00922
00923
00924 do_update = fnpiv > 0 &&
00925 (((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
00926 DEBUG2 (("do_update "ID"\n", do_update))
00927
00928 }
00929 else
00930 {
00931
00932 do_update = fnpiv > 0 ;
00933 fnzeros = 0 ;
00934 DEBUG2 (("IN-OUT do_update forced true: "ID"\n", do_update)) ;
00935
00936
00937
00938 fnrows_new [IN][OUT] = cdeg_in - 1 ;
00939 fncols_new [IN][OUT] = rdeg [IN][OUT] - 1 ;
00940
00941 }
00942
00943 DEBUG2 (("option IN OUT: nr "ID" nc "ID" cost "ID"("ID") relax "ID
00944 "\n", nr_in, nc, thiscost, extra_zeros, do_extend)) ;
00945
00946 if (bestcost == EMPTY || thiscost < bestcost)
00947 {
00948
00949 Work->pivot_case = IN_OUT ;
00950 bestcost = thiscost ;
00951 Work->do_extend = do_extend ;
00952 Work->do_update = do_update ;
00953 new_fnzeros = fnzeros ;
00954 }
00955 }
00956 }
00957
00958
00959
00960
00961
00962 search_pivcol_out = (bestcost != 0 && pivcol [OUT] != EMPTY) ;
00963 if (Symbolic->prefer_diagonal)
00964 {
00965 search_pivcol_out = search_pivcol_out && (pivrow [IN][IN] == EMPTY) ;
00966 }
00967
00968 if (search_pivcol_out)
00969 {
00970
00971 #ifndef NDEBUG
00972 DEBUG2 (("out_col column "ID" NOT in front at position = "ID"\n",
00973 pivcol [OUT], Fcpos [pivcol [OUT]])) ;
00974 UMF_dump_rowcol (1, Numeric, Work, pivcol [OUT], !Symbolic->fixQ) ;
00975 DEBUG2 (("fncols "ID" fncols_max "ID"\n", fncols, Work->fncols_max)) ;
00976 ASSERT (fncols < Work->fncols_max) ;
00977 #endif
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987
00988
00989
00990
00991 cdeg_out = 0 ;
00992
00993 #ifndef NDEBUG
00994
00995 DEBUG4 (("COL ASSEMBLE: cdeg 0\nREDUCE COL out "ID"\n", pivcol [OUT])) ;
00996 if (UMF_debug > 0 || MAX (n_row, n_col) < 1000)
00997 {
00998 for (i = 0 ; i < MAX (n_row, n_col) ; i++)
00999 {
01000 ASSERT (Wp [i] < 0) ;
01001 }
01002 }
01003 DEBUG4 (("max_cdeg: "ID"\n", max_cdeg)) ;
01004 #endif
01005
01006 ASSERT (pivcol [OUT] >= 0 && pivcol [OUT] < n_col) ;
01007 ASSERT (NON_PIVOTAL_COL (pivcol [OUT])) ;
01008
01009 tpi = Col_tuples [pivcol [OUT]] ;
01010 if (tpi)
01011 {
01012 tp = (Tuple *) (Memory + tpi) ;
01013 tp1 = tp ;
01014 tp2 = tp ;
01015 tpend = tp + Col_tlen [pivcol [OUT]] ;
01016 for ( ; tp < tpend ; tp++)
01017 {
01018 e = tp->e ;
01019 ASSERT (e > 0 && e <= Work->nel) ;
01020 if (!E [e]) continue ;
01021 f = tp->f ;
01022 p = Memory + E [e] ;
01023 ep = (Element *) p ;
01024 p += UNITS (Element, 1) ;
01025 Cols = (Int *) p ;
01026 if (Cols [f] == EMPTY) continue ;
01027 ASSERT (pivcol [OUT] == Cols [f]) ;
01028
01029 Rows = Cols + ep->ncols ;
01030 nrows = ep->nrows ;
01031 p += UNITS (Int, ep->ncols + nrows) ;
01032 C = ((Entry *) p) + f * nrows ;
01033
01034 for (i = 0 ; i < nrows ; i++)
01035 {
01036 row = Rows [i] ;
01037 if (row >= 0)
01038 {
01039 ASSERT (row < n_row) ;
01040 pos = Wp [row] ;
01041 if (pos < 0)
01042 {
01043
01044 ASSERT (cdeg_out < n_row) ;
01045 if (cdeg_out >= max_cdeg)
01046 {
01047
01048 DEBUGm4 (("cdeg out failure\n")) ;
01049 return (UMFPACK_ERROR_different_pattern) ;
01050 }
01051 Wp [row] = cdeg_out ;
01052 Wm [cdeg_out] = row ;
01053 Wx [cdeg_out++] = C [i] ;
01054 }
01055 else
01056 {
01057
01058
01059 ASSEMBLE (Wx [pos], C [i]) ;
01060 }
01061 }
01062 }
01063 *tp2++ = *tp ;
01064 }
01065 Col_tlen [pivcol [OUT]] = tp2 - tp1 ;
01066 }
01067
01068
01069
01070 #ifndef NDEBUG
01071 DEBUG4 (("Reduced column: cdeg out "ID"\n", cdeg_out)) ;
01072 for (i = 0 ; i < cdeg_out ; i++)
01073 {
01074 DEBUG4 ((" "ID" "ID" "ID, i, Wm [i], Wp [Wm [i]])) ;
01075 EDEBUG4 (Wx [i]) ;
01076 DEBUG4 (("\n")) ;
01077 ASSERT (i == Wp [Wm [i]]) ;
01078 }
01079 #endif
01080
01081
01082
01083
01084
01085
01086 status = UMF_row_search (Numeric, Work, Symbolic,
01087 0, cdeg_out, Wm, Wp,
01088 pivrow [OUT], rdeg [OUT], Woi, Woo, pivrow [IN], Wx,
01089 pivcol [OUT], freebie) ;
01090
01091
01092
01093
01094
01095 if (status == UMFPACK_ERROR_different_pattern)
01096 {
01097
01098 return (UMFPACK_ERROR_different_pattern) ;
01099 }
01100
01101
01102
01103
01104
01105 for (i = 0 ; i < cdeg_out ; i++)
01106 {
01107 Wp [Wm [i]] = EMPTY ;
01108 }
01109
01110
01111
01112
01113
01114 if (status == UMFPACK_WARNING_singular_matrix)
01115 {
01116
01117
01118
01119
01120
01121
01122
01123 DEBUGm4 (("Note: pivcol [OUT]: "ID" discard\n", pivcol [OUT])) ;
01124 ASSERT ((Work->fnrows == 0 && Work->fncols == 0)) ;
01125 ASSERT (pivcol [IN] == EMPTY) ;
01126
01127
01128 ASSERT (pivcol [OUT] == Work->Candidates [jcand [OUT]]) ;
01129 remove_candidate (jcand [OUT], Work, Symbolic) ;
01130 Work->ndiscard++ ;
01131
01132
01133 DEBUG1 (("Prune tuples of dead outcol: "ID"\n", pivcol [OUT])) ;
01134 Col_tlen [pivcol [OUT]] = 0 ;
01135 UMF_mem_free_tail_block (Numeric, Col_tuples [pivcol [OUT]]) ;
01136 Col_tuples [pivcol [OUT]] = 0 ;
01137
01138
01139 return (UMFPACK_WARNING_singular_matrix) ;
01140 }
01141
01142
01143
01144 if (freebie [IN])
01145 {
01146
01147 Woi = Fcols ;
01148 rdeg [OUT][IN] = rdeg [IN][IN] ;
01149 DEBUG4 (("Freebie in, row "ID"\n", pivrow [IN][IN])) ;
01150 }
01151
01152 if (freebie [OUT])
01153 {
01154
01155 Woo = Wio ;
01156 rdeg [OUT][OUT] = rdeg [IN][OUT] ;
01157 DEBUG4 (("Freebie out, row "ID"\n", pivrow [IN][OUT])) ;
01158 }
01159
01160
01161
01162
01163
01164 if (pivrow [OUT][IN] != EMPTY)
01165 {
01166
01167
01168
01169
01170 ASSERT (fnrows > 0 && fncols >= 0) ;
01171
01172 did_rowmerge = (cdeg_out == 0) ;
01173 if (did_rowmerge)
01174 {
01175
01176
01177 DEBUGm4 (("did row merge OUT col, IN row\n")) ;
01178 Wm [0] = pivrow [OUT][IN] ;
01179 CLEAR (Wx [0]) ;
01180 cdeg_out = 1 ;
01181 ASSERT (nr_out == EMPTY) ;
01182 }
01183
01184 nc = rdeg [OUT][IN] - fncols ;
01185 ASSERT (nc >= 1) ;
01186
01187
01188 nr_out = 0 ;
01189 #ifndef NDEBUG
01190 debug_ok = FALSE ;
01191 #endif
01192 for (i = 0 ; i < cdeg_out ; i++)
01193 {
01194 row = Wm [i] ;
01195 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
01196 if (Frpos [row] < 0 || Frpos [row] >= fnrows) nr_out++ ;
01197 #ifndef NDEBUG
01198
01199 if (row == pivrow [OUT][IN])
01200 {
01201 ASSERT (Frpos [row] >= 0) ;
01202 debug_ok = TRUE ;
01203 }
01204 #endif
01205 }
01206 ASSERT (debug_ok) ;
01207
01208 thiscost =
01209
01210 (nr_out * fncols) +
01211
01212 ((nc-1) * (cdeg_out-1)) ;
01213
01214
01215
01216 extra_rows = ((fnrows-1) + nr_out) - (cdeg_out - 1) ;
01217 ASSERT (extra_rows >= 0) ;
01218 ASSERT (fnrows + nr_out == extra_rows + cdeg_out) ;
01219 extra_zeros = (nc-1) * extra_rows ;
01220
01221 ASSERT (fnrows + nr_out == cdeg_out + extra_rows) ;
01222 ASSERT (fncols + nc == rdeg [OUT][IN]) ;
01223
01224
01225 fnrows_new [OUT][IN] = (fnrows-1) + nr_out ;
01226 fncols_new [OUT][IN] = fncols + (nc-1) ;
01227 relaxed_front = fnrows_new [OUT][IN] * fncols_new [OUT][IN] ;
01228
01229
01230
01231
01232
01233 if (did_rowmerge)
01234 {
01235 do_extend = FALSE ;
01236 }
01237 else
01238 {
01239
01240 if (extra_zeros == 0)
01241 {
01242 do_extend = TRUE ;
01243 }
01244 else
01245 {
01246 do_extend = ((double) extra_zeros) <
01247 (relax1 * (double) relaxed_front) ;
01248 }
01249 }
01250
01251 if (do_extend)
01252 {
01253
01254 thiscost += extra_zeros ;
01255
01256 DEBUG2 (("Evaluating option OUT-IN:\n")) ;
01257 DEBUG2 ((" Work->fnzeros "ID" fnpiv "ID" nr_out "ID" nc "ID"\n",
01258 Work->fnzeros, fnpiv, nr_out, nc)) ;
01259 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
01260
01261
01262
01263 fnzeros = Work->fnzeros + fnpiv * (nr_out + nc) ;
01264
01265 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
01266
01267 new_LUsize = (fnpiv+1) * (fnrows + nr_out + fncols + nc) ;
01268
01269 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
01270
01271
01272 do_update = fnpiv > 0 &&
01273 (((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
01274 DEBUG2 (("do_update "ID"\n", do_update))
01275 }
01276 else
01277 {
01278
01279 do_update = fnpiv > 0 ;
01280 fnzeros = 0 ;
01281 DEBUG2 (("OUT-IN do_update forced true: "ID"\n", do_update)) ;
01282
01283
01284
01285 fnrows_new [OUT][IN] = cdeg_out - 1 ;
01286 fncols_new [OUT][IN] = rdeg [OUT][IN] - 1 ;
01287 }
01288
01289 DEBUG2 (("option OUT IN : nr "ID" nc "ID" cost "ID"("ID") relax "ID
01290 "\n", nr_out, nc, thiscost, extra_zeros, do_extend)) ;
01291
01292 if (bestcost == EMPTY || thiscost < bestcost)
01293 {
01294
01295 Work->pivot_case = OUT_IN ;
01296 bestcost = thiscost ;
01297 Work->do_extend = do_extend ;
01298 Work->do_update = do_update ;
01299 new_fnzeros = fnzeros ;
01300 }
01301 }
01302
01303
01304
01305
01306
01307 if (pivrow [OUT][OUT] != EMPTY)
01308 {
01309
01310
01311
01312 ASSERT (fnrows >= 0 && fncols >= 0) ;
01313
01314 did_rowmerge = (cdeg_out == 0) ;
01315 if (did_rowmerge)
01316 {
01317
01318
01319 DEBUGm4 (("did row merge OUT col, OUT row\n")) ;
01320 Wm [0] = pivrow [OUT][OUT] ;
01321 CLEAR (Wx [0]) ;
01322 cdeg_out = 1 ;
01323 ASSERT (nr_out == EMPTY) ;
01324 nr_out = 1 ;
01325 }
01326
01327 if (fnrows == 0 && fncols == 0)
01328 {
01329
01330 ASSERT (fnpiv == 0) ;
01331 nc = rdeg [OUT][OUT] ;
01332 extra_cols = 0 ;
01333 nr_out = cdeg_out ;
01334 extra_rows = 0 ;
01335 extra_zeros = 0 ;
01336
01337 thiscost = (nc-1) * (cdeg_out-1) ;
01338
01339
01340 fnrows_new [OUT][OUT] = nr_out-1 ;
01341 fncols_new [OUT][OUT] = nc-1 ;
01342 relaxed_front = fnrows_new [OUT][OUT] * fncols_new [OUT][OUT] ;
01343 }
01344 else
01345 {
01346
01347
01348 if (nr_out == EMPTY)
01349 {
01350 nr_out = 0 ;
01351 #ifndef NDEBUG
01352 debug_ok = FALSE ;
01353 #endif
01354 for (i = 0 ; i < cdeg_out ; i++)
01355 {
01356 row = Wm [i] ;
01357 ASSERT (row >= 0 && row < n_row) ;
01358 ASSERT (NON_PIVOTAL_ROW (row)) ;
01359 if (Frpos [row] < 0 || Frpos [row] >= fnrows) nr_out++ ;
01360 #ifndef NDEBUG
01361
01362 if (row == pivrow [OUT][OUT])
01363 {
01364 ASSERT (Frpos [row] < 0 || Frpos [row] >= fnrows) ;
01365 debug_ok = TRUE ;
01366 }
01367 #endif
01368 }
01369 ASSERT (debug_ok) ;
01370 }
01371
01372
01373 nc = 0 ;
01374 #ifndef NDEBUG
01375 debug_ok = FALSE ;
01376 #endif
01377 for (i = 0 ; i < rdeg [OUT][OUT] ; i++)
01378 {
01379 col = Woo [i] ;
01380 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
01381 if (Fcpos [col] < 0) nc++ ;
01382 #ifndef NDEBUG
01383
01384 if (col == pivcol [OUT])
01385 {
01386 ASSERT (Fcpos [col] < 0) ;
01387 debug_ok = TRUE ;
01388 }
01389 #endif
01390 }
01391 ASSERT (debug_ok) ;
01392
01393 extra_cols = (fncols + (nc-1)) - (rdeg [OUT][OUT] - 1) ;
01394 extra_rows = (fnrows + (nr_out-1)) - (cdeg_out - 1) ;
01395 ASSERT (extra_rows >= 0) ;
01396 ASSERT (extra_cols >= 0) ;
01397 extra_zeros = ((nc-1) * extra_rows) + ((nr_out-1) * extra_cols);
01398
01399 ASSERT (fnrows + nr_out == cdeg_out + extra_rows) ;
01400 ASSERT (fncols + nc == rdeg [OUT][OUT] + extra_cols) ;
01401
01402 thiscost =
01403
01404 ((nc-1) * (cdeg_out-1)) +
01405
01406 ((nr_out-1) * (fncols - extra_cols)) ;
01407
01408
01409 fnrows_new [OUT][OUT] = fnrows + (nr_out-1) ;
01410 fncols_new [OUT][OUT] = fncols + (nc-1) ;
01411 relaxed_front = fnrows_new [OUT][OUT] * fncols_new [OUT][OUT] ;
01412
01413 }
01414
01415
01416
01417
01418
01419 if (did_rowmerge)
01420 {
01421 do_extend = FALSE ;
01422 }
01423 else
01424 {
01425
01426 if (extra_zeros == 0)
01427 {
01428 do_extend = TRUE ;
01429 }
01430 else
01431 {
01432 do_extend = ((double) extra_zeros) <
01433 (relax1 * (double) relaxed_front) ;
01434 }
01435 }
01436
01437 if (do_extend)
01438 {
01439
01440 thiscost += extra_zeros ;
01441
01442 DEBUG2 (("Evaluating option OUT-OUT:\n")) ;
01443 DEBUG2 (("Work->fnzeros "ID" fnpiv "ID" nr_out "ID" nc "ID"\n",
01444 Work->fnzeros, fnpiv, nr_out, nc)) ;
01445 DEBUG2 (("fncols "ID" fnrows "ID"\n", fncols, fnrows)) ;
01446
01447
01448
01449 fnzeros = Work->fnzeros + fnpiv * (nr_out + nc) ;
01450
01451 DEBUG2 (("fnzeros "ID"\n", fnzeros)) ;
01452
01453 new_LUsize = (fnpiv+1) * (fnrows + nr_out + fncols + nc) ;
01454
01455 DEBUG2 (("new_LUsize "ID"\n", new_LUsize)) ;
01456
01457
01458 do_update = fnpiv > 0 &&
01459 (((double) fnzeros) / ((double) new_LUsize)) > RELAX3 ;
01460 DEBUG2 (("do_update "ID"\n", do_update))
01461 }
01462 else
01463 {
01464
01465 do_update = fnpiv > 0 ;
01466 fnzeros = 0 ;
01467 DEBUG2 (("OUT-OUT do_update forced true: "ID"\n", do_update)) ;
01468
01469
01470
01471 fnrows_new [OUT][OUT] = cdeg_out - 1 ;
01472 fncols_new [OUT][OUT] = rdeg [OUT][OUT] - 1 ;
01473 }
01474
01475 DEBUG2 (("option OUT OUT: nr "ID" nc "ID" cost "ID"\n",
01476 rdeg [OUT][OUT], cdeg_out, thiscost)) ;
01477
01478 if (bestcost == EMPTY || thiscost < bestcost)
01479 {
01480
01481 Work->pivot_case = OUT_OUT ;
01482 bestcost = thiscost ;
01483 Work->do_extend = do_extend ;
01484 Work->do_update = do_update ;
01485 new_fnzeros = fnzeros ;
01486 }
01487 }
01488 }
01489
01490
01491
01492 ASSERT (Work->pivot_case != EMPTY) ;
01493 DEBUG2 (("local search, best option "ID", best cost "ID"\n",
01494 Work->pivot_case, bestcost)) ;
01495
01496
01497
01498
01499
01500 Work->fnzeros = new_fnzeros ;
01501
01502
01503
01504
01505
01506 switch (Work->pivot_case)
01507 {
01508 case IN_IN:
01509 DEBUG2 (("IN-IN option selected\n")) ;
01510 ASSERT (fnrows > 0 && fncols > 0) ;
01511 Work->pivcol_in_front = TRUE ;
01512 Work->pivrow_in_front = TRUE ;
01513 Work->pivcol = pivcol [IN] ;
01514 Work->pivrow = pivrow [IN][IN] ;
01515 Work->ccdeg = nr_in ;
01516 Work->Wrow = Fcols ;
01517 Work->rrdeg = rdeg [IN][IN] ;
01518 jj = jcand [IN] ;
01519 Work->fnrows_new = fnrows_new [IN][IN] ;
01520 Work->fncols_new = fncols_new [IN][IN] ;
01521 break ;
01522
01523 case IN_OUT:
01524 DEBUG2 (("IN-OUT option selected\n")) ;
01525 ASSERT (fnrows >= 0 && fncols > 0) ;
01526 Work->pivcol_in_front = TRUE ;
01527 Work->pivrow_in_front = FALSE ;
01528 Work->pivcol = pivcol [IN] ;
01529 Work->pivrow = pivrow [IN][OUT] ;
01530 Work->ccdeg = nr_in ;
01531 Work->Wrow = Wio ;
01532 Work->rrdeg = rdeg [IN][OUT] ;
01533 jj = jcand [IN] ;
01534 Work->fnrows_new = fnrows_new [IN][OUT] ;
01535 Work->fncols_new = fncols_new [IN][OUT] ;
01536 break ;
01537
01538 case OUT_IN:
01539 DEBUG2 (("OUT-IN option selected\n")) ;
01540 ASSERT (fnrows > 0 && fncols >= 0) ;
01541 Work->pivcol_in_front = FALSE ;
01542 Work->pivrow_in_front = TRUE ;
01543 Work->pivcol = pivcol [OUT] ;
01544 Work->pivrow = pivrow [OUT][IN] ;
01545 Work->ccdeg = cdeg_out ;
01546
01547 Work->Wrow = Woi ;
01548 Work->rrdeg = rdeg [OUT][IN] ;
01549
01550 jj = jcand [OUT] ;
01551 Work->fnrows_new = fnrows_new [OUT][IN] ;
01552 Work->fncols_new = fncols_new [OUT][IN] ;
01553 break ;
01554
01555 case OUT_OUT:
01556 DEBUG2 (("OUT-OUT option selected\n")) ;
01557 ASSERT (fnrows >= 0 && fncols >= 0) ;
01558 Work->pivcol_in_front = FALSE ;
01559 Work->pivrow_in_front = FALSE ;
01560 Work->pivcol = pivcol [OUT] ;
01561 Work->pivrow = pivrow [OUT][OUT] ;
01562 Work->ccdeg = cdeg_out ;
01563
01564 Work->Wrow = Woo ;
01565 Work->rrdeg = rdeg [OUT][OUT] ;
01566 jj = jcand [OUT] ;
01567 Work->fnrows_new = fnrows_new [OUT][OUT] ;
01568 Work->fncols_new = fncols_new [OUT][OUT] ;
01569 break ;
01570
01571 }
01572
01573 ASSERT (IMPLIES (fnrows == 0 && fncols == 0, Work->pivot_case == OUT_OUT)) ;
01574
01575 if (!Work->pivcol_in_front && pivcol [IN] != EMPTY)
01576 {
01577
01578 for (i = fnrows ; i < cdeg_in ; i++)
01579 {
01580 Frpos [Frows [i]] = EMPTY;
01581 }
01582 }
01583
01584
01585
01586
01587
01588
01589
01590
01591
01592 ASSERT (jj >= 0 && jj < Work->nCandidates) ;
01593 ASSERT (Work->pivcol == Work->Candidates [jj]) ;
01594 remove_candidate (jj, Work, Symbolic) ;
01595
01596
01597
01598
01599
01600 DEBUG1 (("Check frontal growth:\n")) ;
01601 DEBUG1 (("fnrows_new "ID" + 1 = "ID", fnr_curr "ID"\n",
01602 Work->fnrows_new, Work->fnrows_new + 1, fnr_curr)) ;
01603 DEBUG1 (("fncols_new "ID" + 1 = "ID", fnc_curr "ID"\n",
01604 Work->fncols_new, Work->fncols_new + 1, fnc_curr)) ;
01605
01606 Work->do_grow = (Work->fnrows_new + 1 > fnr_curr
01607 || Work->fncols_new + 1 > fnc_curr) ;
01608 if (Work->do_grow)
01609 {
01610 DEBUG0 (("\nNeed to grow frontal matrix, force do_update true\n")) ;
01611
01612
01613
01614 if (!Work->do_update && fnpiv > 0)
01615 {
01616
01617
01618 Work->nforced++ ;
01619 Work->do_update = TRUE ;
01620 }
01621 }
01622
01623
01624
01625
01626
01627
01628
01629
01630
01631
01632
01633
01634
01635
01636
01637
01638
01639
01640
01641
01642
01643
01644
01645
01646
01647
01648
01649
01650
01651
01652
01653
01654
01655
01656
01657
01658
01659
01660
01661
01662 #ifndef NDEBUG
01663 DEBUG2 (("\n\nSEARCH DONE: Pivot col "ID" in: ("ID") pivot row "ID" in: ("ID
01664 ") extend: "ID"\n\n", Work->pivcol, Work->pivcol_in_front,
01665 Work->pivrow, Work->pivrow_in_front, Work->do_extend)) ;
01666 UMF_dump_rowcol (1, Numeric, Work, Work->pivcol, !Symbolic->fixQ) ;
01667 DEBUG2 (("Pivot col "ID": fnrows "ID" ccdeg "ID"\n", Work->pivcol, fnrows,
01668 Work->ccdeg)) ;
01669 if (Work->pivcol_in_front)
01670 {
01671 Int found = FALSE ;
01672 DEBUG3 (("Pivcol in front\n")) ;
01673 for (i = 0 ; i < fnrows ; i++)
01674 {
01675 row = Frows [i] ;
01676 DEBUG3 ((ID": row:: "ID" in front ", i, row)) ;
01677 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
01678 ASSERT (Frpos [row] == i) ;
01679 EDEBUG3 (Wy [i]) ;
01680 if (row == Work->pivrow)
01681 {
01682 DEBUG3 ((" <- pivrow")) ;
01683 found = TRUE ;
01684 }
01685 DEBUG3 (("\n")) ;
01686 }
01687 ASSERT (found == Work->pivrow_in_front) ;
01688 found = FALSE ;
01689 for (i = fnrows ; i < fnrows + Work->ccdeg ; i++)
01690 {
01691 row = Frows [i] ;
01692 DEBUG3 ((ID": row:: "ID" (new)", i, row)) ;
01693 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
01694 ASSERT (Frpos [row] == i) ;
01695 EDEBUG3 (Wy [i]) ;
01696 if (row == Work->pivrow)
01697 {
01698 DEBUG3 ((" <- pivrow")) ;
01699 found = TRUE ;
01700 }
01701 DEBUG3 (("\n")) ;
01702 }
01703 ASSERT (found == !Work->pivrow_in_front) ;
01704 }
01705 else
01706 {
01707 if (Work->do_extend)
01708 {
01709 Int found = FALSE ;
01710 DEBUG3 (("Pivcol not in front (extend)\n")) ;
01711 for (i = 0 ; i < fnrows ; i++)
01712 {
01713 row = Frows [i] ;
01714 DEBUG3 ((ID": row:: "ID" in front ", i, row)) ;
01715 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
01716 ASSERT (Frpos [row] == i) ;
01717 if (row == Work->pivrow)
01718 {
01719 DEBUG3 ((" <- pivrow")) ;
01720 found = TRUE ;
01721 }
01722 DEBUG3 (("\n")) ;
01723 }
01724 ASSERT (found == Work->pivrow_in_front) ;
01725 found = FALSE ;
01726 DEBUG3 (("----\n")) ;
01727 for (i = 0 ; i < Work->ccdeg ; i++)
01728 {
01729 row = Wm [i] ;
01730 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
01731 DEBUG3 ((ID": row:: "ID" ", i, row)) ;
01732 EDEBUG3 (Wx [i]) ;
01733 if (Frpos [row] < 0)
01734 {
01735 DEBUG3 ((" (new) ")) ;
01736 }
01737 if (row == Work->pivrow)
01738 {
01739 DEBUG3 ((" <- pivrow")) ;
01740 found = TRUE ;
01741
01742 if (Work->pivrow_in_front) ASSERT (Frpos [row] >= 0) ;
01743 else ASSERT (Frpos [row] < 0) ;
01744 }
01745 DEBUG3 (("\n")) ;
01746 }
01747 ASSERT (found) ;
01748 }
01749 else
01750 {
01751 Int found = FALSE ;
01752 DEBUG3 (("Pivcol not in front (no extend)\n")) ;
01753 for (i = 0 ; i < Work->ccdeg ; i++)
01754 {
01755 row = Wm [i] ;
01756 ASSERT (row >= 0 && row < n_row && NON_PIVOTAL_ROW (row)) ;
01757 DEBUG3 ((ID": row:: "ID" ", i, row)) ;
01758 EDEBUG3 (Wx [i]) ;
01759 DEBUG3 ((" (new) ")) ;
01760 if (row == Work->pivrow)
01761 {
01762 DEBUG3 ((" <- pivrow")) ;
01763 found = TRUE ;
01764 }
01765 DEBUG3 (("\n")) ;
01766 }
01767 ASSERT (found) ;
01768 }
01769 }
01770 #endif
01771
01772
01773
01774
01775
01776
01777
01778
01779
01780
01781
01782
01783
01784
01785
01786
01787
01788
01789
01790
01791
01792
01793
01794
01795
01796
01797
01798
01799
01800
01801
01802
01803
01804
01805
01806
01807
01808
01809
01810
01811
01812
01813
01814
01815 #ifndef NDEBUG
01816 UMF_dump_rowcol (0, Numeric, Work, Work->pivrow, TRUE) ;
01817 DEBUG2 (("Pivot row "ID":\n", Work->pivrow)) ;
01818 if (Work->pivrow_in_front)
01819 {
01820 Int found = FALSE ;
01821 for (i = 0 ; i < fncols ; i++)
01822 {
01823 col = Fcols [i] ;
01824 DEBUG3 ((" col:: "ID" in front\n", col)) ;
01825 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
01826 ASSERT (Fcpos [col] == i * fnr_curr) ;
01827 if (col == Work->pivcol) found = TRUE ;
01828 }
01829 ASSERT (found == Work->pivcol_in_front) ;
01830 found = FALSE ;
01831 ASSERT (IMPLIES (Work->pivcol_in_front, Fcols == Work->Wrow)) ;
01832 for (i = fncols ; i < Work->rrdeg ; i++)
01833 {
01834 col = Work->Wrow [i] ;
01835 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
01836 ASSERT (Fcpos [col] < 0) ;
01837 if (col == Work->pivcol) found = TRUE ;
01838 else DEBUG3 ((" col:: "ID" (new)\n", col)) ;
01839 }
01840 ASSERT (found == !Work->pivcol_in_front) ;
01841 }
01842 else
01843 {
01844 if (Work->do_extend)
01845 {
01846 Int found = FALSE ;
01847 for (i = 0 ; i < fncols ; i++)
01848 {
01849 col = Fcols [i] ;
01850 DEBUG3 ((" col:: "ID" in front\n", col)) ;
01851 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
01852 ASSERT (Fcpos [col] == i * fnr_curr) ;
01853 if (col == Work->pivcol) found = TRUE ;
01854 }
01855 ASSERT (found == Work->pivcol_in_front) ;
01856 found = FALSE ;
01857 for (i = 0 ; i < Work->rrdeg ; i++)
01858 {
01859 col = Work->Wrow [i] ;
01860 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
01861 if (Fcpos [col] >= 0) continue ;
01862 if (col == Work->pivcol) found = TRUE ;
01863 else DEBUG3 ((" col:: "ID" (new, extend)\n", col)) ;
01864 }
01865 ASSERT (found == !Work->pivcol_in_front) ;
01866 }
01867 else
01868 {
01869 Int found = FALSE ;
01870 for (i = 0 ; i < Work->rrdeg ; i++)
01871 {
01872 col = Work->Wrow [i] ;
01873 ASSERT (col >= 0 && col < n_col && NON_PIVOTAL_COL (col)) ;
01874 if (col == Work->pivcol) found = TRUE ;
01875 else DEBUG3 ((" col:: "ID" (all new)\n", col)) ;
01876 }
01877 ASSERT (found) ;
01878 }
01879 }
01880 #endif
01881
01882
01883
01884
01885
01886 if (Work->do_extend)
01887 {
01888 Work->do_scan2row = (fncols > 0) ;
01889 Work->do_scan2col = (fnrows > 0) ;
01890 }
01891 else
01892 {
01893 Work->do_scan2row = (fncols > 0) && Work->pivrow_in_front ;
01894 Work->do_scan2col = (fnrows > 0) && Work->pivcol_in_front ;
01895 }
01896
01897
01898
01899 DEBUG2 (("LOCAL SEARCH DONE: pivot column "ID" pivot row: "ID"\n",
01900 Work->pivcol, Work->pivrow)) ;
01901 DEBUG2 (("do_extend: "ID"\n", Work->do_extend)) ;
01902 DEBUG2 (("do_update: "ID"\n", Work->do_update)) ;
01903 DEBUG2 (("do_grow: "ID"\n", Work->do_grow)) ;
01904
01905
01906
01907
01908
01909 if (Symbolic->prefer_diagonal
01910 && Work->pivcol < Work->n_col - Symbolic->nempty_col)
01911 {
01912 Diagonal_map = Work->Diagonal_map ;
01913 Diagonal_imap = Work->Diagonal_imap ;
01914 ASSERT (Diagonal_map != (Int *) NULL) ;
01915 ASSERT (Diagonal_imap != (Int *) NULL) ;
01916
01917 row2 = Diagonal_map [Work->pivcol] ;
01918 col2 = Diagonal_imap [Work->pivrow] ;
01919
01920 if (row2 < 0)
01921 {
01922
01923 Work->noff_diagonal++ ;
01924 row2 = UNFLIP (row2) ;
01925 }
01926
01927 ASSERT (Diagonal_imap [row2] == Work->pivcol) ;
01928 ASSERT (UNFLIP (Diagonal_map [col2]) == Work->pivrow) ;
01929
01930 if (row2 != Work->pivrow)
01931 {
01932
01933
01934
01935
01936 DEBUG0 (("Unsymmetric pivot\n")) ;
01937 Diagonal_map [Work->pivcol] = FLIP (Work->pivrow) ;
01938 Diagonal_imap [Work->pivrow] = Work->pivcol ;
01939
01940 Diagonal_map [col2] = FLIP (row2) ;
01941 Diagonal_imap [row2] = col2 ;
01942
01943 }
01944 ASSERT (n_row == n_col) ;
01945 #ifndef NDEBUG
01946 UMF_dump_diagonal_map (Diagonal_map, Diagonal_imap, Symbolic->n1,
01947 Symbolic->n_col, Symbolic->nempty_col) ;
01948 #endif
01949 }
01950
01951 return (UMFPACK_OK) ;
01952 }