/*
 * -- SuperLU routine (version 3.0) --
 * Univ. of California Berkeley, Xerox Palo Alto Research Center,
 * and Lawrence Berkeley National Lab.
 * October 15, 2003
 *
 */
#include <stdio.h>
#include "mex.h"
#include "slu_ddefs.h"

#ifdef V5
#define  MatlabMatrix mxArray
#else    /* V4 */
#define  MatlabMatrix Matrix
#endif


/* Aliases for input and output arguments */
#define A_in		prhs[0]
#define b_in    	prhs[1]
#define Pc_in		prhs[2]
#define x_out		plhs[0]

#define verbose (SPUMONI>0)
#define babble  (SPUMONI>1)
#define burble  (SPUMONI>2)

void mexFunction(
    int          nlhs,           /* number of expected outputs */
    MatlabMatrix *plhs[],        /* matrix pointer array returning outputs */
    int          nrhs,           /* number of inputs */
#ifdef V5
    const MatlabMatrix *prhs[]   /* matrix pointer array for inputs. */
#else /* V4 */
    MatlabMatrix *prhs[]         /* matrix pointer array for inputs */
#endif
    )
{
    int SPUMONI;             /* ... as should the sparse monitor flag */
#ifdef V5
    double FlopsInSuperLU;   /* ... as should the flop counter. */
#else /* V4 */
    Real FlopsInSuperLU;     /* ... as should the flop counter */
#endif
    extern flops_t LUFactFlops(SuperLUStat_t*);
    extern flops_t LUSolveFlops(SuperLUStat_t*);
    
    /* Arguments to dgssv(). */
    SuperMatrix A;
    SuperMatrix B;
    SuperMatrix L, U;
    int	   	m, n, nnz;
    int         numrhs;
    double 	*vb, *x;
    double      *val;
    int       	*rowind;
    int		*colptr;
    int    	*perm_r, *perm_c;
    int		info;
    MatlabMatrix *X, *Y;            /* args to calls back to Matlab */
    int         i, mexerr;
    superlu_options_t options;
    SuperLUStat_t stat;

    /* Check number of arguments passed from Matlab. */
    if (nrhs != 3) {
	mexErrMsgTxt("LUSOLVE requires 3 input arguments.");
    } else if (nlhs != 1) {
      	mexErrMsgTxt("LUSOLVE requires 1 output argument.");
    }   

    /* Read the Sparse Monitor Flag */
    X = mxCreateString("spumoni");
    mexerr = mexCallMATLAB(1, &Y, 1, &X, "sparsfun");
    SPUMONI = mxGetScalar(Y);
#ifdef V5
    mxDestroyArray(Y);
    mxDestroyArray(X);
#else
    mxFreeMatrix(Y);
    mxFreeMatrix(X);
#endif

    m = mxGetM(A_in);
    n = mxGetN(A_in);
    numrhs = mxGetN(b_in);
    if ( babble ) printf("m=%d, n=%d, numrhs=%d\n", m, n, numrhs);
    vb = mxGetPr(b_in);
#ifdef V5
    x_out = mxCreateDoubleMatrix(m, numrhs, mxREAL);
#else
    x_out = mxCreateFull(m, numrhs, REAL);
#endif
    x = mxGetPr(x_out);
    perm_r = (int *) mxCalloc(m, sizeof(int));
    perm_c = mxGetIr(Pc_in); 

    val = mxGetPr(A_in);
    rowind = mxGetIr(A_in);
    colptr = mxGetJc(A_in);
    nnz = colptr[n];
    dCreate_CompCol_Matrix(&A, m, n, nnz, val, rowind, colptr,
			   SLU_NC, SLU_D, SLU_GE);
    dCopy_Dense_Matrix(m, numrhs, vb, m, x, m);
    dCreate_Dense_Matrix(&B, m, numrhs, x, m, SLU_DN, SLU_D, SLU_GE);

    FlopsInSuperLU = 0;

    set_default_options(&options);
    options.ColPerm = MY_PERMC;
    StatInit(&stat);

    /* Call simple driver */
    if ( verbose )
      mexPrintf("Call LUSOLVE, use SUPERLU to factor first ...\n");
    dgssv(&options, &A, perm_c, perm_r, &L, &U, &B, &stat, &info);

#if 0 /* FLOPS is not available in the new Matlab. */
    /* Tell Matlab how many flops we did. */
    FlopsInSuperLU += LUFactFlops(&stat) + LUSolveFlops(&stat);
    if ( verbose ) mexPrintf("LUSOLVE flops: %.f\n", FlopsInSuperLU);
    mexerr = mexCallMATLAB(1, &X, 0, NULL, "flops");
    *(mxGetPr(X)) += FlopsInSuperLU;
    mexerr = mexCallMATLAB(1, &Y, 1, &X, "flops");
#ifdef V5
    mxDestroyArray(Y);
    mxDestroyArray(X);
#else
    mxFreeMatrix(Y);
    mxFreeMatrix(X);
#endif
#endif

    /* Construct Matlab solution matrix. */
    if ( !info ) {
        Destroy_SuperNode_Matrix(&L);
        Destroy_CompCol_Matrix(&U);
        if ( babble ) printf("Destroy L & U from SuperLU...\n");
    } else {
	mexErrMsgTxt("Error returned from C dgssv().");
    }

    mxFree(perm_r);
    StatFree(&stat);

    return;
 
}