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

/* Usage: [xlnz,nnzl,xsuper,xlindx,lindx,snode,split,tmpsiz,perm,invp]... */
/*                                     = symfct(P, perm, invp, cachsz) */
#if !defined(MIN)
#define  MIN(A, B)   ((A) < (B) ? (A) : (B))
#endif
#if !defined(MAX)
#define  MAX(A, B)   ((A) > (B) ? (A) : (B))
#endif

#include "mex.h"

void convertin(int neqns, int anzmax, int *mxadj, int *madjncy, double *mperm, 
	       double *minvp, int *xadj, int *adjncy, int *perm, int *invp);

void int2real(int neqns, int nsuper, int nsub, int *xlnz, int nnzl, int *xsuper, 
	      int *xlindx, int *lindx, int *snode, int *split, int *perm, int *invp, 
	      int tmpsiz, double *mxlnz, double *mnnzl, double *mxsuper, double *mxlindx, 
	      double *mlindx, double *msnode, double *msplit, double *mperm, double *minvp, 
	      double *mtmpsiz);

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

void mexFunction(const int nlhs, mxArray *plhs[],
		 const int nrhs, const mxArray *prhs[])
{
/* Local variables */
    int flag, nsub, m, n, neqns, cachsz, iwsiz, anzmax, nsuper, tmpsiz, nnzl;
    int    *pxadj, *pxlindx, *pinvp, *pperm, *pxlnz, *padjncy;
    int    *psnode, *pcolcnt, *piwork, *psplit, *pxsuper, *plindx;
    int    *pmadjncy, *pmxadj;
    double *pminvp, *pmperm, *pmxlindx, *pmlindx, *pmtmpsiz;
    double *pmsplit, *pmnnzl, *pmxlnz, *pmxsuper, *pmsnode;


/* CHECK FOR PROPER NUMBER OF ARGUMENTS */

    if (nrhs != 4) {
	mexErrMsgTxt("SYMFCT requires 4 input arguments");
    } else if (nlhs != 10) {
	mexErrMsgTxt("SYMFCT requires 10 output arguments");
    }

/* CHECK FOR DIMENSIONS OF INPUT ARGUMENTS */

    neqns = mxGetM(prhs[0]);
    if (neqns != mxGetN(prhs[0])) {
	mexErrMsgTxt("Input matrix must be square");
    }

    m = mxGetM(prhs[1]);
    n = mxGetN(prhs[1]);
    if (MAX(m,n) != neqns || MIN(m,n) != 1) {
	mexErrMsgTxt("SYMFCT requires PERM to be neqns x 1");
    }

    m = mxGetM(prhs[2]);
    n = mxGetN(prhs[2]);
    if (MAX(m,n) != neqns || MIN(m,n) != 1) {
	mexErrMsgTxt("SYMFCT requires INVP to be neqns x 1");
    }

    m = mxGetM(prhs[3]);
    n = mxGetN(prhs[3]);
    if (m != 1 || n != 1) {
	mexErrMsgTxt("SYMFCT requires CACHSZ to be 1 x 1");
    }

/* DEREFERENCE INPUT ARGUMENTS TO GET ARRAY POINTERS */

    pmxadj = mxGetJc(prhs[0]);
    pmadjncy = mxGetIr(prhs[0]);
    pmperm = mxGetPr(prhs[1]);
    pminvp = mxGetPr(prhs[2]);

/* SPECIFY THE DIMENSIION OF WORKING VECTORS */

    anzmax = mxGetNzmax(prhs[0]);
    iwsiz = neqns * 7 + 3;

/* CREATE WORKING PARAMETERS */

    pxadj   = (int *)mxCalloc(neqns + 1, sizeof(int));
    padjncy = (int *)mxCalloc(anzmax, sizeof(int));
    pperm   = (int *)mxCalloc(neqns, sizeof(int));
    pinvp   = (int *)mxCalloc(neqns, sizeof(int));
    pxlnz   = (int *)mxCalloc(neqns + 1, sizeof(int));
    pxsuper = (int *)mxCalloc(neqns + 1, sizeof(int));
    pxlindx = (int *)mxCalloc(neqns + 1, sizeof(int));
    psnode  = (int *)mxCalloc(neqns, sizeof(int));
    psplit  = (int *)mxCalloc(neqns, sizeof(int));
    pcolcnt = (int *)mxCalloc(neqns, sizeof(int));
    piwork  = (int *)mxCalloc(iwsiz, sizeof(int));

/* INPUT DATA CONVERSION */
     convertin(neqns, anzmax, pmxadj, pmadjncy, pmperm, pminvp, pxadj,
	       padjncy, pperm, pinvp);

/* CALL THE ACTUAL FORTRAN SUBROUTINES */
     sfinit_(&neqns, anzmax, pxadj, padjncy, pperm, pinvp, pcolcnt, &nnzl, &nsub, &nsuper,
	      psnode, pxsuper, &iwsiz, piwork, &flag);

/* CHECK ERROR FLAG */

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

    plindx = (int *)mxCalloc(nsub << 1, sizeof(int));
    symfct_(&neqns,  anzmax, pxadj, padjncy, pperm, pinvp, pcolcnt, &nsuper,
	     pxsuper, psnode, &nsub, pxlindx, plindx, pxlnz, &iwsiz, piwork, &flag);

/* CHECK ERROR FLAG */

    if (flag == -1) {
	mexErrMsgTxt("Insufficient integer working space in symfct.");
    }
    if (flag == -2) {
	mexErrMsgTxt("Inconsistancy in the input in symfct.");
    }

    cachsz = (int) mxGetScalar(prhs[3]);
    bfinit_(&neqns, &nsuper, pxsuper, psnode, pxlindx, plindx, &cachsz, &tmpsiz, psplit);

/* CREATE MATRICES FOR RETURN ARGUMENTS */

    plhs[0] = mxCreateDoubleMatrix(neqns + 1, 1, mxREAL);
    plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
    plhs[2] = mxCreateDoubleMatrix(nsuper + 1, 1, mxREAL);
    plhs[3] = mxCreateDoubleMatrix(nsuper + 1, 1, mxREAL);
    plhs[4] = mxCreateDoubleMatrix(nsub, 1, mxREAL);
    plhs[5] = mxCreateDoubleMatrix(neqns, 1, mxREAL);
    plhs[6] = mxCreateDoubleMatrix(neqns, 1, mxREAL);
    plhs[7] = mxCreateDoubleMatrix(1, 1, mxREAL);
    plhs[8] = prhs[1];
    plhs[9] = prhs[2];

/* DEREFERENCE OUTPUT ARGUMENTS TO GET REAL PART POINTERS */

    pmxlnz   = mxGetPr(plhs[0]);
    pmnnzl   = mxGetPr(plhs[1]);
    pmxsuper = mxGetPr(plhs[2]);
    pmxlindx = mxGetPr(plhs[3]);
    pmlindx  = mxGetPr(plhs[4]);
    pmsnode  = mxGetPr(plhs[5]);
    pmsplit  = mxGetPr(plhs[6]);
    pmtmpsiz = mxGetPr(plhs[7]);

    int2real(neqns, nsuper, nsub, pxlnz, nnzl, pxsuper, pxlindx, plindx, psnode, psplit,
             pperm, pinvp, tmpsiz, pmxlnz, pmnnzl, pmxsuper, pmxlindx, pmlindx, pmsnode, 
	     pmsplit, pmperm, pminvp, pmtmpsiz);

/* Release working arrays*/
	mxFree(piwork);
	mxFree(pinvp);
	mxFree(pperm);
	mxFree(pxadj);
	mxFree(padjncy);
	mxFree(psplit);
	mxFree(pxlindx);
	mxFree(pxsuper);
	mxFree(psnode);
	mxFree(pcolcnt);
} /* mexFunction */

