/* *********************************************************************** */

/*   Version:        Beta-2.2 */
/*   Last modified:  January 13, 1995 */
/*   Author:         Yin Zhang */
/*                   Department of Mathematics and Statistics */
/*                   University of Maryland Baltimore County */

/* *********************************************************************** */

/* Usage: [sol] = blkslv(xlnz,xsuper,xlindx,lindx,lnz,rrhs) */

#if !defined(MAX)
#define  MAX(A, B)   ((A) > (B) ? (A) : (B))
#endif

#include "mex.h"

void convertin(int neqns, int maxsub, int nsuper, double *mxlnz, double *mxsuper, 
	       int *xlnz, int *xsuper, double *mxlindx, double *mlindx, int *xlindx, 
	       int *lindx);

void copysol(int neqns, double *rhs, double *sol);


void mexFunction(const int nlhs, mxArray *plhs[],
		 const int nrhs, const mxArray *prhs[])
{
    /* System generated locals */
    int i__1, i__2;

    /* Local variables */
    int nsuper, neqns, maxsub;
    int *pxlnz, *pxlindx, *plindx, *pxsuper;
    double *pmxlnz, *pmxsuper, *pmxlindx, *pmlindx, *prrhs, *plnz, *psol;


/* KEEP THE ABOVE SUBROUTINE, ARGUMENT, AND FUNCTION DECLARATIONS */
/* FOR USE IN ALL YOUR FORTRAN MEX FILES. */
/* --------------------------------------------------------------------- 
*/


/* CHECK FOR PROPER NUMBER OF ARGUMENTS */

    if (nrhs != 6) {
	mexErrMsgTxt("SUPSLV requires 6 input arguments");
    } else if (nlhs != 1) {
	mexErrMsgTxt("SUPSLV requires 1 output argument");
    }

/* SPECIFY THE DIMENSIION OF WORKING VECTORS */

/* Computing MAX */
    i__1 = mxGetM(prhs[5]), i__2 = mxGetN(prhs[5]);
    neqns = MAX(i__1,i__2);
/* Computing MAX */
    i__1 = mxGetM(prhs[1]), i__2 = mxGetN(prhs[1]);
    nsuper = MAX(i__1,i__2) - 1;
/* Computing MAX */
    i__1 = mxGetM(prhs[3]), i__2 = mxGetN(prhs[3]);
    maxsub = MAX(i__1,i__2);

/* DEREFERENCE ARGUMENTS TO GET ARRAY POINTERS */

    pmxlnz = mxGetPr(prhs[0]);
    pmxsuper = mxGetPr(prhs[1]);
    pmxlindx = mxGetPr(prhs[2]);
    pmlindx = mxGetPr(prhs[3]);
    plnz = mxGetPr(prhs[4]);
    prrhs = mxGetPr(prhs[5]);

/* CREATE WORKING PARAMETERS */

    pxlnz = (int *)mxCalloc(neqns + 1, sizeof(int));
    pxsuper = (int *)mxCalloc(nsuper + 1, sizeof(int));
    pxlindx = (int *)mxCalloc(nsuper + 1, sizeof(int));
    plindx = (int *)mxCalloc(maxsub, sizeof(int));

/* INPUT DATA TRANSFORMATION */
	convertin(neqns, maxsub, nsuper, pmxlnz, pmxsuper, pxlnz, pxsuper, 
		  pmxlindx, pmlindx, pxlindx, plindx);

/* CREATE OUTPUT PARAMETER */

    plhs[0] = mxCreateDoubleMatrix(neqns, 1, mxREAL);
    psol = mxGetPr(plhs[0]);

/* DO THE ACTUAL COMPUTATIONS IN A FORTRAN SUBROUTINE */
	blkslv_(&nsuper, pxsuper, pxlindx, plindx, pxlnz, plnz, prrhs);

/* Copy solution to output vector */
	copysol(neqns, prrhs, psol);

/* Release working arrays */
	mxFree(pxlnz);
	mxFree(pxsuper);
	mxFree(pxlindx);
	mxFree(plindx);

} /* mexFunction*/

/* ---------------------------------------------------------------- */
/* Convert from type real*8 in Matlab to type integer in Fortran */
/* ---------------------------------------------------------------- */
void convertin(int neqns, int maxsub, int nsuper, double *mxlnz, double *mxsuper, 
	       int *xlnz, int *xsuper, double *mxlindx, double *mlindx, int *xlindx, 
	       int *lindx)
{
    int i;

    for (i = 0; i < neqns + 1; ++i) {
	xlnz[i] = (int) mxlnz[i];
/* L10: */
    }
    for (i = 0; i < nsuper + 1; ++i) {
	xsuper[i] = (int) mxsuper[i];
	xlindx[i] = (int) mxlindx[i];
/* L20: */
    }
    for (i = 0; i < maxsub; ++i) {
	lindx[i] = (int) mlindx[i];
/* L30: */
    }
} /* convertin */

/* ---------------------------------------------------------------- */
/* Copy solution to output vector */
/* ---------------------------------------------------------------- */
void copysol(int neqns, double *rhs, double *sol)
{
    int i;

    for (i = 0; i < neqns; ++i) {
	sol[i] = rhs[i];
/* L10: */
    }
} /* copysol_*/

