/* ========================================================================== */
/* === 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