/* ========================================================================== */
/* === symamd mexFunction =================================================== */
/* ========================================================================== */

/*
    Usage:

	P = symamd (A) ;

    Returns a permutation vector P such that the Cholesky factorization of
    A (P,P) tends to be sparser than that of A.  This routine provides the same
    functionality as SYMMMD, but tends to be much faster and tends to return a
    better permutation vector.  Note that the SYMMMD m-file in
    Matlab 5.2 also performs a symmetric elimination tree post-ordering.  This
    mexFunction does not do this post-ordering.  This mexFunction is a
    replacement for the p = sparsfun ('symmmd', A) operation.

    A must be square, and is assummed to have a symmetric nonzero pattern.
    Only the nonzero pattern of the lower triangular portion of A is accessed.
    This routine constructs a matrix M such that the nonzero pattern of M'M is
    equal to A (assuming A has symmetric pattern), and then performs a column
    ordering of M using COLAMD.

    Uses the parameter spparms ('wh_frac') to determine if any "dense" columns
    of A are to be ignored.

    Calls/uses the following Matlab 5.2 routines/definitions/include files:

	mex.h			matrix.h		mexFunction
	mexCallMATLAB		mexErrMsgTxt		mexPrintf
	mxArray			mxCalloc		mxCreateDoubleMatrix
	mxCreateString		mxDestroyArray		mxFree
	mxGetIr			mxGetJc			mxGetM
	mxGetN			mxGetNumberOfDimensions	mxGetPr
	mxGetScalar		mxIsSparse		mxREAL

    Calls/uses the following ANSI C library routines/definitions/include files:

	stdlib.h

    Calls/uses the following routines/defintions in colamd.c and colamd.h:

	colamd.h		colamd_set_defaults	colamd_recommended
	colamd			COLAMD_KNOBS
	COLAMD_DENSE_ROW	COLAMD_DENSE_COL	COLAMD_DEFRAG_COUNT

    Authors:

	The authors of the code itself are Stefan I. Larimore and Timothy A.
	Davis (davis@cise.ufl.edu), University of Florida.  The algorithm was
	developed in collaboration with John Gilbert, Xerox PARC, and Esmond
	Ng, Oak Ridge National Laboratory.

    Date:

	August 3, 1998.  Version 1.0.  Written under Matlab 5.2.0.3084, and
	tested in Solaris 2.6.

    Acknowledgements:

	This work was supported by the National Science Foundation, under
	grants DMS-9504974 and DMS-9803599.

    Notice:

	Copyright (c) 1998 by the University of Florida.  All Rights Reserved.

	THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
	EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.

	Permission is hereby granted to use or copy this program for any
	purpose, provided the above notices are retained on all copies.
	User documentation of any code that uses this code must cite the
	Authors, the Copyright, and "Used by permission."  If this code is
	accessible from within Matlab, then typing "help colamd" or "colamd"
	(with no arguments) must cite the Authors.  Permission to modify the
	code and to distribute modified code is granted, provided the above
	notices are retained, and a notice that the code was modified is
	included with the above copyright notice.  You must also retain the
	Availability information below, of the original version.

	This software is provided free of charge.

    Availability:

	This file is located at

		http://www.cise.ufl.edu/~davis/colamd/symamdmex.c

	The colamd.h and colamd.c files are required, located in the same
	directory.  All codes are purely ANSI C compliant.
*/

/* ========================================================================== */
/* === Include files and prototypes ========================================= */
/* ========================================================================== */

#include "colamd.h"
#include "mex.h"
#include "matrix.h"
#include <stdlib.h>

static void symamd_help (void) ;

/* ========================================================================== */
/* === symamd mexFunction =================================================== */
/* ========================================================================== */

