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

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

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

/* Usage: [lnz] = blkfct(xlnz,xsuper,snode,split,xlindx,lindx,... */
/*                       lnz,tmpsiz,level) */

#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, 
				double *msnode, double *msplit, double *mxlindx, double *mlindx, 
				int *xlnz, int *xsuper, int *snode, int *split, int *xlindx, int *lindx);

extern int mmpy1_(int *m, int *n, int *q, int *xpnt,double *x, double *y, int *ldy);
extern int mmpy2_(int *m, int *n, int *q, int *xpnt,double *x, double *y, int *ldy);
extern int mmpy4_(int *m, int *n, int *q, int *xpnt,double *x, double *y, int *ldy);
extern int mmpy8_(int *m, int *n, int *q, int *xpnt,double *x, double *y, int *ldy);
extern int smxpy1_(int *m, int *n, double *y, int *apnt, double *a);
extern int smxpy2_(int *m, int *n, double *y, int *apnt, double *a);
extern int smxpy4_(int *m, int *n, double *y, int *apnt, double *a);
extern int smxpy8_(int *m, int *n, double *y, int *apnt, double *a);

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

    /* Local variables */
    int flag, neqns, iwsiz, nsuper, tmpsiz, maxsub, level;
    int *pxlnz, *psnode, *plindx, *psplit, *piwork, *pxlindx, *pxsuper;
    double *ptmpvec, *tmp;
    double *pmxlnz, *pmxlindx, *pmxsuper, *pmsnode, *pmlindx, *plnz, *pmsplit;

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


/* CHECK FOR PROPER NUMBER OF ARGUMENTS */

    /* Function Body */
    if (nrhs != 9) {
	mexErrMsgTxt("BLKFCT requires 9 input arguments");
    } else if (nlhs != 1) {
	mexErrMsgTxt("BLKFCT requires 1 output arguments");
    }

/* SPECIFY THE DIMENSIION OF WORKING VECTORS */

/* 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[2]), i__2 = mxGetN(prhs[2]);
    neqns = MAX(i__1,i__2);
/* Computing MAX */
    i__1 = mxGetM(prhs[5]), i__2 = mxGetN(prhs[5]);
    maxsub = MAX(i__1,i__2);
    tmpsiz = (int) mxGetScalar(prhs[7]);
    level = (int) mxGetScalar(prhs[8]);
    iwsiz = (neqns << 1) + (nsuper << 1);

/* DEREFERENCE ARGUMENTS TO GET ARRAY POINTERS */

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

/* CREATE WORKING PARAMETERS */

    pxlnz = (int *)mxCalloc(neqns + 1, sizeof(int));
    pxsuper = (int *)mxCalloc(nsuper + 1, sizeof(int));
    psnode = (int *)mxCalloc(neqns, sizeof(int));
    psplit = (int *)mxCalloc(neqns, sizeof(int));
    pxlindx = (int *)mxCalloc(nsuper + 1, sizeof(int));
    plindx = (int *)mxCalloc(maxsub, sizeof(int));
    piwork = (int *)mxCalloc(iwsiz, sizeof(int));
    ptmpvec = (double *)mxCalloc(tmpsiz, sizeof(double));

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

/* CREATE OUTPUT PARAMETERS */

    plhs[0] = prhs[6];

/* DO THE ACTUAL COMPUTATIONS IN A FORTRAN SUBROUTINE */

    if (level == 1) {
		blkfct_(&neqns, &nsuper, pxsuper, psnode, psplit, pxlindx, plindx, pxlnz,
               plnz, &iwsiz, piwork, &tmpsiz, ptmpvec, &flag, &mmpy1_, &smxpy1_);
    }
    if (level == 2) {
		blkfct_(&neqns, &nsuper, pxsuper, psnode, psplit, pxlindx, plindx, pxlnz,
               plnz, &iwsiz, piwork, &tmpsiz, ptmpvec, &flag, &mmpy2_, &smxpy2_);
    }
    if (level == 4) {
		blkfct_(&neqns, &nsuper, pxsuper, psnode, psplit, pxlindx, plindx, pxlnz,
               plnz, &iwsiz, piwork, &tmpsiz, ptmpvec, &flag, &mmpy4_, &smxpy4_);

    }
    if (level == 8) {
		blkfct_(&neqns, &nsuper, pxsuper, psnode, psplit, pxlindx, plindx, pxlnz,
               plnz, &iwsiz, piwork, &tmpsiz, ptmpvec, &flag, &mmpy8_, &smxpy8_);
    }

/* CHECK ERROR FLAG */

    if (flag < -1) {
	mexErrMsgTxt("Insufficient working storage in blkfct.");
    }

/* Release working arrays */
	mxFree(pxlnz);
	mxFree(pxsuper);
	mxFree(psnode);
	mxFree(psplit);
	mxFree(pxlindx);
	mxFree(plindx);
	mxFree(piwork);
	mxFree(ptmpvec);

} /* 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, 
				double *msnode, double *msplit, double *mxlindx, double *mlindx, 
				int *xlnz, int *xsuper, int *snode, int *split, int *xlindx, int *lindx)
{
    int i;

    for (i = 0; i < neqns; ++i) {
	xlnz[i] = (int) mxlnz[i];
	snode[i] = (int) msnode[i];
	split[i] = (int) msplit[i];
/* L10: */
    }
    xlnz[neqns] = (int) mxlnz[neqns];

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

