/* ========================================================================== */
/* === UMFPACK mexFunction ================================================== */
/* ========================================================================== */

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

/*
    MATLAB interface for umfpack.

    Factor or solve a sparse linear system, returning either the solution
    x to Ax=b or A'x'=b', or the factorization LU=P(R\A)Q or LU=PAQ.  A must be
    sparse, with nonzero dimensions, but it may be complex, singular, and/or
    rectangular.  b must be a dense n-by-1 vector (real or complex).
    L is unit lower triangular, U is upper triangular, and R is diagonal.
    P and Q are permutation matrices (permutations of an identity matrix).

    The matrix A is scaled, by default.  Each row i is divided by r (i), where
    r (i) is the sum of the absolute values of the entries in that row.  The
    scaled matrix has an infinity norm of 1.  The scale factors r (i) are
    returned in a diagonal sparse matrix.  If the factorization is:

	[L, U, P, Q, R] = umfpack (A) ;

    then the factorization is

	L*U = P * (R \ A) * Q

    This is safer than returning a matrix R such that L*U = P*R*A*Q, because
    it avoids the division by small entries.  If r(i) is subnormal, multiplying
    by 1/r(i) would result in an IEEE Infinity, but dividing by r(i) is safe.

    The factorization

	[L, U, P, Q] = umfpack (A) ;

    returns LU factors such that L*U = P*A*Q, with no scaling.

    See umfpack.m, umfpack_details.m, and umfpack.h for details.

    Note that this mexFunction accesses only the user-callable UMFPACK routines.
    Thus, is also provides another example of how user C code can access
    UMFPACK.

    If NO_TRANSPOSE_FORWARD_SLASH is not defined at compile time, then the
    forward slash (/) operator acts almost like x = b/A in MATLAB 6.1.  It is
    solved by factorizing the array transpose, and then x = (A.'\b.').' is
    solved.  This is the default behavior (for historical reasons), since
    factorizing A can behave perform much differently than factorizing its
    transpose.

    If NO_TRANSPOSE_FORWARD_SLASH is defined at compile time, then the forward
    slash operator does not act like x=b/A in MATLAB 6.1.  It is solved by
    factorizing A, and then solving via the transposed L and U matrices.
    The solution is still x = (A.'\b.').', except that A is factorized instead
    of A.'.

    Modified for v4.3.1, Jan 10, 2005: default has been changed to
    NO_TRANSPOSE_FORWARD_SLASH, to test iterative refinement for b/A.
    v4.4: added method for computing the determinant.
*/
#define NO_TRANSPOSE_FORWARD_SLASH  /* default has changed for v4.3.1 */

#include "umfpack.h"
#include "mex.h"
#include "matrix.h"
#include <string.h>
#include <math.h>
#include <float.h>

#define MIN(a,b) (((a) < (b)) ? (a) : (b))
#define MAX(a,b) (((a) > (b)) ? (a) : (b))
#define STRING_MATCH(s1,s2) (strcmp ((s1), (s2)) == 0)
#ifndef TRUE
#define TRUE (1)
#endif
#ifndef FALSE
#define FALSE (0)
#endif

/* ========================================================================== */
/* === error ================================================================ */
/* ========================================================================== */

/* Return an error message */

static void error
(
    char *s,
    int A_is_complex,
    int nargout,
    mxArray *pargout [ ],
    double Control [UMFPACK_CONTROL],
    double Info [UMFPACK_INFO],
    int status,
    int do_info
)
{
    int i ;
    double *OutInfo ;
    if (A_is_complex)
    {
	umfpack_zi_report_status (Control, status) ;
	umfpack_zi_report_info (Control, Info) ;
    }
    else
    {
	umfpack_di_report_status (Control, status) ;
	umfpack_di_report_info (Control, Info) ;
    }
    if (do_info > 0)
    {
	/* return Info */
	pargout [do_info] = mxCreateDoubleMatrix (1, UMFPACK_INFO, mxREAL) ;
	OutInfo = mxGetPr (pargout [do_info]) ;
	for (i = 0 ; i < UMFPACK_INFO ; i++)
	{
	    OutInfo [i] = Info [i] ;
	}
    }
    mexErrMsgTxt (s) ;
}


/* ========================================================================== */
/* === UMFPACK ============================================================== */
/* ========================================================================== */