void mexFunction
(
    /* === Parameters ======================================================= */

    int nlhs,			/* number of left-hand sides */
    mxArray *plhs [],		/* left-hand side matrices */
    int nrhs,			/* number of right--hand sides */
    const mxArray *prhs []	/* right-hand side matrices */
)
{
    /* === Local variables ================================================== */

    int Mlen ;			/* size of M */
    int m_row ;			/* number of rows in M */
    int *Mp ;			/* column pointers for M */
    int *count ;		/* column counts for M, and temp. col pointers*/
    int *Ai ;			/* row indices of input matrix A */
    int *Ap ;			/* column pointers of input matrix A */
    int j ;			/* loop counter */
    int k ;			/* a row index of M */
    int *M ;			/* row indices of M, and colamd workspace */
    int ptr ;			/* a column "pointer" */
    int mnz ;			/* number of nonzeros in M */
    int n_col ;			/* number of columns of A */
    int n_row ;			/* number of rows of A */
    int full ;			/* TRUE if input matrix full, FALSE if sparse */
    double knobs [COLAMD_KNOBS] ; /* colamd user-controllable parameters */
    double *out_perm ;		/* output permutation vector */
    int i ;			/* loop counter */
    mxArray *Ainput ;		/* input matrix handle */
    mxArray *string ;		/* name of matlab function mexCallMATLAB */
    mxArray *parameter ;	/* return matrix from mexCallMATLAB */
    int spumoni ;		/* matlab variable for verbosity */

    /* === Get spumoni ====================================================== */

    string = mxCreateString ("spumoni") ;
    mexCallMATLAB (1, &parameter, 1, &string, "spparms") ;
    spumoni = (int) mxGetScalar (parameter) ;
    mxDestroyArray (parameter) ;
    mxDestroyArray (string) ;

    /* === Check inputs ===================================================== */

    if (nlhs < 0 || nlhs > 1 || nrhs != 1)
    {
	symamd_help () ;
	mexErrMsgTxt (
	"symamd: incorrect number of input and/or output arguments.") ;
    }

    /* === Set default knobs for colamd ===================================== */

    colamd_set_defaults (knobs) ;

    /* use spparms ('wh_frac') for the dense column setting */
    string = mxCreateString ("wh_frac") ;
    mexCallMATLAB (1, &parameter, 1, &string, "spparms") ;
    knobs [COLAMD_DENSE_COL] = mxGetScalar (parameter) ;

    /* there are no dense rows in M */
    knobs [COLAMD_DENSE_ROW] = 1.0 ;

    /* print knob settings if spumoni > 0 */
    if (spumoni > 0)
    {
	symamd_help () ;
	mexPrintf ("symamd:  dense row/col fraction: %f\n",
	    knobs [COLAMD_DENSE_COL]) ;
    }

    /* === If A is full, convert to a sparse matrix ========================= */

    Ainput = (mxArray *) prhs [0] ;
    if (mxGetNumberOfDimensions (Ainput) != 2)
    {
	symamd_help () ;
	mexErrMsgTxt ("symamd: input matrix must be 2-dimensional.") ;
    }
    full = !mxIsSparse (Ainput) ;
    if (full)
    {
	mexCallMATLAB (1, &Ainput, 1, (mxArray **) prhs, "sparse") ;
    }

    /* === Allocate workspace for colamd and construct M ==================== */

    /* get size of input matrix A */
    n_col = mxGetN (Ainput) ;
    n_row = mxGetM (Ainput) ;
    if (n_col != n_row)
    {
	symamd_help () ;
	mexErrMsgTxt ("symamd: matrix must be square.") ;
    }

    Ai = mxGetIr (Ainput) ;
    Ap = mxGetJc (Ainput) ;

    count = (int *) mxCalloc (n_col+1, sizeof (int)) ;
    Mp = (int *) mxCalloc (n_col+1, sizeof (int)) ;

    for (j = 0 ; j < n_col ; j++)
    {
	for (ptr = Ap [j] ; ptr < Ap [j+1] ; ptr++)
	{
	    i = Ai [ptr] ;
	    if (i > j)
	    {
		/* row k of M will contain column indices i and j */
		count [i]++ ;
		count [j]++ ;
	    }
	}
    }

    Mp [0] = 0 ;
    for (j = 1 ; j <= n_col ; j++)
    {
	Mp [j] = Mp [j-1] + count [j-1] ;
    }
    for (j = 0 ; j < n_col ; j++)
    {
	count [j] = Mp [j] ;
    }

    mnz = Mp [n_col] ;
    m_row = mnz / 2 ;
    Mlen = colamd_recommended (mnz, m_row, n_col) ;
    M = (int *) mxCalloc (Mlen, sizeof (int)) ;

    k = 0 ;
    for (j = 0 ; j < n_col ; j++)
    {
	for (ptr = Ap [j] ; ptr < Ap [j+1] ; ptr++)
	{
	    i = Ai [ptr] ;
	    if (i > j)
	    {
		/* row k of M contains column indices i and j */
		M [count [i]++] = k ;
		M [count [j]++] = k ;
		k++ ;
	    }
	}
    }

    mxFree (count) ;
    if (full)
    {
	mxDestroyArray (Ainput) ;
    }

    /* === Order the columns (destroys M) =================================== */

    if (!colamd (m_row, n_col, Mlen, M, Mp, knobs))
    {
	symamd_help () ;
	mexErrMsgTxt
	("symamd: internal error!  Please report this to the authors.\n") ;
    }

    /* === Return the permutation vector ==================================== */

    plhs [0] = mxCreateDoubleMatrix (1, n_col, mxREAL) ;
    out_perm = mxGetPr (plhs [0]) ;
    for (i = 0 ; i < n_col ; i++)
    {
	/* colamd is 0-based, but Matlab expects this to be 1-based */
	out_perm [i] = Mp [i] + 1 ;
    }
    mxFree (Mp) ;

    /* === Print the statistics ============================================= */

    if (spumoni > 0)
    {
	mexPrintf ("symamd:  dense columns ignored:         %d\n",
		M [COLAMD_DENSE_COL]) ;
	mexPrintf ("symamd:  number of garbage collections: %d\n",
		M [COLAMD_DEFRAG_COUNT]) ;
    }
    mxFree (M) ;
}


