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

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

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

/* Usage: [perm, invp] = ordmmd(P) where diag(P) = 0 */

#include "mex.h"

void convertin(int neqns, int nzmax, int *mxadj, 
			   int *madjncy, int *xadj, int *adjncy);
void  int2real(int neqns, int *perm, int *invp, 
			   double *mperm, double *minvp);
/* ------------------------------------------------------------
   GETADJ - Copies off-diagonal entries from C-style sparse
     matrix (cjc,cir) to Fortran style sparse matrix (forjc,forir).
     On input, n is number of columns.
   ------------------------------------------------------------ */

void mexFunction(const int nlhs, mxArray *plhs[], 
				 const int nrhs, const mxArray *prhs[])
{
/* Local variables */
    int flag, neqns, iwsiz, nofsub, nzmax;
    int *pmadjncy, *pxadj;
    int *pperm, *pinvp, *pmxadj, *piwork, *padjncy;
    double *pminvp, *pmperm;


/* CHECK FOR PROPER NUMBER OF ARGUMENTS */

    if (nrhs != 1) {
	mexErrMsgTxt("ORDMMD requires 1 input arguments");
    } else if (nlhs != 2) {
	mexErrMsgTxt("ORDMMD requires 2 output arguments");
    }

/* CHECK FOR DIMENSIONS OF INPUT ARGUMENTS */

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

/* DEREFERENCE INPUT ARGUMENTS TO GET ARRAY POINTERS */

    pmxadj = mxGetJc(prhs[0]);
    pmadjncy = mxGetIr(prhs[0]);

/* SPECIFY THE DIMENSIION OF WORKING VECTORS */

    nzmax = mxGetNzmax(prhs[0]);
    iwsiz = neqns << 2;

/* CREATE WORKING PARAMETERS */

    pxadj = (int *)mxCalloc(neqns + 1, sizeof(int));
    padjncy = (int *)mxCalloc(nzmax, sizeof(int));
    pperm =(int *)mxCalloc(neqns,sizeof(int));
    pinvp = (int *)mxCalloc(neqns, sizeof(int));
    piwork = (int *)mxCalloc(iwsiz, sizeof(int));

/* INPUT DATA CONVERSION */
	convertin(neqns, nzmax, pmxadj, pmadjncy, pxadj, padjncy);

/* DO THE ACTUAL COMPUTATIONS IN A FORTRAN SUBROUTINE */
	ordmmd_(&neqns, pxadj, padjncy, pinvp, pperm, &iwsiz, piwork, &nofsub, &flag);


/* CHECK ERROR FLAG */

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

/* CREATE MATRICES FOR RETURN ARGUMENTS */

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

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

    pmperm = mxGetPr(plhs[0]);
    pminvp = mxGetPr(plhs[1]);
	int2real(neqns, pperm, pinvp, pmperm, pminvp);

/* Release working arrays */
	mxFree(piwork);
	mxFree(pinvp);
	mxFree(pperm);
	mxFree(pxadj);
	mxFree(padjncy);

} /* mexFunction */

/* ------------------------------------------------------ */
/* Convert mxadj and madjncy of range [0:n-1] */
/*        to  xadj and  adjncy of range [1:n]; */
/* ------------------------------------------------------ */
void convertin(int neqns, int nzmax, int *mxadj, 
			   int *madjncy, int *xadj, int *adjncy)
{
    int i;

    for (i = 0; i < neqns + 1; ++i) {
	xadj[i] = mxadj[i] + 1;
    }
    for (i = 0; i < nzmax; ++i) {
	adjncy[i] = madjncy[i] + 1;
    }
} /* convertin */

/* ------------------------------------------------- */
/* Convert outputs from integer type in C */
/*                   to real*8  type in Matlab */
/* ------------------------------------------------- */
void  int2real(int neqns, int *perm, int *invp, 
			   double *mperm, double *minvp)
{
    int i;

    for (i = 0; i < neqns; ++i) {
	mperm[i] = (double)perm[i];
	minvp[i] = (double)invp[i];
    }
} /* int2real */