/* ------------------------------------------------------ */
/* Convert mxadj and madjncy of range [0:n-1] */
/*        to  xadj and  adjncy of range [1:n]; */
/* Convert mperm and minvp of type real*8 in Matlab */
/*        to  perm and  invp of type integer in Fortran */
/* ------------------------------------------------------ */
void convertin(int neqns, int anzmax, int *mxadj, int *madjncy, double *mperm, 
	       double *minvp, int *xadj, int *adjncy, int *perm, int *invp)
{
    /* Local variables */
    int i;

    for (i = 0; i < neqns; ++i) {
	xadj[i] = mxadj[i] + 1;
	perm[i] = (int) mperm[i];
	invp[i] = (int) minvp[i];
/* L10: */
    }
    xadj[neqns] = mxadj[neqns] + 1;

    for (i = 0; i < anzmax; ++i) {
	adjncy[i] = madjncy[i] + 1;
/* L20: */
    }
} /* convertin */

/* ------------------------------------------------- */
/* Convert outputs from integer type in Fortran */
/*                   to real*8  type in Matlab */
/* ------------------------------------------------- */
void int2real(int neqns, int nsuper, int nsub, int *xlnz, int nnzl, int *xsuper, 
	      int *xlindx, int *lindx, int *snode, int *split, int *perm, int *invp, 
	      int tmpsiz, double *mxlnz, double *mnnzl, double *mxsuper, double *mxlindx, 
	      double *mlindx, double *msnode, double *msplit, double *mperm, double *minvp, 
	      double *mtmpsiz)
{

    int i;

    *mnnzl = (double) (nnzl);
    *mtmpsiz = (double) (tmpsiz);

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

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

