/* ========================================================================== */
/* === klu_analyze ========================================================== */
/* ========================================================================== */
/* Order the matrix using BTF (or not), and then AMD, COLAMD, the natural
* ordering, or the user-provided-function on the blocks. Does not support
* using a given ordering (use klu_analyze_given for that case). */
#include "klu_internal.h"
/* ========================================================================== */
/* === worker =============================================================== */
/* ========================================================================== */
static int worker /* returns KLU_OK or < 0 if error */
(
/* inputs, not modified */
int n, /* A is n-by-n */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
int nblocks, /* # of blocks */
int Pbtf [ ], /* BTF row permutation */
int Qbtf [ ], /* BTF col permutation */
int R [ ], /* size n+1, but only Rbtf [0..nblocks] is used */
int ordering, /* what ordering to use */
/* output only, not defined on input */
int P [ ], /* size n */
int Q [ ], /* size n */
double Lnz [ ], /* size n, but only Lnz [0..nblocks-1] is used */
/* workspace, not defined on input or output */
int Pamd [ ], /* size maxblock */
int Cp [ ], /* size maxblock+1 */
int Ci [ ], /* size nz */
int Ep [ ], /* size maxblock+1 */
int Pinv [ ], /* size maxblock */
/* input/output */
klu_symbolic *Symbolic,
klu_common *Common
)
{
double amd_Info [AMD_INFO], lnz, lnz1, flops ;
int *RowCount, *W, k1, k2, nk, k, block, oldcol, pend, row, newcol,
result, pc, p, newrow, maxnz, nzoff, ok ;
/* ---------------------------------------------------------------------- */
/* initializations */
/* ---------------------------------------------------------------------- */
/* use Pamd as workspace for RowCount and W */
RowCount = Pamd ;
W = Pamd ;
/* compute the inverse of Pbtf */
#ifndef NDEBUG
for (k = 0 ; k < n ; k++)
{
P [k] = EMPTY ;
Q [k] = EMPTY ;
Pinv [k] = EMPTY ;
}
#endif
for (k = 0 ; k < n ; k++)
{
ASSERT (Pbtf [k] >= 0 && Pbtf [k] < n) ;
Pinv [Pbtf [k]] = k ;
}
#ifndef NDEBUG
for (k = 0 ; k < n ; k++) ASSERT (Pinv [k] != EMPTY) ;
#endif
nzoff = 0 ;
lnz = 0 ;
maxnz = 0 ;
flops = 0 ;
Symbolic->symmetry = EMPTY ; /* only computed by AMD */
/* ---------------------------------------------------------------------- */
/* order each block */
/* ---------------------------------------------------------------------- */
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)) ;
if (nk == 1)
{
/* -------------------------------------------------------------- */
/* singleton case */
/* -------------------------------------------------------------- */
Lnz [block] = 1 ;
P [k1] = Pbtf [k1] ;
Q [k1] = Qbtf [k1] ;
oldcol = Q [k1] ;
pend = Ap [oldcol+1] ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
if (Pinv [Ai [p]] < k1)
{
nzoff++ ;
}
else
{
lnz++ ;
}
}
}
else
{
/* -------------------------------------------------------------- */
/* construct the kth block, C */
/* -------------------------------------------------------------- */
Lnz [block] = EMPTY ;
for (row = 0 ; row < nk ; row++)
{
RowCount [row] = 0 ;
}
pc = 0 ;
for (k = k1 ; k < k2 ; k++)
{
newcol = k-k1 ;
Cp [newcol] = pc ;
oldcol = Qbtf [k] ;
pend = Ap [oldcol+1] ;
for (p = Ap [oldcol] ; p < pend ; p++)
{
newrow = Pinv [Ai [p]] ;
if (newrow < k1)
{
nzoff++ ;
}
else
{
/* (newrow,newcol) is an entry in the block */
ASSERT (newrow < k2) ;
newrow -= k1 ;
Ci [pc++] = newrow ;
RowCount [newrow]++ ;
}
}
}
Cp [nk] = pc ;
maxnz = MAX (maxnz, pc) ;
ASSERT (KLU_valid (nk, Cp, Ci, NULL)) ;
/* -------------------------------------------------------------- */
/* order the block C */
/* -------------------------------------------------------------- */
if (ordering == 0)
{
/* ---------------------------------------------------------- */
/* order the block with AMD (C+C') */
/* ---------------------------------------------------------- */
#if 0
/* since AMD requires a sorted input, compute E = C' */
Ep [0] = 0 ;
for (row = 0 ; row < nk ; row++)
{
Ep [row+1] = Ep [row] + RowCount [row] ;
}
ASSERT (Ep [nk] == pc) ;
for (row = 0 ; row < nk ; row++)
{
W [row] = Ep [row] ;
}
for (col = 0 ; col < nk ; col++)
{
pend = Cp [col+1] ;
for (p = Cp [col] ; p < pend ; p++)
{
ASSERT (W [Ci [p]] >= 0 && W [Ci [p]] < pc) ;
Ei [W [Ci [p]]++] = col ;
}
}
ASSERT (KLU_valid (nk, Ep, Ei, NULL)) ;
#endif
/* AMD memory management routines */
amd_malloc = Common->malloc_memory ;
amd_free = Common->free_memory ;
amd_calloc = Common->calloc_memory ;
amd_realloc = Common->realloc_memory ;
result = amd_order (nk, Cp, Ci, Pamd, NULL, amd_Info) ;
if (result == AMD_OUT_OF_MEMORY)
{
return (KLU_OUT_OF_MEMORY) ;
}
else if (result == AMD_INVALID)
{
/* fatal error - something is corrupted in the matrix E */
PRINTF (("AMD invalid!\n")) ;
return (KLU_INVALID) ;
}
/* get the ordering statistics from AMD */
lnz1 = (int) (amd_Info [AMD_LNZ]) + nk ;
Lnz [block] = lnz1 ;
lnz += lnz1 ;
if (pc == maxnz)
{
/* get the symmetry of the biggest block */
Symbolic->symmetry = amd_Info [AMD_SYMMETRY] ;
}
flops += 2 * amd_Info [AMD_NMULTSUBS_LU] + amd_Info [AMD_NDIV] ;
}
else if (ordering == 1)
{
/* ---------------------------------------------------------- */
/* order the block with COLAMD (C) */
/* ---------------------------------------------------------- */
int *AA, cstats [COLAMD_STATS] ;
size_t Alen ;
PRINTF (("calling COLAMD\n")) ;
Alen = colamd_recommended (pc, n, n) ;
if (Alen == 0)
{
/* problem is too large */
return (KLU_TOO_LARGE) ;
}
AA = klu_malloc (Alen, sizeof (int), Common) ;
if (AA == (int *) NULL)
{
/* out of memory */
return (KLU_OUT_OF_MEMORY) ;
}
/* copy the matrix C into AA and Ep */
for (k = 0 ; k <= nk ; k++)
{
Ep [k] = Cp [k] ;
}
for (p = 0 ; p < pc ; p++)
{
AA [p] = Ci [p] ;
}
/* order (and destroy) AA, returning col permutation in Ep */
ok = colamd (nk, nk, (int) Alen, AA, Ep, NULL, cstats) ;
PRINTF (("COLAMD done\n")) ;
/* free the workspace */
AA = klu_free (AA, Common) ;
if (!ok)
{
/* fatal error - something is corrupted in the matrix */
PRINTF (("COLAMD invalid!\n")) ;
return (KLU_INVALID) ;
}
/* copy the permutation from Ep to Pamd */
for (k = 0 ; k < nk ; k++)
{
Pamd [k] = Ep [k] ;
}
/* COLAMD does not return an estimate on fill-in or flops */
flops = EMPTY ;
lnz = EMPTY ;
Lnz [block] = EMPTY ;
PRINTF (("done COLAMD\n")) ;
}
else if (ordering == 3 && Common->user_order != NULL)
{
/* ---------------------------------------------------------- */
/* pass the block to the user-provided ordering function */
/* ---------------------------------------------------------- */
Lnz [block] = (Common->user_order) (nk, Cp, Ci, Pamd,
Common->user_data) ;
if (Lnz [block] == 0)
{
PRINTF (("user ordering function failed\n")) ;
return (KLU_INVALID) ;
}
}
else
{
PRINTF (("invalid ordering\n")) ;
return (KLU_INVALID) ;
}
/* -------------------------------------------------------------- */
/* combine the preordering with the BTF ordering */
/* -------------------------------------------------------------- */
PRINTF (("Pamd, 1-based:\n")) ;
for (k = 0 ; k < nk ; k++)
{
ASSERT (k + k1 < n) ;
ASSERT (Pamd [k] + k1 < n) ;
Q [k + k1] = Qbtf [Pamd [k] + k1] ;
}
for (k = 0 ; k < nk ; k++)
{
ASSERT (k + k1 < n) ;
ASSERT (Pamd [k] + k1 < n) ;
P [k + k1] = Pbtf [Pamd [k] + k1] ;
}
}
}
PRINTF (("nzoff %d Ap[n] %d\n", nzoff, Ap [n])) ;
ASSERT (nzoff >= 0 && nzoff <= Ap [n]) ;
/* return estimates of # of nonzeros in L including diagonal */
Symbolic->lnz = lnz ; /* EMPTY if COLAMD used */
Symbolic->unz = lnz ; /* EMPTY if COLAMD used */
Symbolic->nzoff = nzoff ;
Symbolic->est_flops = flops ; /* EMPTY if COLAMD used */
return (KLU_OK) ;
}
/* ========================================================================== */
/* === order_and_analyze ==================================================== */
/* ========================================================================== */
/* Orders the matrix with or with BTF, then orders each block with AMD, COLAMD,
* or the user ordering function. Does not handle the natural or given
* ordering cases. */
static klu_symbolic *order_and_analyze /* returns NULL if error, or a valid
klu_symbolic object if successful */
(
/* inputs, not modified */
int n, /* A is n-by-n */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
/* --------------------- */
klu_common *Common
)
{
klu_symbolic *Symbolic ;
double *Lnz ;
int *Qbtf, *Cp, *Ci, *Pinv, *Pamd, *Ep, *Pbtf, *P, *Q, *R ;
int nblocks, nz, block, maxblock, k1, k2, nk, j, i, p, pend, do_btf, nzdiag,
ordering, k, ok ;
size_t n1, n5, nz1, m1 ;
/* ---------------------------------------------------------------------- */
/* determine if input matrix is valid, and get # of nonzeros */
/* ---------------------------------------------------------------------- */
/* A is n-by-n, with n > 0. Ap [0] = 0 and nz = Ap [n] >= 0 required.
* Ap [j] <= Ap [j+1] must hold for all j = 0 to n-1. Row indices in Ai
* must be in the range 0 to n-1, and no duplicate entries can be present.
* The list of row indices in each column of A need not be sorted.
*/
if (n <= 0 || Ap == NULL || Ai == NULL)
{
/* Ap and Ai must be present, and n must be > 0 */
Common->status = KLU_INVALID ;
return (NULL) ;
}
nz = Ap [n] ;
if (Ap [0] != 0 || nz < 0)
{
/* nz must be >= 0 and Ap [0] must equal zero */
Common->status = KLU_INVALID ;
return (NULL) ;
}
/* check for size_t overflow */
do_btf = Common->btf ;
do_btf = (do_btf) ? TRUE : FALSE ;
ok = TRUE ;
n1 = klu_add_size_t (n, 1, &ok) ;
nz1 = klu_add_size_t (nz, 1, &ok) ;
n5 = do_btf ? (klu_mult_size_t (n, 5, &ok)) : 0 ;
if (!ok)
{
/* problem is too large */
Common->status = KLU_TOO_LARGE ;
return (NULL) ;
}
for (j = 0 ; j < n ; j++)
{
if (Ap [j] > Ap [j+1])
{
/* column pointers must be non-decreasing */
Common->status = KLU_INVALID ;
return (NULL) ;
}
}
P = klu_malloc (n, sizeof (int), Common) ;
if (Common->status < KLU_OK)
{
return (NULL) ;
}
for (i = 0 ; i < n ; i++)
{
P [i] = EMPTY ;
}
nzdiag = 0 ;
for (j = 0 ; j < n ; j++)
{
pend = Ap [j+1] ;
for (p = Ap [j] ; p < pend ; p++)
{
i = Ai [p] ;
if (i < 0 || i >= n || P [i] == j)
{
/* row index out of range, or duplicate entry */
P = klu_free (P, Common) ;
Common->status = KLU_INVALID ;
return (NULL) ;
}
if (i == j)
{
/* count the number of diagonal entries */
nzdiag++ ;
}
/* flag row i as appearing in column j */
P [i] = j ;
}
}
/* ---------------------------------------------------------------------- */
/* allocate the Symbolic object */
/* ---------------------------------------------------------------------- */
Symbolic = klu_malloc (sizeof (klu_symbolic), 1, Common) ;
if (Common->status < KLU_OK)
{
klu_free (P, Common) ;
return (NULL) ;
}
Q = klu_malloc (n, sizeof (int), Common) ;
R = klu_malloc (n1, sizeof (int), Common) ;
Lnz = klu_malloc (n, sizeof (double), Common) ;
Symbolic->n = n ;
Symbolic->nz = nz ;
Symbolic->P = P ;
Symbolic->Q = Q ;
Symbolic->R = R ;
Symbolic->Lnz = Lnz ;
if (Common->status < KLU_OK)
{
klu_free_symbolic (&Symbolic, Common) ;
return (NULL) ;
}
/* ---------------------------------------------------------------------- */
/* allocate workspace for BTF permutation */
/* ---------------------------------------------------------------------- */
Pbtf = klu_malloc (n, sizeof (int), Common) ;
Qbtf = klu_malloc (n, sizeof (int), Common) ;
if (Common->status < KLU_OK)
{
klu_free (Pbtf, Common) ;
klu_free (Qbtf, Common) ;
klu_free_symbolic (&Symbolic, Common) ;
return ((klu_symbolic *) NULL) ;
}
/* ---------------------------------------------------------------------- */
/* get the common parameters for BTF and ordering method */
/* ---------------------------------------------------------------------- */
ordering = Common->ordering ;
ordering = MAX (ordering, 0) ; /* 0: AMD */
ordering = MIN (ordering, 3) ; /* 1: COLAMD */
Symbolic->ordering = ordering ;
Symbolic->do_btf = do_btf ;
Symbolic->structural_rank = EMPTY ;
/* ---------------------------------------------------------------------- */
/* find the block triangular form (if requested) */
/* ---------------------------------------------------------------------- */
if (do_btf)
{
int *Work ;
Work = klu_malloc (n5, sizeof (int), Common) ;
if (Common->status < KLU_OK)
{
/* out of memory */
klu_free (Pbtf, Common) ;
klu_free (Qbtf, Common) ;
klu_free_symbolic (&Symbolic, Common) ;
return ((klu_symbolic *) NULL) ;
}
nblocks = btf_order (n, Ap, Ai, Pbtf, Qbtf, R,
&(Symbolic->structural_rank), Work) ;
Common->structural_rank = Symbolic->structural_rank ;
klu_free (Work, Common) ;
/* unflip Qbtf if the matrix does not have full structural rank */
if (Symbolic->structural_rank < n)
{
for (k = 0 ; k < n ; k++)
{
Qbtf [k] = MAXTRANS_UNFLIP (Qbtf [k]) ;
}
}
/* find the size of the largest block */
maxblock = 1 ;
for (block = 0 ; block < nblocks ; block++)
{
k1 = R [block] ;
k2 = R [block+1] ;
nk = k2 - k1 ;
PRINTF (("block %d size %d\n", block, nk)) ;
maxblock = MAX (maxblock, nk) ;
}
}
else
{
/* BTF not requested */
nblocks = 1 ;
maxblock = n ;
R [0] = 0 ;
R [1] = n ;
for (k = 0 ; k < n ; k++)
{
Pbtf [k] = k ;
Qbtf [k] = k ;
}
}
Symbolic->nblocks = nblocks ;
PRINTF (("maxblock size %d\n", maxblock)) ;
Symbolic->maxblock = maxblock ;
/* TODO: merge adjacent 1-by-1 blocks into an upper triangular block */
/* ---------------------------------------------------------------------- */
/* allocate more workspace, for klu_analyze_worker */
/* ---------------------------------------------------------------------- */
/* if we get here, n+1 is already safe, so maxblock+1 is safe too */
m1 = ((size_t) maxblock) + 1 ;
Pamd = klu_malloc (maxblock, sizeof (int), Common) ;
Cp = klu_malloc (m1, sizeof (int), Common) ;
Ep = klu_malloc (m1, sizeof (int), Common) ;
Ci = klu_malloc (nz1, sizeof (int), Common) ;
Pinv = klu_malloc (n, sizeof (int), Common) ;
/* ---------------------------------------------------------------------- */
/* order each block of the BTF matrix using AMD or COLAMD */
/* ---------------------------------------------------------------------- */
if (Common->status == KLU_OK)
{
PRINTF (("calling klu_analyze_worker\n")) ;
Common->status = worker (n, Ap, Ai, nblocks, Pbtf, Qbtf, R,
ordering, P, Q, Lnz, Pamd, Cp, Ci, Ep, Pinv, Symbolic, Common) ;
PRINTF (("klu_analyze_worker done\n")) ;
}
/* ---------------------------------------------------------------------- */
/* free all workspace */
/* ---------------------------------------------------------------------- */
klu_free (Pamd, Common) ;
klu_free (Cp, Common) ;
klu_free (Ci, Common) ;
klu_free (Ep, Common) ;
klu_free (Pinv, Common) ;
klu_free (Pbtf, Common) ;
klu_free (Qbtf, Common) ;
/* ---------------------------------------------------------------------- */
/* return the symbolic object */
/* ---------------------------------------------------------------------- */
if (Common->status < KLU_OK)
{
klu_free_symbolic (&Symbolic, Common) ;
}
return (Symbolic) ;
}
/* ========================================================================== */
/* === klu_analyze ====================================================== */
/* ========================================================================== */
klu_symbolic *klu_analyze /* returns NULL if error, or a valid
klu_symbolic object if successful */
(
/* inputs, not modified */
int n, /* A is n-by-n */
int Ap [ ], /* size n+1, column pointers */
int Ai [ ], /* size nz, row indices */
/* -------------------- */
klu_common *Common
)
{
/* ---------------------------------------------------------------------- */
/* get the control parameters for BTF and ordering method */
/* ---------------------------------------------------------------------- */
if (Common == NULL)
{
return (NULL) ;
}
Common->status = KLU_OK ;
Common->structural_rank = EMPTY ;
/* ---------------------------------------------------------------------- */
/* order and analyze */
/* ---------------------------------------------------------------------- */
if (Common->ordering == 2)
{
/* natural ordering */
return (klu_analyze_given (n, Ap, Ai, NULL, NULL, Common)) ;
}
else
{
/* order with P and Q */
return (order_and_analyze (n, Ap, Ai, Common)) ;
}
}
syntax highlighted by Code2HTML, v. 0.9.1