/* ========================================================================== */
/* === CHOLMOD/MATLAB/symbfact2 mexFunction ================================= */
/* ========================================================================== */

/* -----------------------------------------------------------------------------
 * CHOLMOD/MATLAB Module.  Copyright (C) 2005-2006, Timothy A. Davis
 * The CHOLMOD/MATLAB Module is licensed under Version 2.0 of the GNU
 * General Public License.  See gpl.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
 * MATLAB(tm) is a Trademark of The MathWorks, Inc.
 * -------------------------------------------------------------------------- */

/* Usage:
 *
 *  count = symbfact2 (A)		returns row counts of R=chol(A)
 *  count = symbfact2 (A,'col')		returns row counts of R=chol(A'*A)
 *  count = symbfact2 (A,'sym')		same as symbfact2(A)
 *  count = symbfact2 (A,'lo')		same as symbfact2(A'), uses tril(A)
 *  count = symbfact2 (A,'row')		returns row counts of R=chol(A*A')
 *
 *	[count, h, parent, post, R] = symbfact2 (...)
 *
 *	h: height of the elimination tree
 *	parent: the elimination tree itself
 *	post: postordering of the elimination tree
 *	R: a 0-1 matrix whose structure is that of chol(A) or chol(A'*A)
 *		for the 'col' case
 *
 *	symbfact2(A) and symbfact2(A,'sym') uses triu(A).
 *	symbfact2(A,'lo') uses tril(A).
 *
 * These forms return L = R' instead of R.  They are faster and take less
 * memory.  They return the same count, h, parent, and post outputs.
 *
 *  [count, h, parent, post, L] = symbfact2 (A,'col','L')
 *  [count, h, parent, post, L] = symbfact2 (A,'sym','L')
 *  [count, h, parent, post, L] = symbfact2 (A,'lo', 'L')
 *  [count, h, parent, post, L] = symbfact2 (A,'row','L')
 */

#include "cholmod_matlab.h"

