/* ========================================================================== */
/* === UMF_store_lu ========================================================= */
/* ========================================================================== */
/* -------------------------------------------------------------------------- */
/* UMFPACK Version 5.0, Copyright (c) 1995-2006 by Timothy A. Davis. CISE, */
/* Univ. of Florida. All Rights Reserved. See ../Doc/License for License. */
/* web: http://www.cise.ufl.edu/research/sparse/umfpack */
/* -------------------------------------------------------------------------- */
/*
Store the LU factors. Called by the kernel.
Returns TRUE if successful, FALSE if out of memory.
*/
#include "umf_internal.h"
#include "umf_mem_alloc_head_block.h"
#include "umf_get_memory.h"
/* ========================================================================== */
#ifdef DROP
GLOBAL Int UMF_store_lu_drop
#else
GLOBAL Int UMF_store_lu
#endif
(
NumericType *Numeric,
WorkType *Work
)
{
/* ---------------------------------------------------------------------- */
/* local variables */
/* ---------------------------------------------------------------------- */
Entry pivot_value ;
#ifdef DROP
double droptol ;
#endif
Entry *D, *Lval, *Uval, *Fl1, *Fl2, *Fu1, *Fu2,
*Flublock, *Flblock, *Fublock ;
Int i, k, fnr_curr, fnrows, fncols, row, col, pivrow, pivcol, *Frows,
*Fcols, *Lpattern, *Upattern, *Lpos, *Upos, llen, ulen, fnc_curr, fnpiv,
uilen, lnz, unz, nb, *Lilen,
*Uilen, *Lip, *Uip, *Li, *Ui, pivcol_position, newLchain, newUchain,
pivrow_position, p, size, lip, uip, lnzi, lnzx, unzx, lnz2i, lnz2x,
unz2i, unz2x, zero_pivot, *Pivrow, *Pivcol, kk,
Lnz [MAXNB] ;
#ifndef NDEBUG
Int *Col_degree, *Row_degree ;
#endif
#ifdef DROP
Int all_lnz, all_unz ;
droptol = Numeric->droptol ;
#endif
/* ---------------------------------------------------------------------- */
/* get parameters */
/* ---------------------------------------------------------------------- */
fnrows = Work->fnrows ;
fncols = Work->fncols ;
fnpiv = Work->fnpiv ;
Lpos = Numeric->Lpos ;
Upos = Numeric->Upos ;
Lilen = Numeric->Lilen ;
Uilen = Numeric->Uilen ;
Lip = Numeric->Lip ;
Uip = Numeric->Uip ;
D = Numeric->D ;
Flublock = Work->Flublock ;
Flblock = Work->Flblock ;
Fublock = Work->Fublock ;
fnr_curr = Work->fnr_curr ;
fnc_curr = Work->fnc_curr ;
Frows = Work->Frows ;
Fcols = Work->Fcols ;
#ifndef NDEBUG
Col_degree = Numeric->Cperm ; /* for NON_PIVOTAL_COL macro */
Row_degree = Numeric->Rperm ; /* for NON_PIVOTAL_ROW macro */
#endif
Lpattern = Work->Lpattern ;
llen = Work->llen ;
Upattern = Work->Upattern ;
ulen = Work->ulen ;
nb = Work->nb ;
#ifndef NDEBUG
DEBUG1 (("\nSTORE LU: fnrows "ID
" fncols "ID"\n", fnrows, fncols)) ;
DEBUG2 (("\nFrontal matrix, including all space:\n"
"fnr_curr "ID" fnc_curr "ID" nb "ID"\n"
"fnrows "ID" fncols "ID" fnpiv "ID"\n",
fnr_curr, fnc_curr, nb, fnrows, fncols, fnpiv)) ;
DEBUG2 (("\nJust the active part:\n")) ;
DEBUG7 (("C block: ")) ;
UMF_dump_dense (Work->Fcblock, fnr_curr, fnrows, fncols) ;
DEBUG7 (("L block: ")) ;
UMF_dump_dense (Work->Flblock, fnr_curr, fnrows, fnpiv);
DEBUG7 (("U' block: ")) ;
UMF_dump_dense (Work->Fublock, fnc_curr, fncols, fnpiv) ;
DEBUG7 (("LU block: ")) ;
UMF_dump_dense (Work->Flublock, nb, fnpiv, fnpiv) ;
DEBUG7 (("Current frontal matrix: (prior to store LU)\n")) ;
UMF_dump_current_front (Numeric, Work, TRUE) ;
#endif
Pivrow = Work->Pivrow ;
Pivcol = Work->Pivcol ;
/* ---------------------------------------------------------------------- */
/* store the columns of L */
/* ---------------------------------------------------------------------- */
for (kk = 0 ; kk < fnpiv ; kk++)
{
/* ------------------------------------------------------------------ */
/* one more pivot row and column is being stored into L and U */
/* ------------------------------------------------------------------ */
k = Work->npiv + kk ;
/* ------------------------------------------------------------------ */
/* find the kth pivot row and pivot column */
/* ------------------------------------------------------------------ */
pivrow = Pivrow [kk] ;
pivcol = Pivcol [kk] ;
#ifndef NDEBUG
ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
DEBUGm4 ((
"\n -------------------------------------------------------------"
"Store LU: step " ID"\n", k)) ;
ASSERT (k < MIN (Work->n_row, Work->n_col)) ;
DEBUG2 (("Store column of L, k = "ID", llen "ID"\n", k, llen)) ;
for (i = 0 ; i < llen ; i++)
{
row = Lpattern [i] ;
ASSERT (row >= 0 && row < Work->n_row) ;
DEBUG2 ((" Lpattern["ID"] "ID" Lpos "ID, i, row, Lpos [row])) ;
if (row == pivrow) DEBUG2 ((" <- pivot row")) ;
DEBUG2 (("\n")) ;
ASSERT (i == Lpos [row]) ;
}
#endif
/* ------------------------------------------------------------------ */
/* remove pivot row from L */
/* ------------------------------------------------------------------ */
/* remove pivot row index from current column of L */
/* if a new Lchain starts, then all entries are removed later */
DEBUG2 (("Removing pivrow from Lpattern, k = "ID"\n", k)) ;
ASSERT (!NON_PIVOTAL_ROW (pivrow)) ;
pivrow_position = Lpos [pivrow] ;
if (pivrow_position != EMPTY)
{
/* place the last entry in the column in the */
/* position of the pivot row index */
ASSERT (pivrow == Lpattern [pivrow_position]) ;
row = Lpattern [--llen] ;
/* ASSERT (NON_PIVOTAL_ROW (row)) ; */
Lpattern [pivrow_position] = row ;
Lpos [row] = pivrow_position ;
Lpos [pivrow] = EMPTY ;
}
/* ------------------------------------------------------------------ */
/* store the pivot value, for the diagonal matrix D */
/* ------------------------------------------------------------------ */
/* kk-th column of LU block */
Fl1 = Flublock + kk * nb ;
/* kk-th column of L in the L block */
Fl2 = Flblock + kk * fnr_curr ;
/* kk-th pivot in frontal matrix located in Flublock [kk, kk] */
pivot_value = Fl1 [kk] ;
D [k] = pivot_value ;
zero_pivot = IS_ZERO (pivot_value) ;
DEBUG4 (("Pivot D["ID"]=", k)) ;
EDEBUG4 (pivot_value) ;
DEBUG4 (("\n")) ;
/* ------------------------------------------------------------------ */
/* count nonzeros in kth column of L */
/* ------------------------------------------------------------------ */
lnz = 0 ;
lnz2i = 0 ;
lnz2x = llen ;
#ifdef DROP
all_lnz = 0 ;
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
double s ;
x = Fl1 [i] ;
if (IS_ZERO (x)) continue ;
all_lnz++ ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
lnz++ ;
if (Lpos [Pivrow [i]] == EMPTY) lnz2i++ ;
}
for (i = 0 ; i < fnrows ; i++)
{
Entry x ;
double s ;
x = Fl2 [i] ;
if (IS_ZERO (x)) continue ;
all_lnz++ ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
lnz++ ;
if (Lpos [Frows [i]] == EMPTY) lnz2i++ ;
}
#else
for (i = kk + 1 ; i < fnpiv ; i++)
{
if (IS_ZERO (Fl1 [i])) continue ;
lnz++ ;
if (Lpos [Pivrow [i]] == EMPTY) lnz2i++ ;
}
for (i = 0 ; i < fnrows ; i++)
{
if (IS_ZERO (Fl2 [i])) continue ;
lnz++ ;
if (Lpos [Frows [i]] == EMPTY) lnz2i++ ;
}
#endif
lnz2x += lnz2i ;
/* determine if we start a new Lchain or continue the old one */
if (llen == 0 || zero_pivot)
{
/* llen == 0 means there is no prior Lchain */
/* D [k] == 0 means the pivot column is empty */
newLchain = TRUE ;
}
else
{
newLchain =
/* storage for starting a new Lchain */
UNITS (Entry, lnz) + UNITS (Int, lnz)
<=
/* storage for continuing a prior Lchain */
UNITS (Entry, lnz2x) + UNITS (Int, lnz2i) ;
}
if (newLchain)
{
/* start a new chain for column k of L */
DEBUG2 (("Start new Lchain, k = "ID"\n", k)) ;
pivrow_position = EMPTY ;
/* clear the prior Lpattern */
for (i = 0 ; i < llen ; i++)
{
row = Lpattern [i] ;
Lpos [row] = EMPTY ;
}
llen = 0 ;
lnzi = lnz ;
lnzx = lnz ;
}
else
{
/* continue the prior Lchain */
DEBUG2 (("Continue Lchain, k = "ID"\n", k)) ;
lnzi = lnz2i ;
lnzx = lnz2x ;
}
/* ------------------------------------------------------------------ */
/* allocate space for the column of L */
/* ------------------------------------------------------------------ */
size = UNITS (Int, lnzi) + UNITS (Entry, lnzx) ;
#ifndef NDEBUG
UMF_allocfail = FALSE ;
if (UMF_gprob > 0)
{
double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
UMF_allocfail = rrr < UMF_gprob ;
if (UMF_allocfail) DEBUGm2 (("Random garbage coll. (store LU)\n"));
}
#endif
p = UMF_mem_alloc_head_block (Numeric, size) ;
if (!p)
{
Int r2, c2 ;
/* Do garbage collection, realloc, and try again. */
/* Note that there are pivot rows/columns in current front. */
if (Work->do_grow)
{
/* full compaction of current frontal matrix, since
* UMF_grow_front will be called next anyway. */
r2 = fnrows ;
c2 = fncols ;
}
else
{
/* partial compaction. */
r2 = MAX (fnrows, Work->fnrows_new + 1) ;
c2 = MAX (fncols, Work->fncols_new + 1) ;
}
DEBUGm3 (("get_memory from umf_store_lu:\n")) ;
if (!UMF_get_memory (Numeric, Work, size, r2, c2, TRUE))
{
DEBUGm4 (("out of memory: store LU (1)\n")) ;
return (FALSE) ; /* out of memory */
}
p = UMF_mem_alloc_head_block (Numeric, size) ;
if (!p)
{
DEBUGm4 (("out of memory: store LU (2)\n")) ;
return (FALSE) ; /* out of memory */
}
/* garbage collection may have moved the current front */
fnc_curr = Work->fnc_curr ;
fnr_curr = Work->fnr_curr ;
Flublock = Work->Flublock ;
Flblock = Work->Flblock ;
Fublock = Work->Fublock ;
Fl1 = Flublock + kk * nb ;
Fl2 = Flblock + kk * fnr_curr ;
}
/* ------------------------------------------------------------------ */
/* store the column of L */
/* ------------------------------------------------------------------ */
lip = p ;
Li = (Int *) (Numeric->Memory + p) ;
p += UNITS (Int, lnzi) ;
Lval = (Entry *) (Numeric->Memory + p) ;
p += UNITS (Entry, lnzx) ;
for (i = 0 ; i < lnzx ; i++)
{
CLEAR (Lval [i]) ;
}
/* store the numerical entries */
if (newLchain)
{
/* flag the first column in the Lchain by negating Lip [k] */
lip = -lip ;
ASSERT (llen == 0) ;
#ifdef DROP
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
double s ;
Int row2, pos ;
x = Fl1 [i] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
row2 = Pivrow [i] ;
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
Li [pos] = row2 ;
Lval [pos] = x ;
}
for (i = 0 ; i < fnrows ; i++)
{
Entry x ;
double s ;
Int row2, pos ;
x = Fl2 [i] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
row2 = Frows [i] ;
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
Li [pos] = row2 ;
Lval [pos] = x ;
}
#else
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
Int row2, pos ;
x = Fl1 [i] ;
if (IS_ZERO (x)) continue ;
row2 = Pivrow [i] ;
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
Li [pos] = row2 ;
Lval [pos] = x ;
}
for (i = 0 ; i < fnrows ; i++)
{
Entry x ;
Int row2, pos ;
x = Fl2 [i] ;
if (IS_ZERO (x)) continue ;
row2 = Frows [i] ;
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
Li [pos] = row2 ;
Lval [pos] = x ;
}
#endif
}
else
{
ASSERT (llen > 0) ;
#ifdef DROP
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
double s ;
Int row2, pos ;
x = Fl1 [i] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
row2 = Pivrow [i] ;
pos = Lpos [row2] ;
if (pos == EMPTY)
{
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
*Li++ = row2 ;
}
Lval [pos] = x ;
}
for (i = 0 ; i < fnrows ; i++)
{
Entry x ;
double s ;
Int row2, pos ;
x = Fl2 [i] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
row2 = Frows [i] ;
pos = Lpos [row2] ;
if (pos == EMPTY)
{
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
*Li++ = row2 ;
}
Lval [pos] = x ;
}
#else
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
Int row2, pos ;
x = Fl1 [i] ;
if (IS_ZERO (x)) continue ;
row2 = Pivrow [i] ;
pos = Lpos [row2] ;
if (pos == EMPTY)
{
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
*Li++ = row2 ;
}
Lval [pos] = x ;
}
for (i = 0 ; i < fnrows ; i++)
{
Entry x ;
Int row2, pos ;
x = Fl2 [i] ;
if (IS_ZERO (x)) continue ;
row2 = Frows [i] ;
pos = Lpos [row2] ;
if (pos == EMPTY)
{
pos = llen++ ;
Lpattern [pos] = row2 ;
Lpos [row2] = pos ;
*Li++ = row2 ;
}
Lval [pos] = x ;
}
#endif
}
DEBUG4 (("llen "ID" lnzx "ID"\n", llen, lnzx)) ;
ASSERT (llen == lnzx) ;
ASSERT (lnz <= llen) ;
DEBUG4 (("lnz "ID" \n", lnz)) ;
#ifdef DROP
DEBUG4 (("all_lnz "ID" \n", all_lnz)) ;
ASSERT (lnz <= all_lnz) ;
Numeric->lnz += lnz ;
Numeric->all_lnz += all_lnz ;
Lnz [kk] = all_lnz ;
#else
Numeric->lnz += lnz ;
Numeric->all_lnz += lnz ;
Lnz [kk] = lnz ;
#endif
Numeric->nLentries += lnzx ;
Work->llen = llen ;
Numeric->isize += lnzi ;
/* ------------------------------------------------------------------ */
/* the pivot column is fully assembled and scaled, and is now the */
/* k-th column of L */
/* ------------------------------------------------------------------ */
Lpos [pivrow] = pivrow_position ; /* not aliased */
Lip [pivcol] = lip ; /* aliased with Col_tuples */
Lilen [pivcol] = lnzi ; /* aliased with Col_tlen */
}
/* ---------------------------------------------------------------------- */
/* store the rows of U */
/* ---------------------------------------------------------------------- */
for (kk = 0 ; kk < fnpiv ; kk++)
{
/* ------------------------------------------------------------------ */
/* one more pivot row and column is being stored into L and U */
/* ------------------------------------------------------------------ */
k = Work->npiv + kk ;
/* ------------------------------------------------------------------ */
/* find the kth pivot row and pivot column */
/* ------------------------------------------------------------------ */
pivrow = Pivrow [kk] ;
pivcol = Pivcol [kk] ;
#ifndef NDEBUG
ASSERT (pivrow >= 0 && pivrow < Work->n_row) ;
ASSERT (pivcol >= 0 && pivcol < Work->n_col) ;
DEBUG2 (("Store row of U, k = "ID", ulen "ID"\n", k, ulen)) ;
for (i = 0 ; i < ulen ; i++)
{
col = Upattern [i] ;
DEBUG2 ((" Upattern["ID"] "ID, i, col)) ;
if (col == pivcol) DEBUG2 ((" <- pivot col")) ;
DEBUG2 (("\n")) ;
ASSERT (col >= 0 && col < Work->n_col) ;
ASSERT (i == Upos [col]) ;
}
#endif
/* ------------------------------------------------------------------ */
/* get the pivot value, for the diagonal matrix D */
/* ------------------------------------------------------------------ */
zero_pivot = IS_ZERO (D [k]) ;
/* ------------------------------------------------------------------ */
/* count the nonzeros in the row of U */
/* ------------------------------------------------------------------ */
/* kk-th row of U in the LU block */
Fu1 = Flublock + kk ;
/* kk-th row of U in the U block */
Fu2 = Fublock + kk * fnc_curr ;
unz = 0 ;
unz2i = 0 ;
unz2x = ulen ;
DEBUG2 (("unz2x is "ID", lnzx "ID"\n", unz2x, lnzx)) ;
/* if row k does not end a Uchain, pivcol not included in ulen */
ASSERT (!NON_PIVOTAL_COL (pivcol)) ;
pivcol_position = Upos [pivcol] ;
if (pivcol_position != EMPTY)
{
unz2x-- ;
DEBUG2 (("(exclude pivcol) unz2x is now "ID"\n", unz2x)) ;
}
ASSERT (unz2x >= 0) ;
#ifdef DROP
all_unz = 0 ;
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
double s ;
x = Fu1 [i*nb] ;
if (IS_ZERO (x)) continue ;
all_unz++ ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
unz++ ;
if (Upos [Pivcol [i]] == EMPTY) unz2i++ ;
}
for (i = 0 ; i < fncols ; i++)
{
Entry x ;
double s ;
x = Fu2 [i] ;
if (IS_ZERO (x)) continue ;
all_unz++ ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
unz++ ;
if (Upos [Fcols [i]] == EMPTY) unz2i++ ;
}
#else
for (i = kk + 1 ; i < fnpiv ; i++)
{
if (IS_ZERO (Fu1 [i*nb])) continue ;
unz++ ;
if (Upos [Pivcol [i]] == EMPTY) unz2i++ ;
}
for (i = 0 ; i < fncols ; i++)
{
if (IS_ZERO (Fu2 [i])) continue ;
unz++ ;
if (Upos [Fcols [i]] == EMPTY) unz2i++ ;
}
#endif
unz2x += unz2i ;
ASSERT (IMPLIES (k == 0, ulen == 0)) ;
/* determine if we start a new Uchain or continue the old one */
if (ulen == 0 || zero_pivot)
{
/* ulen == 0 means there is no prior Uchain */
/* D [k] == 0 means the matrix is singular (pivot row might */
/* not be empty, however, but start a new Uchain to prune zero */
/* entries for the deg > 0 test in UMF_u*solve) */
newUchain = TRUE ;
}
else
{
newUchain =
/* approximate storage for starting a new Uchain */
UNITS (Entry, unz) + UNITS (Int, unz)
<=
/* approximate storage for continuing a prior Uchain */
UNITS (Entry, unz2x) + UNITS (Int, unz2i) ;
/* this would be exact, except for the Int to Unit rounding, */
/* because the Upattern is stored only at the end of the Uchain */
}
/* ------------------------------------------------------------------ */
/* allocate space for the row of U */
/* ------------------------------------------------------------------ */
size = 0 ;
if (newUchain)
{
/* store the pattern of the last row in the prior Uchain */
size += UNITS (Int, ulen) ;
unzx = unz ;
}
else
{
unzx = unz2x ;
}
size += UNITS (Entry, unzx) ;
#ifndef NDEBUG
UMF_allocfail = FALSE ;
if (UMF_gprob > 0)
{
double rrr = ((double) (rand ( ))) / (((double) RAND_MAX) + 1) ;
DEBUG4 (("Check random %e %e\n", rrr, UMF_gprob)) ;
UMF_allocfail = rrr < UMF_gprob ;
if (UMF_allocfail) DEBUGm2 (("Random garbage coll. (store LU)\n"));
}
#endif
p = UMF_mem_alloc_head_block (Numeric, size) ;
if (!p)
{
Int r2, c2 ;
/* Do garbage collection, realloc, and try again. */
/* Note that there are pivot rows/columns in current front. */
if (Work->do_grow)
{
/* full compaction of current frontal matrix, since
* UMF_grow_front will be called next anyway. */
r2 = fnrows ;
c2 = fncols ;
}
else
{
/* partial compaction. */
r2 = MAX (fnrows, Work->fnrows_new + 1) ;
c2 = MAX (fncols, Work->fncols_new + 1) ;
}
DEBUGm3 (("get_memory from umf_store_lu:\n")) ;
if (!UMF_get_memory (Numeric, Work, size, r2, c2, TRUE))
{
/* :: get memory, column of L :: */
DEBUGm4 (("out of memory: store LU (1)\n")) ;
return (FALSE) ; /* out of memory */
}
p = UMF_mem_alloc_head_block (Numeric, size) ;
if (!p)
{
/* :: out of memory, column of U :: */
DEBUGm4 (("out of memory: store LU (2)\n")) ;
return (FALSE) ; /* out of memory */
}
/* garbage collection may have moved the current front */
fnc_curr = Work->fnc_curr ;
fnr_curr = Work->fnr_curr ;
Flublock = Work->Flublock ;
Flblock = Work->Flblock ;
Fublock = Work->Fublock ;
Fu1 = Flublock + kk ;
Fu2 = Fublock + kk * fnc_curr ;
}
/* ------------------------------------------------------------------ */
/* store the row of U */
/* ------------------------------------------------------------------ */
uip = p ;
if (newUchain)
{
/* starting a new Uchain - flag this by negating Uip [k] */
uip = -uip ;
DEBUG2 (("Start new Uchain, k = "ID"\n", k)) ;
pivcol_position = EMPTY ;
/* end the prior Uchain */
/* save the current Upattern, and then */
/* clear it and start a new Upattern */
DEBUG2 (("Ending prior chain, k-1 = "ID"\n", k-1)) ;
uilen = ulen ;
Ui = (Int *) (Numeric->Memory + p) ;
Numeric->isize += ulen ;
p += UNITS (Int, ulen) ;
for (i = 0 ; i < ulen ; i++)
{
col = Upattern [i] ;
ASSERT (col >= 0 && col < Work->n_col) ;
Upos [col] = EMPTY ;
Ui [i] = col ;
}
ulen = 0 ;
}
else
{
/* continue the prior Uchain */
DEBUG2 (("Continue Uchain, k = "ID"\n", k)) ;
ASSERT (k > 0) ;
/* remove pivot col index from current row of U */
/* if a new Uchain starts, then all entries are removed later */
DEBUG2 (("Removing pivcol from Upattern, k = "ID"\n", k)) ;
if (pivcol_position != EMPTY)
{
/* place the last entry in the row in the */
/* position of the pivot col index */
ASSERT (pivcol == Upattern [pivcol_position]) ;
col = Upattern [--ulen] ;
ASSERT (col >= 0 && col < Work->n_col) ;
Upattern [pivcol_position] = col ;
Upos [col] = pivcol_position ;
Upos [pivcol] = EMPTY ;
}
/* this row continues the Uchain. Keep track of how much */
/* to trim from the k-th length to get the length of the */
/* (k-1)st row of U */
uilen = unz2i ;
}
Uval = (Entry *) (Numeric->Memory + p) ;
/* p += UNITS (Entry, unzx), no need to increment p */
for (i = 0 ; i < unzx ; i++)
{
CLEAR (Uval [i]) ;
}
if (newUchain)
{
ASSERT (ulen == 0) ;
#ifdef DROP
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
double s ;
Int col2, pos ;
x = Fu1 [i*nb] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
col2 = Pivcol [i] ;
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
Uval [pos] = x ;
}
for (i = 0 ; i < fncols ; i++)
{
Entry x ;
double s ;
Int col2, pos ;
x = Fu2 [i] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
col2 = Fcols [i] ;
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
Uval [pos] = x ;
}
#else
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
Int col2, pos ;
x = Fu1 [i*nb] ;
if (IS_ZERO (x)) continue ;
col2 = Pivcol [i] ;
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
Uval [pos] = x ;
}
for (i = 0 ; i < fncols ; i++)
{
Entry x ;
Int col2, pos ;
x = Fu2 [i] ;
if (IS_ZERO (x)) continue ;
col2 = Fcols [i] ;
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
Uval [pos] = x ;
}
#endif
}
else
{
ASSERT (ulen > 0) ;
/* store the numerical entries and find new nonzeros */
#ifdef DROP
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
double s ;
Int col2, pos ;
x = Fu1 [i*nb] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
col2 = Pivcol [i] ;
pos = Upos [col2] ;
if (pos == EMPTY)
{
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
}
Uval [pos] = x ;
}
for (i = 0 ; i < fncols ; i++)
{
Entry x ;
double s ;
Int col2, pos ;
x = Fu2 [i] ;
APPROX_ABS (s, x) ;
if (s <= droptol) continue ;
col2 = Fcols [i] ;
pos = Upos [col2] ;
if (pos == EMPTY)
{
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
}
Uval [pos] = x ;
}
#else
for (i = kk + 1 ; i < fnpiv ; i++)
{
Entry x ;
Int col2, pos ;
x = Fu1 [i*nb] ;
if (IS_ZERO (x)) continue ;
col2 = Pivcol [i] ;
pos = Upos [col2] ;
if (pos == EMPTY)
{
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
}
Uval [pos] = x ;
}
for (i = 0 ; i < fncols ; i++)
{
Entry x ;
Int col2, pos ;
x = Fu2 [i] ;
if (IS_ZERO (x)) continue ;
col2 = Fcols [i] ;
pos = Upos [col2] ;
if (pos == EMPTY)
{
pos = ulen++ ;
Upattern [pos] = col2 ;
Upos [col2] = pos ;
}
Uval [pos] = x ;
}
#endif
}
ASSERT (ulen == unzx) ;
ASSERT (unz <= ulen) ;
DEBUG4 (("unz "ID" \n", unz)) ;
#ifdef DROP
DEBUG4 (("all_unz "ID" \n", all_unz)) ;
ASSERT (unz <= all_unz) ;
Numeric->unz += unz ;
Numeric->all_unz += all_unz ;
/* count the "true" flops, based on LU pattern only */
Numeric->flops += DIV_FLOPS * Lnz [kk] /* scale pivot column */
+ MULTSUB_FLOPS * (Lnz [kk] * all_unz) ; /* outer product */
#else
Numeric->unz += unz ;
Numeric->all_unz += unz ;
/* count the "true" flops, based on LU pattern only */
Numeric->flops += DIV_FLOPS * Lnz [kk] /* scale pivot column */
+ MULTSUB_FLOPS * (Lnz [kk] * unz) ; /* outer product */
#endif
Numeric->nUentries += unzx ;
Work->ulen = ulen ;
DEBUG1 (("Work->ulen = "ID" at end of pivot step, k: "ID"\n", ulen, k));
/* ------------------------------------------------------------------ */
/* the pivot row is the k-th row of U */
/* ------------------------------------------------------------------ */
Upos [pivcol] = pivcol_position ; /* not aliased */
Uip [pivrow] = uip ; /* aliased with Row_tuples */
Uilen [pivrow] = uilen ; /* aliased with Row_tlen */
}
/* ---------------------------------------------------------------------- */
/* no more pivots in frontal working array */
/* ---------------------------------------------------------------------- */
Work->npiv += fnpiv ;
Work->fnpiv = 0 ;
Work->fnzeros = 0 ;
return (TRUE) ;
}
syntax highlighted by Code2HTML, v. 0.9.1