/* ========================================================================== */
/* === klu_refactor ========================================================= */
/* ========================================================================== */
/* Factor the matrix, after ordering and analyzing it with klu_analyze, and
* factoring it once with klu_factor. This routine cannot do any numerical
* pivoting. The pattern of the input matrix (Ap, Ai) must be identical to
* the pattern given to klu_factor.
*/
#include "klu_internal.h"
/* ========================================================================== */
/* === compute the kth column of L and U ==================================== */
/* ========================================================================== */
#define COMPUTE_KTH_COLUMN_OF_L_AND_U \
{ \
/* compute kth column of U, and update kth column of A */ \
GET_POINTER (LU, Uip, Ulen, Ui, Ux, k, ulen) ; \
for (up = 0 ; up < ulen ; up++) \
{ \
j = Ui [up] ; \
ujk = X [j] ; \
/* X [j] = 0 ;*/ \
CLEAR (X [j]) ; \
Ux [up] = ujk ; \
GET_POINTER (LU, Lip, Llen, Li, Lx, j, llen) ; \
for (p = 0 ; p < llen ; p++) \
{ \
/* X [Li [p]] -= Lx [p] * ujk ;*/ \
MULT_SUB (X [Li [p]], Lx [p], ujk) ; \
} \
} \
/* get the diagonal entry of U */ \
ukk = X [k] ; \
/* X [k] = 0 ;*/ \
CLEAR (X [k]) ; \
/* singular case */ \
if (IS_ZERO (ukk)) \
{ \
/* matrix is numerically singular */ \
Common->status = KLU_SINGULAR ; \
if (Common->numerical_rank == EMPTY) \
{ \
Common->numerical_rank = k+k1 ; \
Common->singular_col = Q [k+k1] ; \
} \
if (Common->halt_if_singular) \
{ \
/* do not continue the factorization */ \
return (FALSE) ; \
} \
} \
Udiag_entry [k] = ukk ; \
/* gather and divide by pivot to get kth column of L */ \
GET_POINTER (LU, Lip, Llen, Li, Lx, k, llen) ; \
for (p = 0 ; p < llen ; p++) \
{ \
i = Li [p] ; \
DIV (Lx [p], X [i], ukk) ; \
CLEAR (X [i]) ; \
} \
}
/* ========================================================================== */
/* === klu_refactor ========================================================= */
/* ========================================================================== */
int KLU_refactor /* returns TRUE if successful, FALSE otherwise */
(
/* inputs, not modified */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
double Ax [ ],
klu_symbolic *Symbolic,
/* input/output */
klu_numeric *Numeric,
klu_common *Common
)
{
Entry ukk, ujk, s ;
Entry *Singleton, *Offx, *Lx, *Ux, *X, *Az, *Udiag_entry ;
double *Rs ;
int *P, *Q, *R, *Pnum, *Offp, *Offi, *Ui, *Li, *Pinv, *Lip, *Uip, *Llen,
*Ulen ;
int **Lbip, **Ubip, **Lblen, **Ublen ;
Unit **LUbx, **Udiag ;
Unit *LU ;
int k1, k2, nk, k, block, oldcol, pend, oldrow, n, p, newrow, scale,
nblocks, poff, i, j, up, ulen, llen, maxblock ;
if (Common == NULL)
{
return (FALSE) ;
}
Common->status = KLU_OK ;
if (Numeric == NULL)
{
/* invalid Numeric object */
Common->status = KLU_INVALID ;
return (FALSE) ;
}
Common->numerical_rank = EMPTY ;
Common->singular_col = EMPTY ;
Az = (Entry *) Ax ;
/* ---------------------------------------------------------------------- */
/* get the contents of the Symbolic object */
/* ---------------------------------------------------------------------- */
n = Symbolic->n ;
P = Symbolic->P ;
Q = Symbolic->Q ;
R = Symbolic->R ;
nblocks = Symbolic->nblocks ;
maxblock = Symbolic->maxblock ;
/* ---------------------------------------------------------------------- */
/* get the contents of the Numeric object */
/* ---------------------------------------------------------------------- */
Pnum = Numeric->Pnum ;
Offp = Numeric->Offp ;
Offi = Numeric->Offi ;
Offx = (Entry *) Numeric->Offx ;
Singleton = (Entry *) Numeric->Singleton ;
Lbip = Numeric->Lbip ;
Lblen = Numeric->Lblen ;
Ubip = Numeric->Ubip ;
Ublen = Numeric->Ublen ;
LUbx = (Unit **) Numeric->LUbx ;
scale = Common->scale ;
if (scale > 0)
{
/* factorization was not scaled, but refactorization is scaled */
if (Numeric->Rs == NULL)
{
Numeric->Rs = klu_malloc (n, sizeof (double), Common) ;
if (Common->status < KLU_OK)
{
Common->status = KLU_OUT_OF_MEMORY ;
return (FALSE) ;
}
}
}
else
{
/* no scaling for refactorization; ensure Numeric->Rs is freed. This
* does nothing if Numeric->Rs is already NULL. */
Numeric->Rs = klu_free (Numeric->Rs, Common) ;
}
Rs = Numeric->Rs ;
Pinv = Numeric->Pinv ;
X = (Entry *) Numeric->Xwork ;
Common->nrealloc = 0 ;
Udiag = (Unit **) Numeric->Udiag ;
/* ---------------------------------------------------------------------- */
/* get control parameters */
/* ---------------------------------------------------------------------- */
PRINTF (("klu_refactor scale %d: n %d nzoff %d nblocks %d maxblock %d\n",
scale, n, Symbolic->nzoff, nblocks, maxblock)) ;
#ifdef TESTING
/* randomly mangle the numerical values to test the refactor code */
{
int block, k1, k2, nk, p ;
for (block = 0 ; block < nblocks ; block++)
{
k1 = R [block] ;
k2 = R [block+1] ;
nk = k2 - k1 ;
if (nk == 1)
{
/* note that this will cause a failure for singular matrices */
Singleton [block] = 100 * block ;
}
else
{
Entry *Ux, *Lx;
Lip = Lbip [block] ;
Llen = Lblen [block] ;
Uip = Ubip [block] ;
Ulen = Ublen [block] ;
LU = LUbx [block] ;
for (p = 0 ; p < nk ; p++)
{
/*Lx = (Entry *) (LU + Lip [p] + UNITS (int, Llen [p]) ) ;*/
GET_X_POINTER (LU, Lip, Llen, Lx, p) ;
for (k = 0 ; k < Llen [p] ; k++)
{
Lx [k] = 42 + p * 1e-5 ; /*check if this is ok*/
}
/*Ux = (Entry *) (LU + Uip [p] + UNITS (int, Ulen [p]) ) ;*/
GET_X_POINTER (LU, Uip, Ulen, Ux, p) ;
for (k = 0 ; k < Ulen [p] ; k++)
{
Ux [k] = 99 + p * 1e-4 ; /*check if this is ok*/
}
}
}
}
for (p = 0 ; p < Symbolic->nzoff ; p++) Offx [p] = p ;
if (Rs != NULL) for (k = 0 ; k < n ; k++) Rs [k] = 1 + k/1000. ;
}
#endif
/* ---------------------------------------------------------------------- */
/* check the input matrix compute the row scale factors, Rs */
/* ---------------------------------------------------------------------- */
/* do no scale, or check the input matrix, if scale < 0 */
if (scale >= 0)
{
/* check for out-of-range indices, but do not check for duplicates */
if (!KLU_scale (scale, n, Ap, Ai, Ax, Rs, NULL, Common))
{
return (FALSE) ;
}
}
/* ---------------------------------------------------------------------- */
/* clear workspace X */
/* ---------------------------------------------------------------------- */
for (k = 0 ; k < maxblock ; k++)
{
/* X [k] = 0 ;*/
CLEAR (X [k]) ;
}
poff = 0 ;
/* ---------------------------------------------------------------------- */
/* factor each block */
/* ---------------------------------------------------------------------- */
if (scale < 0)
{
/* ------------------------------------------------------------------ */
/* no scaling or error-checking */
/* ------------------------------------------------------------------ */
for (block = 0 ; block < nblocks ; block++)
{
/* -------------------------------------------------------------- */
/* the block is from rows/columns k1 to k2-1 */
/* -------------------------------------------------------------- */
k1 = R [block] ;
k2 = R [block+1] ;
nk = k2 - k1 ;
PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1,k2-1,nk)) ;
ASSERT (nk <= maxblock) ;
if (nk == 1)
{
/* ---------------------------------------------------------- */
/* singleton case */
/* ---------------------------------------------------------- */
oldcol = Q [k1] ;
pend = Ap [oldcol+1] ;
CLEAR (s) ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
if (Pinv [Ai [p]] < k1)
{
Offx [poff++] = Az [p] ;
}
else
{
s = Az [p] ;
}
}
Singleton [block] = s ;
}
else
{
/* ---------------------------------------------------------- */
/* construct and factor the kth block */
/* ---------------------------------------------------------- */
Lip = Lbip [block] ;
Llen = Lblen [block] ;
Uip = Ubip [block] ;
Ulen = Ublen [block] ;
LU = LUbx [block] ;
Udiag_entry = (Entry *) Udiag [block] ;
for (k = 0 ; k < nk ; k++)
{
/* scatter kth column of the block into workspace X */
oldcol = Q [k+k1] ;
pend = Ap [oldcol+1] ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
newrow = Pinv [Ai [p]] ;
if (newrow < k1)
{
/* this is an entry in the off-diagonal part */
Offx [poff++] = Az [p] ;
}
else
{
/* (newrow-k1,k) is an entry in the block */
X [newrow-k1] = Az [p] ;
}
}
COMPUTE_KTH_COLUMN_OF_L_AND_U ;
}
}
}
}
else if (scale == 0)
{
/* ------------------------------------------------------------------ */
/* no scaling, but with error-checking */
/* ------------------------------------------------------------------ */
for (block = 0 ; block < nblocks ; block++)
{
/* -------------------------------------------------------------- */
/* the block is from rows/columns k1 to k2-1 */
/* -------------------------------------------------------------- */
k1 = R [block] ;
k2 = R [block+1] ;
nk = k2 - k1 ;
PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1,k2-1,nk)) ;
ASSERT (nk <= maxblock) ;
if (nk == 1)
{
/* ---------------------------------------------------------- */
/* singleton case */
/* ---------------------------------------------------------- */
if (Offp [k1] != poff)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
oldcol = Q [k1] ;
pend = Ap [oldcol+1] ;
CLEAR (s) ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
newrow = Pinv [Ai [p]] ;
if (newrow < k1)
{
if (Offi [poff] != newrow)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
Offx [poff++] = Az [p] ;
}
else
{
if (newrow != k1)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
s = Az [p] ;
}
}
Singleton [block] = s ;
}
else
{
/* ---------------------------------------------------------- */
/* construct and factor the kth block */
/* ---------------------------------------------------------- */
Lip = Lbip [block] ;
Llen = Lblen [block] ;
Uip = Ubip [block] ;
Ulen = Ublen [block] ;
LU = LUbx [block] ;
Udiag_entry = (Entry *) Udiag [block] ;
for (k = 0 ; k < nk ; k++)
{
/* scatter kth column of the block into workspace X */
if (Offp [k+k1] != poff)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
oldcol = Q [k+k1] ;
pend = Ap [oldcol+1] ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
oldrow = Ai [p] ;
newrow = Pinv [oldrow] ;
if (newrow < k1)
{
/* this is an entry in the off-diagonal part */
ASSERT (poff < Symbolic->nzoff) ;
if (Offi [poff] != newrow)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
Offx [poff++] = Az [p] ;
}
else
{
/* (newrow,k) is an entry in the block */
ASSERT (newrow < k2) ;
newrow -= k1 ;
if (newrow >= nk)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
X [newrow] = Az [p] ;
}
}
COMPUTE_KTH_COLUMN_OF_L_AND_U ;
}
}
}
}
else
{
/* ------------------------------------------------------------------ */
/* scaling and error-checking */
/* ------------------------------------------------------------------ */
for (block = 0 ; block < nblocks ; block++)
{
/* -------------------------------------------------------------- */
/* the block is from rows/columns k1 to k2-1 */
/* -------------------------------------------------------------- */
k1 = R [block] ;
k2 = R [block+1] ;
nk = k2 - k1 ;
PRINTF (("BLOCK %d, k1 %d k2-1 %d nk %d\n", block, k1,k2-1,nk)) ;
ASSERT (nk <= maxblock) ;
if (nk == 1)
{
/* ---------------------------------------------------------- */
/* singleton case */
/* ---------------------------------------------------------- */
if (Offp [k1] != poff)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
oldcol = Q [k1] ;
pend = Ap [oldcol+1] ;
CLEAR (s) ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
oldrow = Ai [p] ;
newrow = Pinv [oldrow] ;
if (newrow < k1)
{
if (Offi [poff] != newrow)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
/* Offx [poff++] = Az [p] / Rs [oldrow] ;*/
SCALE_DIV_ASSIGN (Offx [poff++], Az [p], Rs [oldrow]) ;
}
else
{
if (newrow != k1)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
/* s = Az [p] / Rs [oldrow] */
SCALE_DIV_ASSIGN (s, Az [p], Rs [oldrow]) ;
}
}
Singleton [block] = s ;
}
else
{
/* ---------------------------------------------------------- */
/* construct and factor the kth block */
/* ---------------------------------------------------------- */
Lip = Lbip [block] ;
Llen = Lblen [block] ;
Uip = Ubip [block] ;
Ulen = Ublen [block] ;
LU = LUbx [block] ;
Udiag_entry = (Entry *) Udiag [block] ;
for (k = 0 ; k < nk ; k++)
{
/* scatter kth column of the block into workspace X */
if (Offp [k+k1] != poff)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
oldcol = Q [k+k1] ;
pend = Ap [oldcol+1] ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
oldrow = Ai [p] ;
newrow = Pinv [oldrow] ;
if (newrow < k1)
{
/* this is an entry in the off-diagonal part */
ASSERT (poff < Symbolic->nzoff) ;
if (Offi [poff] != newrow)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
SCALE_DIV_ASSIGN (Offx [poff++], Az [p],
Rs [oldrow]) ;
}
else
{
/* (newrow,k) is an entry in the block */
ASSERT (newrow < k2) ;
newrow -= k1 ;
if (newrow >= nk)
{
Common->status = KLU_INVALID ;
return (FALSE) ;
}
SCALE_DIV_ASSIGN (X [newrow], Az [p], Rs [oldrow]) ;
}
}
COMPUTE_KTH_COLUMN_OF_L_AND_U ;
}
}
}
}
/* ---------------------------------------------------------------------- */
/* permute scale factors Rs according to pivotal row order */
/* ---------------------------------------------------------------------- */
if (scale > 0)
{
for (k = 0 ; k < n ; k++)
{
REAL (X [k]) = Rs [Pnum [k]] ;
}
for (k = 0 ; k < n ; k++)
{
Rs [k] = REAL (X [k]) ;
}
}
ASSERT (Offp [n] == poff) ;
ASSERT (Symbolic->nzoff == poff) ;
PRINTF (("\n------------------- Off diagonal entries, new:\n")) ;
ASSERT (KLU_valid (n, Offp, Offi, Offx)) ;
#ifndef NDEBUG
if (Common->status == KLU_OK)
{
int block, k1, k2, nk ;
PRINTF (("\n ########### KLU_BTF_REFACTOR done, nblocks %d\n",nblocks));
for (block = 0 ; block < nblocks ; block++)
{
k1 = R [block] ;
k2 = R [block+1] ;
nk = k2 - k1 ;
PRINTF ((
"\n================klu_refactor output: k1 %d k2 %d nk %d\n",
k1, k2, nk)) ;
if (nk == 1)
{
PRINTF (("singleton ")) ;
PRINT_ENTRY ( Singleton [block] ) ;
}
else
{
Lip = Lbip [block] ;
Llen = Lblen [block] ;
LU = (Unit *) Numeric->LUbx [block] ;
PRINTF (("\n---- L block %d\n", block)) ;
ASSERT (KLU_valid_LU (nk, TRUE, Lip, Llen, LU)) ;
Uip = Ubip [block] ;
Ulen = Ublen [block] ;
PRINTF (("\n---- U block %d\n", block)) ;
ASSERT (KLU_valid_LU (nk, FALSE, Uip, Ulen, LU)) ;
}
}
}
#endif
return (TRUE) ;
}
syntax highlighted by Code2HTML, v. 0.9.1