void mexFunction
(
    int nargout,
    mxArray *pargout [ ],
    int nargin,
    const mxArray *pargin [ ]
)
{
    double dummy = 0 ;
    double *Lx ;
    int *Parent, *Post, *ColCount, *First, *Level, *Rp, *Ri, *Lp, *Li, *W ;
    cholmod_sparse *A, Amatrix, *F, *Aup, *Alo, *R, *A1, *A2, *L, *S ;
    cholmod_common Common, *cm ;
    int n, i, coletree, j, lnz, p, k, height, c ;
    char buf [LEN] ;

    /* ---------------------------------------------------------------------- */
    /* start CHOLMOD and set defaults */
    /* ---------------------------------------------------------------------- */

    cm = &Common ;
    cholmod_start (cm) ;
    sputil_config (SPUMONI, cm) ;

    /* ---------------------------------------------------------------------- */
    /* get inputs */
    /* ---------------------------------------------------------------------- */

    if (nargout > 5 || nargin < 1 || nargin > 3)
    {
	mexErrMsgTxt (
	    "Usage: [count h parent post R] = symbfact2 (A, mode, Lmode)") ;
    }

    /* ---------------------------------------------------------------------- */
    /* get input matrix A */
    /* ---------------------------------------------------------------------- */

    A = sputil_get_sparse_pattern (pargin [0], &Amatrix, &dummy, cm) ;
    S = (A == &Amatrix) ? NULL : A ;

    /* ---------------------------------------------------------------------- */
    /* get A->stype, default is to use triu(A) */
    /* ---------------------------------------------------------------------- */

    A->stype = 1 ;
    n = A->nrow ;
    coletree = FALSE ;
    if (nargin > 1)
    {
	buf [0] = '\0' ;
	if (mxIsChar (pargin [1]))
	{
	    mxGetString (pargin [1], buf, LEN) ;
	}
	c = buf [0] ;
	if (tolower (c) == 'r')
	{
	    /* unsymmetric case (A*A') if string starts with 'r' */
	    A->stype = 0 ;
	}
	else if (tolower (c) == 'c')
	{
	    /* unsymmetric case (A'*A) if string starts with 'c' */
	    n = A->ncol ;
	    coletree = TRUE ;
	    A->stype = 0 ;
	}
	else if (tolower (c) == 's')
	{
	    /* symmetric upper case (A) if string starts with 's' */
	    A->stype = 1 ;
	}
	else if (tolower (c) == 'l')
	{
	    /* symmetric lower case (A) if string starts with 'l' */
	    A->stype = -1 ;
	}
	else
	{
	    mexErrMsgTxt ("symbfact2: unrecognized mode") ;
	}
    }

    if (A->stype && A->nrow != A->ncol)
    {
	mexErrMsgTxt ("symbfact2: A must be square") ;
    }

    /* ---------------------------------------------------------------------- */
    /* compute the etree, its postorder, and the row/column counts */
    /* ---------------------------------------------------------------------- */

    Parent = cholmod_malloc (n, sizeof (int), cm) ;
    Post = cholmod_malloc (n, sizeof (int), cm) ;
    ColCount = cholmod_malloc (n, sizeof (int), cm) ;
    First = cholmod_malloc (n, sizeof (int), cm) ;
    Level = cholmod_malloc (n, sizeof (int), cm) ;

    /* F = A' */
    F = cholmod_transpose (A, 0, cm) ;

    if (A->stype == 1 || coletree)
    {
	/* symmetric upper case: find etree of A, using triu(A) */
	/* column case: find column etree of A, which is etree of A'*A */
	Aup = A ;
	Alo = F ;
    }
    else
    {
	/* symmetric lower case: find etree of A, using tril(A) */
	/* row case: find row etree of A, which is etree of A*A' */
	Aup = F ;
	Alo = A ;
    }

    cholmod_etree (Aup, Parent, cm) ;

    if (cm->status < CHOLMOD_OK)
    {
	/* out of memory or matrix invalid */
	mexErrMsgTxt ("symbfact2 failed: matrix corrupted!") ;
    }

    if (cholmod_postorder (Parent, n, NULL, Post, cm) != n)
    {
	/* out of memory or Parent invalid */
	mexErrMsgTxt ("symbfact2 postorder failed!") ;
    }

    /* symmetric upper case: analyze tril(F), which is triu(A) */
    /* column case: analyze F*F', which is A'*A */
    /* symmetric lower case: analyze tril(A) */
    /* row case: analyze A*A' */
    cholmod_rowcolcounts (Alo, NULL, 0, Parent, Post, NULL, ColCount,
		First, Level, cm) ;

    if (cm->status < CHOLMOD_OK)
    {
	/* out of memory or matrix invalid */
	mexErrMsgTxt ("symbfact2 failed: matrix corrupted!") ;
    }

    /* ---------------------------------------------------------------------- */
    /* return results to MATLAB: count, h, parent, and post */
    /* ---------------------------------------------------------------------- */

    pargout [0] = sputil_put_int (ColCount, n, 0) ;
    if (nargout > 1)
    {
	/* compute the elimination tree height */
	height = 0 ;
	for (i = 0 ; i < n ; i++)
	{
	    height = MAX (height, Level [i]) ;
	}
	height++ ;
	pargout [1] = mxCreateDoubleScalar ((double) height) ;
    }
    if (nargout > 2)
    {
	pargout [2] = sputil_put_int (Parent, n, 1) ;
    }
    if (nargout > 3)
    {
	pargout [3] = sputil_put_int (Post, n, 1) ;
    }

    /* ---------------------------------------------------------------------- */
    /* construct L, if requested */
    /* ---------------------------------------------------------------------- */

    if (nargout > 4)
    {

	if (A->stype == 1)
	{
	    /* symmetric upper case: use triu(A) only, A2 not needed */
	    A1 = A ;
	    A2 = NULL ;
	}
	else if (A->stype == -1)
	{
	    /* symmetric lower case: use tril(A) only, A2 not needed */
	    A1 = F ;
	    A2 = NULL ;
	}
	else if (coletree)
	{
	    /* column case: analyze F*F' */
	    A1 = F ;
	    A2 = A ;
	}
	else
	{
	    /* row case: analyze A*A' */
	    A1 = A ;
	    A2 = F ;
	}

	/* count the total number of entries in L */
	lnz = 0 ;
	for (j = 0 ; j < n ; j++)
	{
	    lnz += ColCount [j] ;
	}

	/* allocate the output matrix L (pattern-only) */
	L = cholmod_allocate_sparse (n, n, lnz, TRUE, TRUE, 0, CHOLMOD_PATTERN,
		cm) ;
	Lp = L->p ;
	Li = L->i ;

	/* initialize column pointers */
	lnz = 0 ;
	for (j = 0 ; j < n ; j++)
	{
	    Lp [j] = lnz ;
	    lnz += ColCount [j] ;
	}
	Lp [j] = lnz ;

	/* create a copy of the column pointers */
	W = First ;
	for (j = 0 ; j < n ; j++)
	{
	    W [j] = Lp [j] ;
	}

	/* get workspace for computing one row of L */
	R = cholmod_allocate_sparse (n, 1, n, FALSE, TRUE, 0, CHOLMOD_PATTERN,
		cm) ;
	Rp = R->p ;
	Ri = R->i ;

	/* compute L one row at a time */
	for (k = 0 ; k < n ; k++)
	{
	    /* get the kth row of L and store in the columns of L */
	    cholmod_row_subtree (A1, A2, k, Parent, R, cm) ;
	    for (p = 0 ; p < Rp [1] ; p++)
	    {
		Li [W [Ri [p]]++] = k ;
	    }
	    /* add the diagonal entry */
	    Li [W [k]++] = k ;
	}

	/* free workspace */
	cholmod_free_sparse (&R, cm) ;

	/* transpose L to get R, or leave as is */
	if (nargin < 3)
	{
	    /* R = L' */
	    R = cholmod_transpose (L, 0, cm) ;
	    cholmod_free_sparse (&L, cm) ;
	    L = R ;
	}

	/* fill numerical values of L with one's (only MATLAB needs this...) */
	L->x = cholmod_malloc (lnz, sizeof (double), cm) ;
	Lx = L->x ;
	for (p = 0 ; p < lnz ; p++)
	{
	    Lx [p] = 1 ;
	}
	L->xtype = CHOLMOD_REAL ;

	/* return L (or R) to MATLAB */
	pargout [4] = sputil_put_sparse (&L, cm) ;
    }

    /* ---------------------------------------------------------------------- */
    /* free workspace */
    /* ---------------------------------------------------------------------- */

    cholmod_free (n, sizeof (int), Parent, cm) ;
    cholmod_free (n, sizeof (int), Post, cm) ;
    cholmod_free (n, sizeof (int), ColCount, cm) ;
    cholmod_free (n, sizeof (int), First, cm) ;
    cholmod_free (n, sizeof (int), Level, cm) ;
    cholmod_free_sparse (&F, cm) ;
    cholmod_free_sparse (&S, cm) ;
    cholmod_finish (cm) ;
    cholmod_print_common (" ", cm) ;
    /*
    if (cm->malloc_count != ((nargout == 5) ? 3:0)) mexErrMsgTxt ("!") ;
    */
}