void mexFunction
(
    int nargout,		/* number of outputs */
    mxArray *pargout [ ],	/* output arguments */
    int nargin,			/* number of inputs */
    const mxArray *pargin [ ]	/* input arguments */
)
{

    /* ---------------------------------------------------------------------- */
    /* local variables */
    /* ---------------------------------------------------------------------- */

    double Info [UMFPACK_INFO], Control [UMFPACK_CONTROL], dx, dz, dexp ;
    double *Lx, *Lz, *Ux, *Uz, *Ax, *Az, *Bx, *Bz, *Xx, *Xz, *User_Control,
	*p, *q, *OutInfo, *p1, *p2, *p3, *p4, *Ltx, *Ltz, *Rs, *Px, *Qx ;
    void *Symbolic, *Numeric ;
    int *Lp, *Li, *Up, *Ui, *Ap, *Ai, *P, *Q, do_solve, lnz, unz, nn, i,
	transpose, size, do_info, do_numeric, *Front_npivcol, op, k, *Rp, *Ri,
	*Front_parent, *Chain_start, *Chain_maxrows, *Chain_maxcols, nz, status,
	nfronts, nchains, *Ltp, *Ltj, *Qinit, print_level, status2, no_scale,
	*Front_1strow, *Front_leftmostdesc, n_row, n_col, n_inner, sys,
	ignore1, ignore2, ignore3, A_is_complex, B_is_complex, X_is_complex,
	*Pp, *Pi, *Qp, *Qi, do_recip, do_det ;
    mxArray *Amatrix, *Bmatrix, *User_Control_matrix, *User_Qinit ;
    char *operator, *operation ;
    mxComplexity Atype, Xtype ;
    char warning [200] ;

#ifndef NO_TRANSPOSE_FORWARD_SLASH
    int *Cp, *Ci ;
    double *Cx, *Cz ;
#endif

    /* ---------------------------------------------------------------------- */
    /* get inputs A, b, and the operation to perform */
    /* ---------------------------------------------------------------------- */

    User_Control_matrix = (mxArray *) NULL ;
    User_Qinit = (mxArray *) NULL ;

    do_info = 0 ;
    do_solve = FALSE ;
    do_numeric = TRUE ;
    transpose = FALSE ;
    no_scale = FALSE ;
    do_det = FALSE ;

    /* find the operator */
    op = 0 ;
    for (i = 0 ; i < nargin ; i++)
    {
	if (mxIsChar (pargin [i]))
	{
	    op = i ;
	    break ;
	}
    }

    if (op > 0)
    {
	operator = mxArrayToString (pargin [op]) ;

	if (STRING_MATCH (operator, "\\"))
	{

	    /* -------------------------------------------------------------- */
	    /* matrix left divide, x = A\b */
	    /* -------------------------------------------------------------- */

	    /*
		[x, Info] = umfpack (A, '\', b) ;
		[x, Info] = umfpack (A, '\', b, Control) ;
		[x, Info] = umfpack (A, Qinit, '\', b, Control) ;
		[x, Info] = umfpack (A, Qinit, '\', b) ;
	    */

	    operation = "x = A\\b" ;
	    do_solve = TRUE ;
	    Amatrix = (mxArray *) pargin [0] ;
	    Bmatrix = (mxArray *) pargin [op+1] ;

	    if (nargout == 2)
	    {
		do_info = 1 ;
	    }
	    if (op == 2)
	    {
		User_Qinit = (mxArray *) pargin [1] ;
	    }
	    if ((op == 1 && nargin == 4) || (op == 2 && nargin == 5))
	    {
		User_Control_matrix = (mxArray *) pargin [nargin-1] ;
	    }
	    if (nargin < 3 || nargin > 5 || nargout > 2)
	    {
		mexErrMsgTxt ("wrong number of arguments") ;
	    }

	}
	else if (STRING_MATCH (operator, "/"))
	{

	    /* -------------------------------------------------------------- */
	    /* matrix right divide, x = b/A */
	    /* -------------------------------------------------------------- */

	    /*
		[x, Info] = umfpack (b, '/', A) ;
		[x, Info] = umfpack (b, '/', A, Control) ;
		[x, Info] = umfpack (b, '/', A, Qinit) ;
		[x, Info] = umfpack (b, '/', A, Qinit, Control) ;
	    */

	    operation = "x = b/A" ;
	    do_solve = TRUE ;
	    transpose = TRUE ;
	    Amatrix = (mxArray *) pargin [2] ;
	    Bmatrix = (mxArray *) pargin [0] ;

	    if (nargout == 2)
	    {
		do_info = 1 ;
	    }
	    if (nargin == 5)
	    {
		User_Qinit = (mxArray *) pargin [3] ;
		User_Control_matrix = (mxArray *) pargin [4] ;
	    }
	    else if (nargin == 4)
	    {
		/* Control is k-by-1 where k > 1, Qinit is 1-by-n */
		if (mxGetM (pargin [3]) == 1)
		{
		    User_Qinit = (mxArray *) pargin [3] ;
		}
		else
		{
		    User_Control_matrix = (mxArray *) pargin [3] ;
		}
	    }
	    else if (nargin < 3 || nargin > 5 || nargout > 2)
	    {
		mexErrMsgTxt ("wrong number of arguments") ;
	    }

	}
	else if (STRING_MATCH (operator, "symbolic"))
	{

	    /* -------------------------------------------------------------- */
	    /* symbolic factorization only */
	    /* -------------------------------------------------------------- */

	    /*
	    [P Q Fr Ch Info] = umfpack (A, 'symbolic') ;
	    [P Q Fr Ch Info] = umfpack (A, 'symbolic', Control) ;
	    [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic') ;
	    [P Q Fr Ch Info] = umfpack (A, Qinit, 'symbolic', Control) ;
	    */

	    operation = "symbolic factorization" ;
	    do_numeric = FALSE ;
	    Amatrix = (mxArray *) pargin [0] ;

	    if (nargout == 5)
	    {
		do_info = 4 ;
	    }
	    if (op == 2)
	    {
		User_Qinit = (mxArray *) pargin [1] ;
	    }
	    if ((op == 1 && nargin == 3) || (op == 2 && nargin == 4))
	    {
		User_Control_matrix = (mxArray *) pargin [nargin-1] ;
	    }
	    if (nargin < 2 || nargin > 4 || nargout > 5 || nargout < 4)
	    {
		mexErrMsgTxt ("wrong number of arguments") ;
	    }

	}
	else if (STRING_MATCH (operator, "det"))
	{

	    /* -------------------------------------------------------------- */
	    /* compute the determinant */
	    /* -------------------------------------------------------------- */

	    /*
	     * [det] = umfpack (A, 'det') ;
	     * [dmantissa dexp] = umfpack (A, 'det') ;
	     */

	    operation = "determinant" ;
	    do_det = TRUE ;
	    Amatrix = (mxArray *) pargin [0] ;
	    if (nargin > 2 || nargout > 2)
	    {
		mexErrMsgTxt ("wrong number of arguments") ;
	    }

	}
	else
	{
	    mexErrMsgTxt ("operator must be '/', '\\', or 'symbolic'") ;
	}
	mxFree (operator) ;

    }
    else if (nargin > 0)
    {

	/* ------------------------------------------------------------------ */
	/* LU factorization */
	/* ------------------------------------------------------------------ */

	/*
	    with scaling:
	    [L, U, P, Q, R, Info] = umfpack (A) ;
	    [L, U, P, Q, R, Info] = umfpack (A, Qinit) ;

	    scaling determined by Control settings:
	    [L, U, P, Q, R, Info] = umfpack (A, Control) ;
	    [L, U, P, Q, R, Info] = umfpack (A, Qinit, Control) ;

	    with no scaling:
	    [L, U, P, Q] = umfpack (A) ;
	    [L, U, P, Q] = umfpack (A, Control) ;
	    [L, U, P, Q] = umfpack (A, Qinit) ;
	    [L, U, P, Q] = umfpack (A, Qinit, Control) ;
	*/

	operation = "numeric factorization" ;
	Amatrix = (mxArray *) pargin [0] ;

	no_scale = nargout <= 4 ;

	if (nargout == 6)
	{
	    do_info = 5 ;
	}
	if (nargin == 3)
	{
	    User_Qinit = (mxArray *) pargin [1] ;
	    User_Control_matrix = (mxArray *) pargin [2] ;
	}
	else if (nargin == 2)
	{
	    /* Control is k-by-1 where k > 1, Qinit is 1-by-n */
	    if (mxGetM (pargin [1]) == 1)
	    {
		User_Qinit = (mxArray *) pargin [1] ;
	    }
	    else
	    {
		User_Control_matrix = (mxArray *) pargin [1] ;
	    }
	}
	else if (nargin > 3 || nargout > 6 || nargout < 4)
	{
	    mexErrMsgTxt ("wrong number of arguments") ;
	}

    }
    else
    {

	/* ------------------------------------------------------------------ */
	/* return default control settings */
	/* ------------------------------------------------------------------ */

	/*
	    Control = umfpack ;
	    umfpack ;
	*/

	if (nargout > 1)
	{
	    mexErrMsgTxt ("wrong number of arguments") ;
	}

	pargout [0] = mxCreateDoubleMatrix (UMFPACK_CONTROL, 1, mxREAL) ;
	User_Control = mxGetPr (pargout [0]) ;
	umfpack_di_defaults (User_Control) ;

	return ;
    }

    /* ---------------------------------------------------------------------- */
    /* check inputs */
    /* ---------------------------------------------------------------------- */

    if (mxGetNumberOfDimensions (Amatrix) != 2)
    {
	mexErrMsgTxt ("input matrix A must be 2-dimensional") ;
    }
    n_row = mxGetM (Amatrix) ;
    n_col = mxGetN (Amatrix) ;
    nn = MAX (n_row, n_col) ;
    n_inner = MIN (n_row, n_col) ;
    if (do_solve && n_row != n_col)
    {
	mexErrMsgTxt ("input matrix A must square for '\\' or '/'") ;
    }
    if (!mxIsSparse (Amatrix))
    {
	mexErrMsgTxt ("input matrix A must be sparse") ;
    }
    if (n_row == 0 || n_col == 0)
    {
	mexErrMsgTxt ("input matrix A cannot have zero rows or zero columns") ;
    }

    /* The real/complex status of A determines which version to use, */
    /* (umfpack_di_* or umfpack_zi_*). */
    A_is_complex = mxIsComplex (Amatrix) ;
    Atype = A_is_complex ? mxCOMPLEX : mxREAL ;
    Ap = mxGetJc (Amatrix) ;
    Ai = mxGetIr (Amatrix) ;
    Ax = mxGetPr (Amatrix) ;
    Az = mxGetPi (Amatrix) ;

    if (do_solve)
    {

	if (n_row != n_col)
	{
	    mexErrMsgTxt ("A must be square for \\ or /") ;
	}
	if (transpose)
	{
	    if (mxGetM (Bmatrix) != 1 || mxGetN (Bmatrix) != nn)
	    {
		mexErrMsgTxt ("b has the wrong dimensions") ;
	    }
	}
	else
	{
	    if (mxGetM (Bmatrix) != nn || mxGetN (Bmatrix) != 1)
	    {
		mexErrMsgTxt ("b has the wrong dimensions") ;
	    }
	}
	if (mxGetNumberOfDimensions (Bmatrix) != 2)
	{
	    mexErrMsgTxt ("input matrix b must be 2-dimensional") ;
	}
	if (mxIsSparse (Bmatrix))
	{
	    mexErrMsgTxt ("input matrix b cannot be sparse") ;
	}
	if (mxGetClassID (Bmatrix) != mxDOUBLE_CLASS)
	{
	    mexErrMsgTxt ("input matrix b must double precision matrix") ;
	}

	B_is_complex = mxIsComplex (Bmatrix) ;
	Bx = mxGetPr (Bmatrix) ;
	Bz = mxGetPi (Bmatrix) ;

	X_is_complex = A_is_complex || B_is_complex ;
	Xtype = X_is_complex ? mxCOMPLEX : mxREAL ;
    }

    /* ---------------------------------------------------------------------- */
    /* set the Control parameters */
    /* ---------------------------------------------------------------------- */

    if (A_is_complex)
    {
	umfpack_zi_defaults (Control) ;
    }
    else
    {
	umfpack_di_defaults (Control) ;
    }
    if (User_Control_matrix)
    {
	if (mxGetClassID (User_Control_matrix) != mxDOUBLE_CLASS ||
	    mxIsSparse (User_Control_matrix))
	{
	    mexErrMsgTxt ("Control must be a dense real matrix") ;
	}
	size = UMFPACK_CONTROL ;
	size = MIN (size, mxGetNumberOfElements (User_Control_matrix)) ;
	User_Control = mxGetPr (User_Control_matrix) ;
	for (i = 0 ; i < size ; i++)
	{
	    Control [i] = User_Control [i] ;
	}
    }

    if (no_scale)
    {
	/* turn off scaling for [L, U, P, Q] = umfpack (A) ;
	 * ignoring the input value of Control (24) for the usage
	 * [L, U, P, Q] = umfpack (A, Control) ; */
	Control [UMFPACK_SCALE] = UMFPACK_SCALE_NONE ;
    }

    if (mxIsNaN (Control [UMFPACK_PRL]))
    {
	print_level = UMFPACK_DEFAULT_PRL ;
    }
    else
    {
	print_level = (int) Control [UMFPACK_PRL] ;
    }

    Control [UMFPACK_PRL] = print_level ;

    /* ---------------------------------------------------------------------- */
    /* get Qinit, if present */
    /* ---------------------------------------------------------------------- */

    if (User_Qinit)
    {
	if (mxGetM (User_Qinit) != 1 || mxGetN (User_Qinit) != n_col)
	{
	    mexErrMsgTxt ("Qinit must be 1-by-n_col") ;
	}
	if (mxGetNumberOfDimensions (User_Qinit) != 2)
	{
	    mexErrMsgTxt ("input Qinit must be 2-dimensional") ;
	}
	if (mxIsComplex (User_Qinit))
	{
	    mexErrMsgTxt ("input Qinit must not be complex") ;
	}
	if (mxGetClassID (User_Qinit) != mxDOUBLE_CLASS)
	{
	    mexErrMsgTxt ("input Qinit must be a double matrix") ;
	}
	if (mxIsSparse (User_Qinit))
	{
	    mexErrMsgTxt ("input Qinit must be dense") ;
	}
	Qinit = (int *) mxMalloc (n_col * sizeof (int)) ;
	p = mxGetPr (User_Qinit) ;
	for (k = 0 ; k < n_col ; k++)
	{
	    /* convert from 1-based to 0-based indexing */
	    Qinit [k] = ((int) (p [k])) - 1 ;
	}

    }
    else
    {
	/* umfpack_*_qsymbolic will call colamd to get Qinit. This is the */
	/* same as calling umfpack_*_symbolic with Qinit set to NULL*/
	Qinit = (int *) NULL ;
    }

    /* ---------------------------------------------------------------------- */
    /* report the inputs A and Qinit */
    /* ---------------------------------------------------------------------- */

    if (print_level >= 2)
    {
	/* print the operation */
	mexPrintf ("\numfpack: %s\n", operation) ;
    }

    if (A_is_complex)
    {
	umfpack_zi_report_control (Control) ;
	if (print_level >= 3) mexPrintf ("\nA: ") ;
	(void) umfpack_zi_report_matrix (n_row, n_col, Ap, Ai, Ax, Az,
	    1, Control) ;
	if (Qinit)
	{
	    if (print_level >= 3) mexPrintf ("\nQinit: ") ;
	    (void) umfpack_zi_report_perm (n_col, Qinit, Control) ;
	}
    }
    else
    {
	umfpack_di_report_control (Control) ;
	if (print_level >= 3) mexPrintf ("\nA: ") ;
	(void) umfpack_di_report_matrix (n_row, n_col, Ap, Ai, Ax,
	    1, Control) ;
	if (Qinit)
	{
	    if (print_level >= 3) mexPrintf ("\nQinit: ") ;
	    (void) umfpack_di_report_perm (n_col, Qinit, Control) ;
	}
    }

#ifndef NO_TRANSPOSE_FORWARD_SLASH
    /* ---------------------------------------------------------------------- */
    /* create the array transpose for x = b/A */
    /* ---------------------------------------------------------------------- */

    if (transpose)
    {
	/* note that in this case A will be square (nn = n_row = n_col) */
	/* x = (A.'\b.').' will be computed */

	/* make sure Ci and Cx exist, avoid malloc of zero-sized arrays. */
	nz = MAX (Ap [nn], 1) ;

	Cp = (int *) mxMalloc ((nn+1) * sizeof (int)) ;
	Ci = (int *) mxMalloc (nz * sizeof (int)) ;
	Cx = (double *) mxMalloc (nz * sizeof (double)) ;
	if (A_is_complex)
	{
	    Cz = (double *) mxMalloc (nz * sizeof (double)) ;
	    status = umfpack_zi_transpose (nn, nn, Ap, Ai, Ax, Az,
	        (int *) NULL, (int *) NULL, Cp, Ci, Cx, Cz, FALSE) ;
	}
	else
	{
	    status = umfpack_di_transpose (nn, nn, Ap, Ai, Ax,
	        (int *) NULL, (int *) NULL, Cp, Ci, Cx) ;
	}

	if (status != UMFPACK_OK)
	{
	    error ("transpose of A failed", A_is_complex, nargout, pargout,
		Control, Info, status, do_info);
	    return ;
	}

	/* modify pointers so that C will be factorized and solved, not A */
	Ap = Cp ;
	Ai = Ci ;
	Ax = Cx ;
	if (A_is_complex)
	{
	    Az = Cz ;
	}
    }
#endif

    /* ---------------------------------------------------------------------- */
    /* perform the symbolic factorization */
    /* ---------------------------------------------------------------------- */

    if (A_is_complex)
    {
	status = umfpack_zi_qsymbolic (n_row, n_col, Ap, Ai, Ax, Az,
	    Qinit, &Symbolic, Control, Info) ;
    }
    else
    {
	status = umfpack_di_qsymbolic (n_row, n_col, Ap, Ai, Ax,
	    Qinit, &Symbolic, Control, Info) ;
    }

    if (Qinit)
    {
	mxFree (Qinit) ;
    }

    if (status < 0)
    {
	error ("symbolic factorization failed", A_is_complex, nargout, pargout,
	    Control, Info, status, do_info) ;
	return ;
    }

    /* ---------------------------------------------------------------------- */
    /* report the Symbolic object */
    /* ---------------------------------------------------------------------- */

    if (A_is_complex)
    {
	(void) umfpack_zi_report_symbolic (Symbolic, Control) ;
    }
    else
    {
	(void) umfpack_di_report_symbolic (Symbolic, Control) ;
    }

    /* ---------------------------------------------------------------------- */
    /* perform numeric factorization, or just return symbolic factorization */
    /* ---------------------------------------------------------------------- */

    if (do_numeric)
    {

	/* ------------------------------------------------------------------ */
	/* perform the numeric factorization */
	/* ------------------------------------------------------------------ */

	if (A_is_complex)
	{
	    status = umfpack_zi_numeric (Ap, Ai, Ax, Az, Symbolic, &Numeric,
		Control, Info) ;
	}
	else
	{
	    status = umfpack_di_numeric (Ap, Ai, Ax, Symbolic, &Numeric,
		Control, Info) ;
	}

	/* ------------------------------------------------------------------ */
	/* free the symbolic factorization */
	/* ------------------------------------------------------------------ */

	if (A_is_complex)
	{
	    umfpack_zi_free_symbolic (&Symbolic) ;
	}
	else
	{
	    umfpack_di_free_symbolic (&Symbolic) ;
	}

	/* ------------------------------------------------------------------ */
	/* report the Numeric object */
	/* ------------------------------------------------------------------ */

	if (status < 0)
	{
	    error ("numeric factorization failed", A_is_complex, nargout,
		pargout, Control, Info, status, do_info);
	    return ;
	}

	if (A_is_complex)
	{
	    (void) umfpack_zi_report_numeric (Numeric, Control) ;
	}
	else
	{
	    (void) umfpack_di_report_numeric (Numeric, Control) ;
	}

	/* ------------------------------------------------------------------ */
	/* return the solution, determinant, or the factorization */
	/* ------------------------------------------------------------------ */

	if (do_solve)
	{
	    /* -------------------------------------------------------------- */
	    /* solve Ax=b or A'x'=b', and return just the solution x */
	    /* -------------------------------------------------------------- */

#ifndef NO_TRANSPOSE_FORWARD_SLASH
	    if (transpose)
	    {
		/* A.'x.'=b.' gives the same x=b/A as solving A'x'=b' */
		/* since C=A.' was factorized, solve with sys = UMFPACK_A */
		/* since x and b are vectors, x.' and b.' are implicit */
		pargout [0] = mxCreateDoubleMatrix (1, nn, Xtype) ;
	    }
	    else
	    {
		pargout [0] = mxCreateDoubleMatrix (nn, 1, Xtype) ;
	    }
	    sys = UMFPACK_A ;
#else
	    if (transpose)
	    {
		/* If A is real, A'x=b is the same as A.'x=b. */
		/* x and b are vectors, so x and b are the same as x' and b'. */
		/* If A is complex, then A.'x.'=b.' gives the same solution x */
		/* as the complex conjugate transpose.  If we used the A'x=b */
		/* option in umfpack_*_solve, we would have to form b' on */
		/* input and x' on output (negating the imaginary part). */
		/* We can save this work by just using the A.'x=b option in */
		/* umfpack_*_solve.  Then, forming x.' and b.' is implicit, */
		/* since x and b are just vectors anyway. */
		/* In both cases, the system to solve is A.'x=b */
		pargout [0] = mxCreateDoubleMatrix (1, nn, Xtype) ;
		sys = UMFPACK_Aat ;
	    }
	    else
	    {
		pargout [0] = mxCreateDoubleMatrix (nn, 1, Xtype) ;
		sys = UMFPACK_A ;
	    }
#endif

	    /* -------------------------------------------------------------- */
	    /* print the right-hand-side, B */
	    /* -------------------------------------------------------------- */

	    if (print_level >= 3) mexPrintf ("\nright-hand side, b: ") ;
	    if (B_is_complex)
	    {
		(void) umfpack_zi_report_vector (nn, Bx, Bz, Control) ;
	    }
	    else
	    {
		(void) umfpack_di_report_vector (nn, Bx, Control) ;
	    }

	    /* -------------------------------------------------------------- */
	    /* solve the system */
	    /* -------------------------------------------------------------- */

	    Xx = mxGetPr (pargout [0]) ;
	    Xz = mxGetPi (pargout [0]) ;
	    status2 = UMFPACK_OK ;

	    if (A_is_complex)
	    {
		if (!B_is_complex)
		{
		    /* umfpack_zi_solve expects a complex B */
		    Bz = (double *) mxCalloc (nn, sizeof (double)) ;
		}
		status = umfpack_zi_solve (sys, Ap, Ai, Ax, Az, Xx, Xz, Bx, Bz,
		    Numeric, Control, Info) ;
		if (!B_is_complex)
		{
		    mxFree (Bz) ;
		}
	    }
	    else
	    {
		if (B_is_complex)
		{
		    /* Ax=b when b is complex and A is sparse can be split */
		    /* into two systems, A*xr=br and A*xi=bi, where r denotes */
		    /* the real part and i the imaginary part of x and b. */
		    status2 = umfpack_di_solve (sys, Ap, Ai, Ax, Xz, Bz,
		    Numeric, Control, Info) ;
		}
		status = umfpack_di_solve (sys, Ap, Ai, Ax, Xx, Bx,
		    Numeric, Control, Info) ;
	    }

#ifndef NO_TRANSPOSE_FORWARD_SLASH
	    /* -------------------------------------------------------------- */
	    /* free the transposed matrix C */
	    /* -------------------------------------------------------------- */

	    if (transpose)
	    {
	        mxFree (Cp) ;
	        mxFree (Ci) ;
	        mxFree (Cx) ;
	        if (A_is_complex)
	        {
	            mxFree (Cz) ;
	        }
	    }
#endif

	    /* -------------------------------------------------------------- */
	    /* free the Numeric object */
	    /* -------------------------------------------------------------- */

	    if (A_is_complex)
	    {
		umfpack_zi_free_numeric (&Numeric) ;
	    }
	    else
	    {
		umfpack_di_free_numeric (&Numeric) ;
	    }

	    /* -------------------------------------------------------------- */
	    /* check error status */
	    /* -------------------------------------------------------------- */

	    if (status < 0 || status2 < 0)
	    {
		mxDestroyArray (pargout [0]) ;
		error ("solve failed", A_is_complex, nargout, pargout, Control,
			Info, status, do_info) ;
		return ;
	    }

	    /* -------------------------------------------------------------- */
	    /* print the solution, X */
	    /* -------------------------------------------------------------- */

	    if (print_level >= 3) mexPrintf ("\nsolution, x: ") ;
	    if (X_is_complex)
	    {
		(void) umfpack_zi_report_vector (nn, Xx, Xz, Control) ;
	    }
	    else
	    {
		(void) umfpack_di_report_vector (nn, Xx, Control) ;
	    }

	    /* -------------------------------------------------------------- */
	    /* warn about singular or near-singular matrices */
	    /* -------------------------------------------------------------- */

	    /* no warning is given if Control (1) is zero */

	    if (Control [UMFPACK_PRL] >= 1)
	    {
		if (status == UMFPACK_WARNING_singular_matrix)
		{
		    sprintf (warning, "matrix is singular\n"
			"Try increasing Control (%d) and Control (%d).\n"
			"(Suppress this warning with Control (%d) = 0.)\n",
			1+UMFPACK_PIVOT_TOLERANCE,
			1+UMFPACK_SYM_PIVOT_TOLERANCE,
			1+UMFPACK_PRL) ;
		    mexWarnMsgTxt (warning) ;
		}
		else if (Info [UMFPACK_RCOND] < DBL_EPSILON)
		{
		    sprintf (warning, "matrix is nearly singular, rcond = %g\n"
			"Try increasing Control (%d) and Control (%d).\n"
			"(Suppress this warning with Control (%d) = 0.)\n",
			Info [UMFPACK_RCOND],
			1+UMFPACK_PIVOT_TOLERANCE,
			1+UMFPACK_SYM_PIVOT_TOLERANCE,
			1+UMFPACK_PRL) ;
		    mexWarnMsgTxt (warning) ;
		}
	    }

	}
	else if (do_det)
	{

	    /* -------------------------------------------------------------- */
	    /* get the determinant */
	    /* -------------------------------------------------------------- */

	    if (nargout == 2)
	    {
		/* [det dexp] = umfpack (A, 'det') ;
		 * return determinant in the form det * 10^dexp */
		p = &dexp ;
	    }
	    else
	    {
		/* [det] = umfpack (A, 'det') ;
		 * return determinant as a single scalar (overflow or
		 * underflow is much more likely) */
		p = (double *) NULL ;
	    }
	    if (A_is_complex)
	    {
		status = umfpack_zi_get_determinant (&dx, &dz, p,
			Numeric, Info) ;
		umfpack_zi_free_numeric (&Numeric) ;
	    }
	    else
	    {
		status = umfpack_di_get_determinant (&dx, p,
			Numeric, Info) ;
		umfpack_di_free_numeric (&Numeric) ;
		dz = 0 ;
	    }
	    if (status < 0)
	    {
		error ("extracting LU factors failed", A_is_complex, nargout,
		    pargout, Control, Info, status, do_info) ;
	    }
	    if (A_is_complex)
	    {
		pargout [0] = mxCreateDoubleMatrix (1, 1, mxCOMPLEX) ;
		p = mxGetPr (pargout [0]) ;
		*p = dx ;
		p = mxGetPi (pargout [0]) ;
		*p = dz ;
	    }
	    else
	    {
		pargout [0] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
		p = mxGetPr (pargout [0]) ;
		*p = dx ;
	    }
	    if (nargout == 2)
	    {
		pargout [1] = mxCreateDoubleMatrix (1, 1, mxREAL) ;
		p = mxGetPr (pargout [1]) ;
		*p = dexp ;
	    }

	}
	else
	{

	    /* -------------------------------------------------------------- */
	    /* get L, U, P, Q, and r */
	    /* -------------------------------------------------------------- */

	    if (A_is_complex)
	    {
	        status = umfpack_zi_get_lunz (&lnz, &unz, &ignore1, &ignore2,
		    &ignore3, Numeric) ;
	    }
	    else
	    {
	        status = umfpack_di_get_lunz (&lnz, &unz, &ignore1, &ignore2,
		    &ignore3, Numeric) ;
	    }

	    if (status < 0)
	    {
		if (A_is_complex)
		{
		    umfpack_zi_free_numeric (&Numeric) ;
		}
		else
		{
		    umfpack_di_free_numeric (&Numeric) ;
		}
		error ("extracting LU factors failed", A_is_complex, nargout,
		    pargout, Control, Info, status, do_info) ;
		return ;
	    }

	    /* avoid malloc of zero-sized arrays */
	    lnz = MAX (lnz, 1) ;
	    unz = MAX (unz, 1) ;

	    /* get temporary space, for the *** ROW *** form of L */
	    Ltp = (int *) mxMalloc ((n_row+1) * sizeof (int)) ;
	    Ltj = (int *) mxMalloc (lnz * sizeof (int)) ;
	    Ltx = (double *) mxMalloc (lnz * sizeof (double)) ;
	    if (A_is_complex)
	    {
	        Ltz = (double *) mxMalloc (lnz * sizeof (double)) ;
	    }
	    else
	    {
	        Ltz = (double *) NULL ;
	    }

	    /* create permanent copy of the output matrix U */
	    pargout [1] = mxCreateSparse (n_inner, n_col, unz, Atype) ;
	    Up = mxGetJc (pargout [1]) ;
	    Ui = mxGetIr (pargout [1]) ;
	    Ux = mxGetPr (pargout [1]) ;
	    Uz = mxGetPi (pargout [1]) ;

	    /* temporary space for the integer permutation vectors */
	    P = (int *) mxMalloc (n_row * sizeof (int)) ;
	    Q = (int *) mxMalloc (n_col * sizeof (int)) ;

	    /* get scale factors, if requested */
	    status2 = UMFPACK_OK ;
	    if (!no_scale)
	    {
		/* create a diagonal sparse matrix for the scale factors */
		pargout [4] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ;
		Rp = mxGetJc (pargout [4]) ;
		Ri = mxGetIr (pargout [4]) ;
		for (i = 0 ; i < n_row ; i++)
		{
		    Rp [i] = i ;
		    Ri [i] = i ;
		}
		Rp [n_row] = n_row ;
		Rs = mxGetPr (pargout [4]) ;
	    }
	    else
	    {
		Rs = (double *) NULL ;
	    }

	    /* get Lt, U, P, Q, and Rs from the numeric object */
	    if (A_is_complex)
	    {
		status = umfpack_zi_get_numeric (Ltp, Ltj, Ltx, Ltz, Up, Ui, Ux,
		    Uz, P, Q, (double *) NULL, (double *) NULL,
		    &do_recip, Rs, Numeric) ;
		umfpack_zi_free_numeric (&Numeric) ;
	    }
	    else
	    {
		status = umfpack_di_get_numeric (Ltp, Ltj, Ltx, Up, Ui,
		    Ux, P, Q, (double *) NULL,
		    &do_recip, Rs, Numeric) ;
		umfpack_di_free_numeric (&Numeric) ;
	    }

	    /* for the mexFunction, -DNRECIPROCAL must be set,
	     * so do_recip must be FALSE */

	    if (status < 0 || status2 < 0 || do_recip)
	    {
		mxFree (Ltp) ;
		mxFree (Ltj) ;
		mxFree (Ltx) ;
		if (Ltz) mxFree (Ltz) ;
		mxFree (P) ;
		mxFree (Q) ;
		mxDestroyArray (pargout [1]) ;
		error ("extracting LU factors failed", A_is_complex, nargout,
		    pargout, Control, Info, status, do_info) ;
		return ;
	    }

	    /* create sparse permutation matrix for P */
	    pargout [2] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ;
	    Pp = mxGetJc (pargout [2]) ;
	    Pi = mxGetIr (pargout [2]) ;
	    Px = mxGetPr (pargout [2]) ;
	    for (k = 0 ; k < n_row ; k++)
	    {
		Pp [k] = k ;
		Px [k] = 1 ;
		Pi [P [k]] = k ;
	    }
	    Pp [n_row] = n_row ;

	    /* create sparse permutation matrix for Q */
	    pargout [3] = mxCreateSparse (n_col, n_col, n_col, mxREAL) ;
	    Qp = mxGetJc (pargout [3]) ;
	    Qi = mxGetIr (pargout [3]) ;
	    Qx = mxGetPr (pargout [3]) ;
	    for (k = 0 ; k < n_col ; k++)
	    {
		Qp [k] = k ;
		Qx [k] = 1 ;
		Qi [k] = Q [k] ;
	    }
	    Qp [n_col] = n_col ;

	    /* permanent copy of L */
	    pargout [0] = mxCreateSparse (n_row, n_inner, lnz, Atype) ;
	    Lp = mxGetJc (pargout [0]) ;
	    Li = mxGetIr (pargout [0]) ;
	    Lx = mxGetPr (pargout [0]) ;
	    Lz = mxGetPi (pargout [0]) ;

	    /* convert L from row form to column form */
	    if (A_is_complex)
	    {
		/* non-conjugate array transpose */
	        status = umfpack_zi_transpose (n_inner, n_row, Ltp, Ltj, Ltx,
		    Ltz, (int *) NULL, (int *) NULL, Lp, Li, Lx, Lz, FALSE) ;
	    }
	    else
	    {
	        status = umfpack_di_transpose (n_inner, n_row, Ltp, Ltj, Ltx,
		    (int *) NULL, (int *) NULL, Lp, Li, Lx) ;
	    }

	    mxFree (Ltp) ;
	    mxFree (Ltj) ;
	    mxFree (Ltx) ;
	    if (Ltz) mxFree (Ltz) ;

	    if (status < 0)
	    {
		mxFree (P) ;
		mxFree (Q) ;
		mxDestroyArray (pargout [0]) ;
		mxDestroyArray (pargout [1]) ;
		mxDestroyArray (pargout [2]) ;
		mxDestroyArray (pargout [3]) ;
		error ("constructing L failed", A_is_complex, nargout, pargout,
		    Control, Info, status, do_info) ;
		return ;
	    }

	    /* -------------------------------------------------------------- */
	    /* print L, U, P, and Q */
	    /* -------------------------------------------------------------- */

	    if (A_is_complex)
	    {
		if (print_level >= 3) mexPrintf ("\nL: ") ;
	        (void) umfpack_zi_report_matrix (n_row, n_inner, Lp, Li,
		    Lx, Lz, 1, Control) ;
		if (print_level >= 3) mexPrintf ("\nU: ") ;
	        (void) umfpack_zi_report_matrix (n_inner, n_col,  Up, Ui,
		    Ux, Uz, 1, Control) ;
		if (print_level >= 3) mexPrintf ("\nP: ") ;
	        (void) umfpack_zi_report_perm (n_row, P, Control) ;
		if (print_level >= 3) mexPrintf ("\nQ: ") ;
	        (void) umfpack_zi_report_perm (n_col, Q, Control) ;
	    }
	    else
	    {
		if (print_level >= 3) mexPrintf ("\nL: ") ;
	        (void) umfpack_di_report_matrix (n_row, n_inner, Lp, Li,
		    Lx, 1, Control) ;
		if (print_level >= 3) mexPrintf ("\nU: ") ;
	        (void) umfpack_di_report_matrix (n_inner, n_col,  Up, Ui,
		    Ux, 1, Control) ;
		if (print_level >= 3) mexPrintf ("\nP: ") ;
	        (void) umfpack_di_report_perm (n_row, P, Control) ;
		if (print_level >= 3) mexPrintf ("\nQ: ") ;
	        (void) umfpack_di_report_perm (n_col, Q, Control) ;
	    }

	    mxFree (P) ;
	    mxFree (Q) ;

	}

    }
    else
    {

	/* ------------------------------------------------------------------ */
	/* return the symbolic factorization */
	/* ------------------------------------------------------------------ */

	Q = (int *) mxMalloc (n_col * sizeof (int)) ;
	P = (int *) mxMalloc (n_row * sizeof (int)) ;
	Front_npivcol = (int *) mxMalloc ((nn+1) * sizeof (int)) ;
	Front_parent = (int *) mxMalloc ((nn+1) * sizeof (int)) ;
	Front_1strow = (int *) mxMalloc ((nn+1) * sizeof (int)) ;
	Front_leftmostdesc = (int *) mxMalloc ((nn+1) * sizeof (int)) ;
	Chain_start = (int *) mxMalloc ((nn+1) * sizeof (int)) ;
	Chain_maxrows = (int *) mxMalloc ((nn+1) * sizeof (int)) ;
	Chain_maxcols = (int *) mxMalloc ((nn+1) * sizeof (int)) ;

	if (A_is_complex)
	{
	    status = umfpack_zi_get_symbolic (&ignore1, &ignore2, &ignore3,
	        &nz, &nfronts, &nchains, P, Q, Front_npivcol,
	        Front_parent, Front_1strow, Front_leftmostdesc,
	        Chain_start, Chain_maxrows, Chain_maxcols, Symbolic) ;
	    umfpack_zi_free_symbolic (&Symbolic) ;
	}
	else
	{
	    status = umfpack_di_get_symbolic (&ignore1, &ignore2, &ignore3,
	        &nz, &nfronts, &nchains, P, Q, Front_npivcol,
	        Front_parent, Front_1strow, Front_leftmostdesc,
	        Chain_start, Chain_maxrows, Chain_maxcols, Symbolic) ;
	    umfpack_di_free_symbolic (&Symbolic) ;
	}

	if (status < 0)
	{
	    mxFree (P) ;
	    mxFree (Q) ;
	    mxFree (Front_npivcol) ;
	    mxFree (Front_parent) ;
	    mxFree (Front_1strow) ;
	    mxFree (Front_leftmostdesc) ;
	    mxFree (Chain_start) ;
	    mxFree (Chain_maxrows) ;
	    mxFree (Chain_maxcols) ;
	    error ("extracting symbolic factors failed", A_is_complex, nargout,
		pargout, Control, Info, status, do_info) ;
	    return ;
	}

	/* create sparse permutation matrix for P */
	pargout [0] = mxCreateSparse (n_row, n_row, n_row, mxREAL) ;
	Pp = mxGetJc (pargout [0]) ;
	Pi = mxGetIr (pargout [0]) ;
	Px = mxGetPr (pargout [0]) ;
	for (k = 0 ; k < n_row ; k++)
	{
	    Pp [k] = k ;
	    Px [k] = 1 ;
	    Pi [P [k]] = k ;
	}
	Pp [n_row] = n_row ;

	/* create sparse permutation matrix for Q */
	pargout [1] = mxCreateSparse (n_col, n_col, n_col, mxREAL) ;
	Qp = mxGetJc (pargout [1]) ;
	Qi = mxGetIr (pargout [1]) ;
	Qx = mxGetPr (pargout [1]) ;
	for (k = 0 ; k < n_col ; k++)
	{
	    Qp [k] = k ;
	    Qx [k] = 1 ;
	    Qi [k] = Q [k] ;
	}
	Qp [n_col] = n_col ;

	/* create Fr */
	pargout [2] = mxCreateDoubleMatrix (nfronts+1, 4, mxREAL) ;

	p1 = mxGetPr (pargout [2]) ;
	p2 = p1 + nfronts + 1 ;
	p3 = p2 + nfronts + 1 ;
	p4 = p3 + nfronts + 1 ;
	for (i = 0 ; i <= nfronts ; i++)
	{
	    /* convert parent, 1strow, and leftmostdesc to 1-based */
	    p1 [i] = (double) (Front_npivcol [i]) ;
	    p2 [i] = (double) (Front_parent [i] + 1) ;
	    p3 [i] = (double) (Front_1strow [i] + 1) ;
	    p4 [i] = (double) (Front_leftmostdesc [i] + 1) ;
	}

	/* create Ch */
	pargout [3] = mxCreateDoubleMatrix (nchains+1, 3, mxREAL) ;
	p1 = mxGetPr (pargout [3]) ;
	p2 = p1 + nchains + 1 ;
	p3 = p2 + nchains + 1 ;
	for (i = 0 ; i < nchains ; i++)
	{
	    p1 [i] = (double) (Chain_start [i] + 1) ;	/* convert to 1-based */
	    p2 [i] = (double) (Chain_maxrows [i]) ;
	    p3 [i] = (double) (Chain_maxcols [i]) ;
	}
	p1 [nchains] = Chain_start [nchains] + 1 ;
	p2 [nchains] = 0 ;
	p3 [nchains] = 0 ;

	mxFree (P) ;
	mxFree (Q) ;
	mxFree (Front_npivcol) ;
	mxFree (Front_parent) ;
	mxFree (Front_1strow) ;
	mxFree (Front_leftmostdesc) ;
	mxFree (Chain_start) ;
	mxFree (Chain_maxrows) ;
	mxFree (Chain_maxcols) ;
    }

    /* ---------------------------------------------------------------------- */
    /* report Info */
    /* ---------------------------------------------------------------------- */

    if (A_is_complex)
    {
	umfpack_zi_report_info (Control, Info) ;
    }
    else
    {
	umfpack_di_report_info (Control, Info) ;
    }

    if (do_info > 0)
    {
	/* return Info */
	pargout [do_info] = mxCreateDoubleMatrix (1, UMFPACK_INFO, mxREAL) ;
	OutInfo = mxGetPr (pargout [do_info]) ;
	for (i = 0 ; i < UMFPACK_INFO ; i++)
	{
	    OutInfo [i] = Info [i] ;
	}
    }
}