/* ========================================================================== */
/* === symamd_help ========================================================== */
/* ========================================================================== */

/* NOTICE:  You do not have permission to delete the Authors' names from this */
/* function.  Typing "help symamd" or just "symamd" must cite the Authors. */

static void symamd_help (void)
{
    mexPrintf ("\n") ;
    mexPrintf (" SYMAMD Symmetric approximate minimum degree permutation.\n") ;
    mexPrintf ("    P = SYMAMD (S), for a symmetric positive definite") ;
    mexPrintf (" matrix S,\n    returns the permutation vector p such that") ;
    mexPrintf (" S(p,p) tends to have a\n    sparser Cholesky factor than") ;
    mexPrintf (" S.  Sometimes SYMAMD works well\n    for symmetric") ;
    mexPrintf (" indefinite matrices too.  SYMAMD tends to be faster\n") ;
    mexPrintf ("    than SYMMMD and tends to return a better ordering.\n") ;
    mexPrintf ("    SYMAMD uses spparms ('wh_frac') to determine how to\n") ;
    mexPrintf ("    treat dense rows and columns.\n") ;
    mexPrintf ("\n    See also SYMAMDTREE, COLAMD, COLAMDTREE, COLMMD,") ;
    mexPrintf (" SYMMMD, SYMRCM, COLPERM.\n\n") ;
    mexPrintf ("    Authors:  Stefan I. Larimore and Timothy A. Davis,") ;
    mexPrintf (" University of Florida,\n") ;
    mexPrintf ("    (davis@cise.ufl.edu); in collaboration with John Gilbert") ;
    mexPrintf (", Xerox PARC, and\n    Esmond Ng, Oak Ridge National") ;
    mexPrintf (" Laboratory.  This work was supported by the\n    National") ;
    mexPrintf (" Science Foundation, under grants DMS-9504974 and") ;
    mexPrintf (" DMS-9803599.\n    COLAMD and SYMAMD are available") ;
    mexPrintf (" at http://www.cise.ufl.edu/~davis/colamd.\n\n") ;
}