/* ========================================================================== */
/* === Cholesky/cholmod_resymbol ============================================ */
/* ========================================================================== */
/* -----------------------------------------------------------------------------
* CHOLMOD/Cholesky Module. Copyright (C) 2005-2006, Timothy A. Davis
* The CHOLMOD/Cholesky Module is licensed under Version 2.1 of the GNU
* Lesser General Public License. See lesser.txt for a text of the license.
* CHOLMOD is also available under other licenses; contact authors for details.
* http://www.cise.ufl.edu/research/sparse
* -------------------------------------------------------------------------- */
/* Recompute the symbolic pattern of L. Entries not in the symbolic pattern
* are dropped. L->Perm can be used (or not) to permute the input matrix A.
*
* These routines are used after a supernodal factorization is converted into
* a simplicial one, to remove zero entries that were added due to relaxed
* supernode amalgamation. They can also be used after a series of downdates
* to remove entries that would no longer be present if the matrix were
* factorized from scratch. A downdate (cholmod_updown) does not remove any
* entries from L.
*
* workspace: Flag (nrow), Head (nrow+1),
* if symmetric: Iwork (2*nrow)
* if unsymmetric: Iwork (2*nrow+ncol).
* Allocates up to 2 copies of its input matrix A (pattern only).
*/
#ifndef NCHOLESKY
#include "cholmod_internal.h"
#include "cholmod_cholesky.h"
/* ========================================================================== */
/* === cholmod_resymbol ===================================================== */
/* ========================================================================== */
/* Remove entries from L that are not in the factorization of P*A*P', P*A*A'*P',
* or P*F*F'*P' (depending on A->stype and whether fset is NULL or not).
*
* cholmod_resymbol is the same as cholmod_resymbol_noperm, except that it
* first permutes A according to L->Perm. A can be upper/lower/unsymmetric,
* in contrast to cholmod_resymbol_noperm (which can be lower or unsym). */
int CHOLMOD(resymbol)
(
/* ---- input ---- */
cholmod_sparse *A, /* matrix to analyze */
Int *fset, /* subset of 0:(A->ncol)-1 */
size_t fsize, /* size of fset */
int pack, /* if TRUE, pack the columns of L */
/* ---- in/out --- */
cholmod_factor *L, /* factorization, entries pruned on output */
/* --------------- */
cholmod_common *Common
)
{
cholmod_sparse *H, *F, *G ;
Int stype, nrow, ncol ;
size_t s ;
int ok = TRUE ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (FALSE) ;
RETURN_IF_NULL (A, FALSE) ;
RETURN_IF_NULL (L, FALSE) ;
RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
Common->status = CHOLMOD_OK ;
if (L->is_super)
{
/* cannot operate on a supernodal factorization */
ERROR (CHOLMOD_INVALID, "cannot operate on supernodal L") ;
return (FALSE) ;
}
if (L->n != A->nrow)
{
/* dimensions must agree */
ERROR (CHOLMOD_INVALID, "A and L dimensions do not match") ;
return (FALSE) ;
}
/* ---------------------------------------------------------------------- */
/* allocate workspace */
/* ---------------------------------------------------------------------- */
stype = A->stype ;
nrow = A->nrow ;
ncol = A->ncol ;
/* s = 2*nrow + (stype ? 0 : ncol) */
s = CHOLMOD(mult_size_t) (nrow, 2, &ok) ;
s = CHOLMOD(add_size_t) (s, (stype ? 0 : ncol), &ok) ;
if (!ok)
{
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
return (FALSE) ;
}
CHOLMOD(allocate_work) (nrow, s, 0, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (FALSE) ;
}
/* ---------------------------------------------------------------------- */
/* permute the input matrix if necessary */
/* ---------------------------------------------------------------------- */
H = NULL ;
G = NULL ;
if (stype > 0)
{
if (L->ordering == CHOLMOD_NATURAL)
{
/* F = triu(A)' */
/* workspace: Iwork (nrow) */
G = CHOLMOD(ptranspose) (A, 0, NULL, NULL, 0, Common) ;
}
else
{
/* F = triu(A(p,p))' */
/* workspace: Iwork (2*nrow) */
G = CHOLMOD(ptranspose) (A, 0, L->Perm, NULL, 0, Common) ;
}
F = G ;
}
else if (stype < 0)
{
if (L->ordering == CHOLMOD_NATURAL)
{
F = A ;
}
else
{
/* G = triu(A(p,p))' */
/* workspace: Iwork (2*nrow) */
G = CHOLMOD(ptranspose) (A, 0, L->Perm, NULL, 0, Common) ;
/* H = G' */
/* workspace: Iwork (nrow) */
H = CHOLMOD(ptranspose) (G, 0, NULL, NULL, 0, Common) ;
F = H ;
}
}
else
{
if (L->ordering == CHOLMOD_NATURAL)
{
F = A ;
}
else
{
/* G = A(p,f)' */
/* workspace: Iwork (nrow if no fset; MAX (nrow,ncol) if fset)*/
G = CHOLMOD(ptranspose) (A, 0, L->Perm, fset, fsize, Common) ;
/* H = G' */
/* workspace: Iwork (ncol) */
H = CHOLMOD(ptranspose) (G, 0, NULL, NULL, 0, Common) ;
F = H ;
}
}
/* No need to check for failure here. cholmod_resymbol_noperm will return
* FALSE if F is NULL. */
/* ---------------------------------------------------------------------- */
/* resymbol */
/* ---------------------------------------------------------------------- */
ok = CHOLMOD(resymbol_noperm) (F, fset, fsize, pack, L, Common) ;
/* ---------------------------------------------------------------------- */
/* free the temporary matrices, if they exist */
/* ---------------------------------------------------------------------- */
CHOLMOD(free_sparse) (&H, Common) ;
CHOLMOD(free_sparse) (&G, Common) ;
return (ok) ;
}
/* ========================================================================== */
/* === cholmod_resymbol_noperm ============================================== */
/* ========================================================================== */
/* Redo symbolic LDL' or LL' factorization of I + F*F' or I+A, where F=A(:,f).
*
* L already exists, but is a superset of the true dynamic pattern (simple
* column downdates and row deletions haven't pruned anything). Just redo the
* symbolic factorization and drop entries that are no longer there. The
* diagonal is not modified. The number of nonzeros in column j of L
* (L->nz[j]) can decrease. The column pointers (L->p[j]) remain unchanged if
* pack is FALSE or if L is not monotonic. Otherwise, the columns of L are
* packed in place.
*
* For the symmetric case, the columns of the lower triangular part of A
* are accessed by column. NOTE that this the transpose of the general case.
*
* For the unsymmetric case, F=A(:,f) is accessed by column.
*
* A need not be sorted, and can be packed or unpacked. If L->Perm is not
* identity, then A must already be permuted according to the permutation used
* to factorize L. The advantage of using this routine is that it does not
* need to create permuted copies of A first.
*
* This routine can be called if L is only partially factored via cholmod_rowfac
* since all it does is prune. If an entry is in F*F' or A, but not in L, it
* isn't added to L.
*
* L must be simplicial LDL' or LL'; it cannot be supernodal or symbolic.
*
* The set f is held in fset and fsize.
* fset = NULL means ":" in MATLAB. fset is ignored.
* fset != NULL means f = fset [0..fset-1].
* fset != NULL and fsize = 0 means f is the empty set.
* There can be no duplicates in fset.
* Common->status is set to CHOLMOD_INVALID if fset is invalid.
*
* workspace: Flag (nrow), Head (nrow+1),
* if symmetric: Iwork (2*nrow)
* if unsymmetric: Iwork (2*nrow+ncol).
* Unlike cholmod_resymbol, this routine does not allocate any temporary
* copies of its input matrix.
*/
int CHOLMOD(resymbol_noperm)
(
/* ---- input ---- */
cholmod_sparse *A, /* matrix to analyze */
Int *fset, /* subset of 0:(A->ncol)-1 */
size_t fsize, /* size of fset */
int pack, /* if TRUE, pack the columns of L */
/* ---- in/out --- */
cholmod_factor *L, /* factorization, entries pruned on output */
/* --------------- */
cholmod_common *Common
)
{
double *Lx, *Lz ;
Int i, j, k, row, parent, p, pend, pdest, ncol, apacked, sorted, nrow, nf,
use_fset, mark, jj, stype, xtype ;
Int *Ap, *Ai, *Anz, *Li, *Lp, *Lnz, *Flag, *Head, *Link, *Anext, *Iwork ;
size_t s ;
int ok = TRUE ;
/* ---------------------------------------------------------------------- */
/* check inputs */
/* ---------------------------------------------------------------------- */
RETURN_IF_NULL_COMMON (FALSE) ;
RETURN_IF_NULL (A, FALSE) ;
RETURN_IF_NULL (L, FALSE) ;
RETURN_IF_XTYPE_INVALID (A, CHOLMOD_PATTERN, CHOLMOD_ZOMPLEX, FALSE) ;
RETURN_IF_XTYPE_INVALID (L, CHOLMOD_REAL, CHOLMOD_ZOMPLEX, FALSE) ;
ncol = A->ncol ;
nrow = A->nrow ;
stype = A->stype ;
ASSERT (IMPLIES (stype != 0, nrow == ncol)) ;
if (stype > 0)
{
/* symmetric, with upper triangular part, not supported */
ERROR (CHOLMOD_INVALID, "symmetric upper not supported ") ;
return (FALSE) ;
}
if (L->is_super)
{
/* cannot operate on a supernodal or symbolic factorization */
ERROR (CHOLMOD_INVALID, "cannot operate on supernodal L") ;
return (FALSE) ;
}
if (L->n != A->nrow)
{
/* dimensions must agree */
ERROR (CHOLMOD_INVALID, "A and L dimensions do not match") ;
return (FALSE) ;
}
Common->status = CHOLMOD_OK ;
/* ---------------------------------------------------------------------- */
/* allocate workspace */
/* ---------------------------------------------------------------------- */
/* s = 2*nrow + (stype ? 0 : ncol) */
s = CHOLMOD(mult_size_t) (nrow, 2, &ok) ;
if (stype != 0)
{
s = CHOLMOD(add_size_t) (s, ncol, &ok) ;
}
if (!ok)
{
ERROR (CHOLMOD_TOO_LARGE, "problem too large") ;
return (FALSE) ;
}
CHOLMOD(allocate_work) (nrow, s, 0, Common) ;
if (Common->status < CHOLMOD_OK)
{
return (FALSE) ; /* out of memory */
}
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
/* ---------------------------------------------------------------------- */
/* get inputs */
/* ---------------------------------------------------------------------- */
Ai = A->i ;
Ap = A->p ;
Anz = A->nz ;
apacked = A->packed ;
sorted = A->sorted ;
Li = L->i ;
Lx = L->x ;
Lz = L->z ;
Lp = L->p ;
Lnz = L->nz ;
xtype = L->xtype ;
/* If L is monotonic on input, then it can be packed or
* unpacked on output, depending on the pack input parameter. */
/* cannot pack a non-monotonic matrix */
if (!(L->is_monotonic))
{
pack = FALSE ;
}
ASSERT (L->nzmax >= (size_t) (Lp [L->n])) ;
pdest = 0 ;
PRINT1 (("\n\n===================== Resymbol pack %d Apacked %d\n",
pack, A->packed)) ;
ASSERT (CHOLMOD(dump_sparse) (A, "ReSymbol A:", Common) >= 0) ;
DEBUG (CHOLMOD(dump_factor) (L, "ReSymbol initial L (i, x):", Common)) ;
/* ---------------------------------------------------------------------- */
/* get workspace */
/* ---------------------------------------------------------------------- */
Flag = Common->Flag ; /* size nrow */
Head = Common->Head ; /* size nrow+1 */
Iwork = Common->Iwork ;
Link = Iwork ; /* size nrow (i/i/l) [ */
Lnz = Iwork + nrow ; /* size nrow (i/i/l), if L not packed */
Anext = Iwork + 2*((size_t) nrow) ; /* size ncol (i/i/l), unsym. only */
for (j = 0 ; j < nrow ; j++)
{
Link [j] = EMPTY ;
}
/* use Lnz in L itself */
Lnz = L->nz ;
ASSERT (Lnz != NULL) ;
/* ---------------------------------------------------------------------- */
/* for the unsymmetric case, queue each column of A (:,f) */
/* ---------------------------------------------------------------------- */
/* place each column of the basis set on the link list corresponding to */
/* the smallest row index in that column */
if (stype == 0)
{
use_fset = (fset != NULL) ;
if (use_fset)
{
nf = fsize ;
/* This is the only O(ncol) loop in cholmod_resymbol.
* It is required only to check the fset. */
for (j = 0 ; j < ncol ; j++)
{
Anext [j] = -2 ;
}
for (jj = 0 ; jj < nf ; jj++)
{
j = fset [jj] ;
if (j < 0 || j > ncol || Anext [j] != -2)
{
/* out-of-range or duplicate entry in fset */
ERROR (CHOLMOD_INVALID, "fset invalid") ;
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
return (FALSE) ;
}
/* flag column j as having been seen */
Anext [j] = EMPTY ;
}
/* the fset is now valid */
ASSERT (CHOLMOD(dump_perm) (fset, nf, ncol, "fset", Common)) ;
}
else
{
nf = ncol ;
}
for (jj = 0 ; jj < nf ; jj++)
{
j = (use_fset) ? (fset [jj]) : jj ;
/* column j is the fset; find the smallest row (if any) */
p = Ap [j] ;
pend = (apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
if (pend > p)
{
k = Ai [p] ;
if (!sorted)
{
for ( ; p < pend ; p++)
{
k = MIN (k, Ai [p]) ;
}
}
/* place column j on link list k */
ASSERT (k >= 0 && k < nrow) ;
Anext [j] = Head [k] ;
Head [k] = j ;
}
}
}
/* ---------------------------------------------------------------------- */
/* recompute symbolic LDL' factorization */
/* ---------------------------------------------------------------------- */
for (k = 0 ; k < nrow ; k++)
{
#ifndef NDEBUG
PRINT1 (("\n\n================== Initial column k = "ID"\n", k)) ;
for (p = Lp [k] ; p < Lp [k] + Lnz [k] ; p++)
{
PRINT1 ((" row: "ID" value: ", Li [p])) ;
PRINT1 (("\n")) ;
}
PRINT1 (("Recomputing LDL, column k = "ID"\n", k)) ;
#endif
/* ------------------------------------------------------------------ */
/* compute column k of I+F*F' or I+A */
/* ------------------------------------------------------------------ */
/* flag the diagonal entry */
mark = CHOLMOD(clear_flag) (Common) ;
Flag [k] = mark ;
PRINT1 ((" row: "ID" (diagonal)\n", k)) ;
if (stype != 0)
{
/* merge column k of A into Flag (lower triangular part only) */
p = Ap [k] ;
pend = (apacked) ? (Ap [k+1]) : (p + Anz [k]) ;
for ( ; p < pend ; p++)
{
i = Ai [p] ;
if (i > k)
{
Flag [i] = mark ;
}
}
}
else
{
/* for each column j whos first row index is in row k */
for (j = Head [k] ; j != EMPTY ; j = Anext [j])
{
/* merge column j of A into Flag */
PRINT1 ((" ---- A column "ID"\n", j)) ;
p = Ap [j] ;
pend = (apacked) ? (Ap [j+1]) : (p + Anz [j]) ;
PRINT1 ((" length "ID" adding\n", pend-p)) ;
for ( ; p < pend ; p++)
{
#ifndef NDEBUG
ASSERT (Ai [p] >= k && Ai [p] < nrow) ;
if (Flag [Ai [p]] < mark) PRINT1 ((" row "ID"\n", Ai [p])) ;
#endif
Flag [Ai [p]] = mark ;
}
}
/* clear the kth link list */
Head [k] = EMPTY ;
}
/* ------------------------------------------------------------------ */
/* compute pruned pattern of kth column of L = union of children */
/* ------------------------------------------------------------------ */
/* for each column j of L whose parent is k */
for (j = Link [k] ; j != EMPTY ; j = Link [j])
{
/* merge column j of L into Flag */
PRINT1 ((" ---- L column "ID"\n", k)) ;
ASSERT (j < k) ;
ASSERT (Lnz [j] > 0) ;
p = Lp [j] ;
pend = p + Lnz [j] ;
ASSERT (Li [p] == j && Li [p+1] == k) ;
p++ ; /* skip past the diagonal entry */
for ( ; p < pend ; p++)
{
/* add to pattern */
ASSERT (Li [p] >= k && Li [p] < nrow) ;
Flag [Li [p]] = mark ;
}
}
/* ------------------------------------------------------------------ */
/* prune the kth column of L */
/* ------------------------------------------------------------------ */
PRINT1 (("Final column of L:\n")) ;
p = Lp [k] ;
pend = p + Lnz [k] ;
if (pack)
{
/* shift column k upwards */
Lp [k] = pdest ;
}
else
{
/* leave column k in place, just reduce Lnz [k] */
pdest = p ;
}
for ( ; p < pend ; p++)
{
ASSERT (pdest < pend) ;
ASSERT (pdest <= p) ;
row = Li [p] ;
ASSERT (row >= k && row < nrow) ;
if (Flag [row] == mark)
{
/* keep this entry */
Li [pdest] = row ;
if (xtype == CHOLMOD_REAL)
{
Lx [pdest] = Lx [p] ;
}
else if (xtype == CHOLMOD_COMPLEX)
{
Lx [2*pdest ] = Lx [2*p ] ;
Lx [2*pdest+1] = Lx [2*p+1] ;
}
else if (xtype == CHOLMOD_ZOMPLEX)
{
Lx [pdest] = Lx [p] ;
Lz [pdest] = Lz [p] ;
}
pdest++ ;
}
}
/* ------------------------------------------------------------------ */
/* prepare this column for its parent */
/* ------------------------------------------------------------------ */
Lnz [k] = pdest - Lp [k] ;
PRINT1 ((" L("ID") length "ID"\n", k, Lnz [k])) ;
ASSERT (Lnz [k] > 0) ;
/* parent is the first entry in the column after the diagonal */
parent = (Lnz [k] > 1) ? (Li [Lp [k] + 1]) : EMPTY ;
PRINT1 (("parent ("ID") = "ID"\n", k, parent)) ;
ASSERT ((parent > k && parent < nrow) || (parent == EMPTY)) ;
if (parent != EMPTY)
{
Link [k] = Link [parent] ;
Link [parent] = k ;
}
}
/* done using Iwork for Link, Lnz (if needed), and Anext ] */
/* ---------------------------------------------------------------------- */
/* convert L to packed, if requested */
/* ---------------------------------------------------------------------- */
if (pack)
{
/* finalize Lp */
Lp [nrow] = pdest ;
/* Shrink L to be just large enough. It cannot fail. */
/* workspace: none */
ASSERT ((size_t) (Lp [nrow]) <= L->nzmax) ;
CHOLMOD(reallocate_factor) (Lp [nrow], L, Common) ;
ASSERT (Common->status >= CHOLMOD_OK) ;
}
/* ---------------------------------------------------------------------- */
/* clear workspace */
/* ---------------------------------------------------------------------- */
CHOLMOD(clear_flag) (Common) ;
DEBUG (CHOLMOD(dump_factor) (L, "ReSymbol final L (i, x):", Common)) ;
ASSERT (CHOLMOD(dump_work) (TRUE, TRUE, 0, Common)) ;
return (TRUE) ;
}
#endif
syntax highlighted by Code2HTML, v. 0.9.1