/* ========================================================================== */
/* === UMF_analyze ========================================================== */
/* ========================================================================== */

/* -------------------------------------------------------------------------- */
/* UMFPACK Version 4.4, Copyright (c) 2005 by Timothy A. Davis.  CISE Dept,   */
/* Univ. of Florida.  All Rights Reserved.  See ../Doc/License for License.   */
/* web: http://www.cise.ufl.edu/research/sparse/umfpack                       */
/* -------------------------------------------------------------------------- */

/*
    Symbolic LL' factorization of A'*A, to get upper bounds on the size of
    L and U for LU = PAQ, and to determine the frontal matrices and
    (supernodal) column elimination tree.  No fill-reducing column pre-ordering
    is used.

    Returns TRUE if successful, FALSE if out of memory.  UMF_analyze can only
    run out of memory if anzmax (which is Ap [n_row]) is too small.

    Uses workspace of size O(nonzeros in A).  On input, the matrix A is
    stored in row-form at the tail end of Ai.  It is destroyed on output.
    The rows of A must be sorted by increasing first column index.
    The matrix is assumed to be valid.

    Empty rows and columns have already been removed.

*/

#include "umf_internal.h"
#include "umf_apply_order.h"
#include "umf_fsize.h"

/* ========================================================================== */

GLOBAL Int UMF_analyze
(
    Int n_row,		/* A is n_row-by-n_col */
    Int n_col,
    Int Ai [ ],		/* Ai [Ap [0]..Ap[n_row]-1]: column indices */
			/* destroyed on output.  Note that this is NOT the */
			/* user's Ai that was passed to UMFPACK_*symbolic */
			/* size of Ai, Ap [n_row] = anzmax >= anz + n_col */
			/* Ap [0] must be => n_col.  The space to the */
			/* front of Ai is used as workspace. */

    Int Ap [ ],		/* of size MAX (n_row, n_col) + 1 */
			/* Ap [0..n_row]: row pointers */
			/* Row i is in Ai [Ap [i] ... Ap [i+1]-1] */

			/* rows must have smallest col index first, or be */
			/* in sorted form.  Used as workspace of size n_col */
			/* and destroyed. */

			/* Note that this is NOT the */
			/* user's Ap that was passed to UMFPACK_*symbolic */

    Int Up [ ],		/* workspace of size n_col, and output column perm.
			 * for column etree postorder. */

    Int fixQ,

    /* temporary workspaces: */
    Int W [ ],		/* W [0..n_col-1] */
    Int Link [ ],	/* Link [0..n_col-1] */

    /* output: information about each frontal matrix: */
    Int Front_ncols [ ],	/* size n_col */
    Int Front_nrows [ ],	/* of size n_col */
    Int Front_npivcol [ ],	/* of size n_col */
    Int Front_parent [ ],	/* of size n_col */
    Int *nfr_out,

    Int *p_ncompactions		/* number of compactions in UMF_analyze */
)
{
    /* ====================================================================== */
    /* ==== local variables ================================================= */
    /* ====================================================================== */

    Int j, j3, col, k, row, parent, j2, pdest, p, p2, thickness, npivots, nfr,
	i, *Winv, kk, npiv, jnext, krow, knext, pfirst, jlast, ncompactions,
	*Front_stack, *Front_order, *Front_child, *Front_sibling,
	Wflag, npivcol, fallrows, fallcols, fpiv, frows, fcols, *Front_size ;

    nfr = 0 ;
    DEBUG0 (("UMF_analyze: anzmax "ID" anrow "ID" ancol "ID"\n",
	Ap [n_row], n_row, n_col)) ;

    /* ====================================================================== */
    /* ==== initializations ================================================= */
    /* ====================================================================== */

#pragma ivdep
    for (j = 0 ; j < n_col ; j++)
    {
	Link [j] = EMPTY ;
	W [j] = EMPTY ;
	Up [j] = EMPTY ;

	/* Frontal matrix data structure: */
	Front_npivcol [j] = 0 ;		/* number of pivot columns */
	Front_nrows [j] = 0 ;		/* number of rows, incl. pivot rows */
	Front_ncols [j] = 0 ;		/* number of cols, incl. pivot cols */
	Front_parent [j] = EMPTY ;	/* parent front */
	/* Note that only non-pivotal columns are stored in a front (a "row" */
	/* of U) during elimination. */
    }

    /* the rows must be sorted by increasing min col */
    krow = 0 ;
    pfirst = Ap [0] ;
    jlast = EMPTY ;
    jnext = EMPTY ;
    Wflag = 0 ;

    /* this test requires the size of Ai to be >= n_col + nz */
    ASSERT (pfirst >= n_col) ;	/* Ai must be large enough */

    /* pdest points to the first free space in Ai */
    pdest = 0 ;
    ncompactions = 0 ;

    /* ====================================================================== */
    /* === compute symbolic LL' factorization (unsorted) ==================== */
    /* ====================================================================== */

    for (j = 0 ; j < n_col ; j = jnext)
    {
	DEBUG1 (("\n\n============Front "ID" starting. nfr = "ID"\n", j, nfr)) ;

	/* ================================================================== */
	/* === garbage collection =========================================== */
	/* ================================================================== */

	if (pdest + (n_col-j) > pfirst)
	{
	    /* we might run out ... compact the rows of U */

#ifndef NDEBUG
	    DEBUG0 (("UMF_analyze COMPACTION, j="ID" pfirst="ID"\n",
		j, pfirst)) ;
	    for (row = 0 ; row < j ; row++)
	    {
		if (Up [row] != EMPTY)
		{
		    /* this is a live row of U */
		    DEBUG1 (("Live row: "ID" cols: ", row)) ;
		    p = Up [row] ;
		    ASSERT (Front_ncols [row] > Front_npivcol [row]) ;
		    p2 = p + (Front_ncols [row] - Front_npivcol [row]) ;
		    for ( ; p < p2 ; p++)
		    {
			DEBUG1 ((ID, Ai [p])) ;
			ASSERT (p < pfirst) ;
			ASSERT (Ai [p] > row && Ai [p] < n_col) ;
		    }
		    DEBUG1 (("\n")) ;
		}
	    }
	    DEBUG1 (("\nStarting to compact:\n")) ;
#endif

	    pdest = 0 ;
	    ncompactions++ ;
	    for (row = 0 ; row < j ; row++)
	    {
		if (Up [row] != EMPTY)
		{
		    /* this is a live row of U */
		    DEBUG1 (("Live row: "ID" cols: ", row)) ;
		    ASSERT (row < n_col) ;
		    p = Up [row] ;
		    ASSERT (Front_ncols [row] > Front_npivcol [row]) ;
		    p2 = p + (Front_ncols [row] - Front_npivcol [row]) ;
		    Up [row] = pdest ;
		    for ( ; p < p2 ; p++)
		    {
			DEBUG1 ((ID, Ai [p])) ;
			ASSERT (p < pfirst) ;
			ASSERT (Ai [p] > row && Ai [p] < n_col) ;
			Ai [pdest++] = Ai [p] ;
			ASSERT (pdest <= pfirst) ;
		    }
		    DEBUG1 (("\n")) ;
		}
	    }

#ifndef NDEBUG
	    DEBUG1 (("\nAFTER COMPACTION, j="ID" pfirst="ID"\n", j, pfirst)) ;
	    for (row = 0 ; row < j ; row++)
	    {
		if (Up [row] != EMPTY)
		{
		    /* this is a live row of U */
		    DEBUG1 (("Live row: "ID" cols: ", row)) ;
		    p = Up [row] ;
		    ASSERT (Front_ncols [row] > Front_npivcol [row]) ;
		    p2 = p + (Front_ncols [row] - Front_npivcol [row]) ;
		    for ( ; p < p2 ; p++)
		    {
			DEBUG1 ((ID, Ai [p])) ;
			ASSERT (p < pfirst) ;
			ASSERT (Ai [p] > row && Ai [p] < n_col) ;
		    }
		    DEBUG1 (("\n")) ;
		}
	    }
#endif

	}

	if (pdest + (n_col-j) > pfirst)
	{
	    /* :: out of memory in umf_analyze :: */
	    /* it can't happen, if pfirst >= n_col */
	    return (FALSE) ;	/* internal error! */
	}

	/* ------------------------------------------------------------------ */
	/* is the last front a child of this one? */
	/* ------------------------------------------------------------------ */

	if (jlast != EMPTY && Link [j] == jlast)
	{
	    /* yes - create row j by appending to jlast */
	    DEBUG1 (("GOT:last front is child of this one: j "ID" jlast "ID"\n",
		j, jlast)) ;
	    ASSERT (jlast >= 0 && jlast < j) ;

	    Up [j] = Up [jlast] ;
	    Up [jlast] = EMPTY ;

	    /* find the parent, delete column j, and update W */
	    parent = n_col ;
	    for (p = Up [j] ; p < pdest ; )
	    {
		j3 = Ai [p] ;
		DEBUG1 (("Initial row of U: col "ID" ", j3)) ;
		ASSERT (j3 >= 0 && j3 < n_col) ;
		DEBUG1 (("W: "ID" \n", W [j3])) ;
		ASSERT (W [j3] == Wflag) ;
		if (j == j3)
		{
		    DEBUG1 (("Found column j at p = "ID"\n", p)) ;
		    Ai [p] = Ai [--pdest] ;
		}
		else
		{
		    if (j3 < parent)
		    {
			parent = j3 ;
		    }
		    p++ ;
		}
	    }

	    /* delete jlast from the link list of j */
	    Link [j] = Link [jlast] ;

	    ASSERT (Front_nrows [jlast] > Front_npivcol [jlast]) ;
	    thickness = (Front_nrows [jlast] - Front_npivcol [jlast]) ;
	    DEBUG1 (("initial thickness: "ID"\n", thickness)) ;

	}
	else
	{
	    Up [j] = pdest ;
	    parent = n_col ;
	    /* thickness: number of (nonpivotal) rows in frontal matrix j */
	    thickness = 0 ;
	    Wflag = j ;
	}

	/* ================================================================== */
	/* === compute row j of A*A' ======================================== */
	/* ================================================================== */

	/* ------------------------------------------------------------------ */
	/* flag the diagonal entry in row U, but do not add to pattern */
	/* ------------------------------------------------------------------ */

	ASSERT (pdest <= pfirst) ;
	W [j] = Wflag ;

	DEBUG1 (("\nComputing row "ID" of A'*A\n", j)) ;
	DEBUG2 (("	col: "ID" (diagonal)\n", j)) ;

	/* ------------------------------------------------------------------ */
	/* find the rows the contribute to this column j */
	/* ------------------------------------------------------------------ */

	jnext = n_col ;
	for (knext = krow ; knext < n_row ; knext++)
	{
	    ASSERT (Ap [knext] < Ap [knext+1]) ;
	    ASSERT (Ap [knext] >= pfirst && Ap [knext] <= Ap [n_row]) ;
	    jnext = Ai [Ap [knext]] ;
	    ASSERT (jnext >= j) ;
	    if (jnext != j)
	    {
		break ;
	    }
	}

	/* rows krow ... knext-1 all have first column index of j */
	/* (or are empty) */

	/* row knext has first column index of jnext */
	/* if knext = n_row, then jnext is n_col */
	if (knext == n_row)
	{
	    jnext = n_col ;
	}

	ASSERT (jnext > j) ;
	ASSERT (jnext <= n_col) ;

	/* ------------------------------------------------------------------ */
	/* for each nonzero A (k,j) in column j of A do: */
	/* ------------------------------------------------------------------ */

	for (k = krow ; k < knext ; k++)
	{
	    p = Ap [k] ;
	    p2 = Ap [k+1] ;
	    ASSERT (p < p2) ;

	    /* merge row k of A into W */
	    DEBUG2 (("	---- A row "ID" ", k)) ;
	    ASSERT (k >= 0 && k < n_row) ;
	    ASSERT (Ai [p] == j) ;
	    DEBUG2 (("  p "ID" p2 "ID"\n        cols:", p, p2)) ;
	    ASSERT (p  >= pfirst && p  < Ap [n_row]) ;
	    ASSERT (p2 >  pfirst && p2 <= Ap [n_row]) ;
	    for ( ; p < p2 ; p++)
	    {
		/* add to pattern if seen for the first time */
		col = Ai [p] ;
		ASSERT (col >= j && col < n_col) ;
		DEBUG3 ((" "ID, col)) ;
		if (W [col] != Wflag)
		{
		    Ai [pdest++] = col ;
		    ASSERT (pdest <= pfirst) ;
		    /* flag this column has having been seen for row j */
		    W [col] = Wflag ;
		    if (col < parent)
		    {
			parent = col ;
		    }
		}
	    }
	    DEBUG2 (("\n")) ;
	    thickness++ ;
	}

#ifndef NDEBUG
	DEBUG3 (("\nRow "ID" of A'A:\n", j)) ;
	for (p = Up [j] ; p < pdest ; p++)
	{
	    DEBUG3 ((" "ID, Ai [p])) ;
	}
	DEBUG3 (("\n")) ;
#endif

	/* ------------------------------------------------------------------ */
	/* delete rows up to but not including knext */
	/* ------------------------------------------------------------------ */

	krow = knext ;
	pfirst = Ap [knext] ;

	/* we can now use Ai [0..pfirst-1] as workspace for rows of U */

	/* ================================================================== */
	/* === compute jth row of U ========================================= */
	/* ================================================================== */

	/* for each nonzero U (k,j) in column j of U (1:j-1,:) do */
	for (k = Link [j] ; k != EMPTY ; k = Link [k])
	{
	    /* merge row k of U into W */
	    DEBUG2 (("	---- U row "ID, k)) ;
	    ASSERT (k >= 0 && k < n_col) ;
	    ASSERT (Up [k] != EMPTY) ;
	    p = Up [k] ;
	    ASSERT (Front_ncols [k] > Front_npivcol [k]) ;
	    p2 = p + (Front_ncols [k] - Front_npivcol [k]) ;
	    DEBUG2 (("  p "ID" p2 "ID"\n        cols:", p, p2)) ;
	    ASSERT (p <= pfirst) ;
	    ASSERT (p2 <= pfirst) ;
	    for ( ; p < p2 ; p++)
	    {
		/* add to pattern if seen for the first time */
		col = Ai [p] ;
		ASSERT (col >= j && col < n_col) ;
		DEBUG3 ((" "ID, col)) ;
		if (W [col] != Wflag)
		{
		    Ai [pdest++] = col ;
		    ASSERT (pdest <= pfirst) ;
		    /* flag this col has having been seen for row j */
		    W [col] = Wflag ;
		    if (col < parent)
		    {
			parent = col ;
		    }
		}
	    }
	    DEBUG2 (("\n")) ;

	    /* mark the row k as deleted */
	    Up [k] = EMPTY ;

	    ASSERT (Front_nrows [k] > Front_npivcol [k]) ;
	    thickness += (Front_nrows [k] - Front_npivcol [k]) ;
	    ASSERT (Front_parent [k] == j) ;
	}

#ifndef NDEBUG
	DEBUG3 (("\nRow "ID" of U prior to supercolumn detection:\n", j));
	for (p = Up [j] ; p < pdest ; p++)
	{
	    DEBUG3 ((" "ID, Ai [p])) ;
	}
	DEBUG3 (("\n")) ;
	DEBUG1 (("thickness, prior to supercol detect: "ID"\n", thickness)) ;
#endif

	/* ================================================================== */
	/* === quicky mass elimination ====================================== */
	/* ================================================================== */

	/* this code detects some supernodes, but it might miss */
	/* some because the elimination tree (created on the fly) */
	/* is not yet post-ordered, and because the pattern of A'*A */
	/* is also computed on the fly. */

	/* j2 is incremented because the pivot columns are not stored */

	for (j2 = j+1 ; j2 < jnext ; j2++)
	{
	    ASSERT (j2 >= 0 && j2 < n_col) ;
	    if (W [j2] != Wflag || Link [j2] != EMPTY)
	    {
		break ;
	    }
	}

	/* the loop above terminated with j2 at the first non-supernode */
	DEBUG1 (("jnext = "ID"\n", jnext)) ;
	ASSERT (j2 <= jnext) ;
	jnext = j2 ;
	j2-- ;
	DEBUG1 (("j2 = "ID"\n", j2)) ;
	ASSERT (j2 < n_col) ;

	npivots = j2-j+1 ;
	DEBUG1 (("Number of pivot columns: "ID"\n", npivots)) ;

	/* rows j:j2 have the same nonzero pattern, except for columns j:j2-1 */

	if (j2 > j)
	{
	    /* supernode detected, prune the pattern of new row j */
	    ASSERT (parent == j+1) ;
	    ASSERT (j2 < n_col) ;
	    DEBUG1 (("Supernode detected, j "ID" to j2 "ID"\n", j, j2)) ;

	    parent = n_col ;
	    p2 = pdest ;
	    pdest = Up [j] ;
	    for (p = Up [j] ; p < p2 ; p++)
	    {
		col = Ai [p] ;
		ASSERT (col >= 0 && col < n_col) ;
		ASSERT (W [col] == Wflag) ;
		if (col > j2)
		{
		    /* keep this col in the pattern of the new row j */
		    Ai [pdest++] = col ;
		    if (col < parent)
		    {
			parent = col ;
		    }
		}
	    }
	}

	DEBUG1 (("Parent ["ID"] = "ID"\n", j, parent)) ;
	ASSERT (parent > j2) ;

	if (parent == n_col)
	{
	    /* this front has no parent - it is the root of a subtree */
	    parent = EMPTY ;
	}

#ifndef NDEBUG
	DEBUG3 (("\nFinal row "ID" of U after supercolumn detection:\n", j)) ;
	for (p = Up [j] ; p < pdest ; p++)
	{
	    ASSERT (Ai [p] >= 0 && Ai [p] < n_col) ;
	    DEBUG3 ((" "ID" ("ID")", Ai [p], W [Ai [p]])) ;
	    ASSERT (W [Ai [p]] == Wflag) ;
	}
	DEBUG3 (("\n")) ;
#endif

	/* ================================================================== */
	/* === frontal matrix =============================================== */
	/* ================================================================== */

	/* front has Front_npivcol [j] pivot columns */
	/* entire front is Front_nrows [j] -by- Front_ncols [j] */
	/* j is first column in the front */

	npivcol = npivots ;
	fallrows = thickness ;
	fallcols = npivots + pdest - Up [j] ;

	/* number of pivots in the front (rows and columns) */
	fpiv = MIN (npivcol, fallrows) ;

	/* size of contribution block */
	frows = fallrows - fpiv ;
	fcols = fallcols - fpiv ;

	if (frows == 0 || fcols == 0)
	{
	    /* front has no contribution block and thus needs no parent */
	    DEBUG1 (("Frontal matrix evaporation\n")) ;
	    Up [j] = EMPTY ;
	    parent = EMPTY ;
	}

	Front_npivcol [j] = npivots ;
	Front_nrows [j] = fallrows ;
	Front_ncols [j] = fallcols ;
	Front_parent [j] = parent ;
	ASSERT (npivots > 0) ;

	/* Front_parent [j] is the first column of the parent frontal matrix */

	DEBUG1 (("\n\n==== Front "ID", nfr "ID" pivot columns "ID":"ID
	    " all front: "ID"-by-"ID" Parent: "ID"\n", j, nfr, j,j+npivots-1,
	    Front_nrows [j], Front_ncols [j], Front_parent [j])) ;
	nfr++ ;

	/* ================================================================== */
	/* === prepare this row for its parent ============================== */
	/* ================================================================== */

	if (parent != EMPTY)
	{
	    Link [j] = Link [parent] ;
	    Link [parent] = j ;
	}

	ASSERT (jnext > j) ;

	jlast = j ;
    }

    /* ====================================================================== */
    /* === postorder the fronts ============================================= */
    /* ====================================================================== */

    *nfr_out = nfr ;

    Front_order = W ;	/* use W for Front_order [ */

    if (fixQ)
    {
	/* do not postorder the fronts if Q is fixed */
	DEBUG1 (("\nNo postorder (Q is fixed)\n")) ;
	k = 0 ;
	/* Pragma added May 14, 2003.  The Intel compiler icl 6.0 (an old
	 * version) incorrectly vectorizes this loop. */
#pragma novector
	for (j = 0 ; j < n_col ; j++)
	{
	    if (Front_npivcol [j] > 0)
	    {
		Front_order [j] = k++ ;
		DEBUG1 (("Front order of j: "ID" is:"ID"\n", j,
		    Front_order [j])) ;
	    }
	    else
	    {
		Front_order [j] = EMPTY ;
	    }
	}
    }
    else
    {

	/* use Ap for Front_child and use Link for Front_sibling [ */
	Front_child = Ap ;
	Front_sibling = Link ;

	/* use Ai for Front_stack, size of Ai is >= 2*n_col */
	Front_stack = Ai ;
	Front_size = Front_stack + n_col ;

	UMF_fsize (n_col, Front_size, Front_nrows, Front_ncols,
	    Front_parent, Front_npivcol) ;

	AMD_postorder (n_col, Front_parent, Front_npivcol, Front_size,
	    Front_order, Front_child, Front_sibling, Front_stack) ;

	/* done with Front_child, Front_sibling, Front_size, and Front_stack ]*/

	/* ------------------------------------------------------------------ */
	/* construct the column permutation (return in Up) */
	/* ------------------------------------------------------------------ */

	/* Front_order [i] = k means that front i is kth front in the new order.
	 * i is in the range 0 to n_col-1, and k is in the range 0 to nfr-1 */

	/* Use Ai as workspace for Winv [ */
	Winv = Ai ;
	for (k = 0 ; k < nfr ; k++)
	{
	    Winv [k] = EMPTY ;
	}

	/* compute the inverse of Front_order, so that Winv [k] = i */
	/* if Front_order [i] = k */

	DEBUG1 (("\n\nComputing output column permutation:\n")) ;
	for (i = 0 ; i < n_col ; i++)
	{
	    k = Front_order [i] ;
	    if (k != EMPTY)
	    {
		DEBUG1 (("Front "ID" new order: "ID"\n", i, k)) ;
		ASSERT (k >= 0 && k < nfr) ;
		ASSERT (Winv [k] == EMPTY) ;
		Winv [k] = i ;
	    }
	}

	/* Use Up as output permutation */
	kk = 0 ;
	for (k = 0 ; k < nfr ; k++)
	{
	    i = Winv [k] ;
	    DEBUG1 (("Old Front "ID" New Front "ID" npivots "ID" nrows "ID
		" ncols "ID"\n",
		i, k, Front_npivcol [i], Front_nrows [i], Front_ncols [i])) ;
	    ASSERT (i >= 0 && i < n_col) ;
	    ASSERT (Front_npivcol [i] > 0) ;
	    for (npiv = 0 ; npiv < Front_npivcol [i] ; npiv++)
	    {
		Up [kk] = i + npiv ;
		DEBUG1 (("    Cperm ["ID"] = "ID"\n", kk, Up [kk])) ;
		kk++ ;
	    }
	}
	ASSERT (kk == n_col) ;

	/* Winv no longer needed ] */
    }

    /* ---------------------------------------------------------------------- */
    /* apply the postorder traversal to renumber the frontal matrices */
    /* (or pack them in same order, if fixQ) */
    /* ---------------------------------------------------------------------- */

    /* use Ai as workspace */

    UMF_apply_order (Front_npivcol, Front_order, Ai, n_col, nfr) ;
    UMF_apply_order (Front_nrows,   Front_order, Ai, n_col, nfr) ;
    UMF_apply_order (Front_ncols,   Front_order, Ai, n_col, nfr) ;
    UMF_apply_order (Front_parent,  Front_order, Ai, n_col, nfr) ;

    /* fix the parent to refer to the new numbering */
    for (i = 0 ; i < nfr ; i++)
    {
	parent = Front_parent [i] ;
	if (parent != EMPTY)
	{
	    ASSERT (parent >= 0 && parent < n_col) ;
	    ASSERT (Front_order [parent] >= 0 && Front_order [parent] < nfr) ;
	    Front_parent [i] = Front_order [parent] ;
	}
    }

    /* Front_order longer needed ] */

#ifndef NDEBUG
    DEBUG1 (("\nFinal frontal matrices:\n")) ;
    for (i = 0 ; i < nfr ; i++)
    {
	DEBUG1 (("Final front "ID": npiv "ID" nrows "ID" ncols "ID" parent "
	    ID"\n", i, Front_npivcol [i], Front_nrows [i],
	    Front_ncols [i], Front_parent [i])) ;
    }
#endif

    *p_ncompactions = ncompactions ;
    return (TRUE) ;
}