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

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

/* *********************************************************************** */
/* Usage: [lnz] = inpnv(Pdiag,P,invp,perm,xlnz,xsuper,xlindx,lindx,nnzl) */

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

#include "mex.h"

void convertin(int neqns, int maxsub, int nzmax, int nsuper, int *mip, int *mjp, 
			   double *minvp, double *mperm, int *ip, int *jp, int *invp, int *perm, 
			   double *mxlnz, double *mxsuper, double *mxlindx, double *mlindx, 
			   int *xlnz, int *xsuper, int *xlindx, int *lindx);
void addiag(int n, int *jp, int *ip, double *p, double *pdiag, int *xadjf, int *adjf, 
			double *anzf);

void mexFunction(const int nlhs, mxArray *plhs[],
				 const int nrhs, const mxArray *prhs[])
{
/* Local variables */
    int i__1, i__2;
    int nnzl, neqns, maxsub, nzmax, nsuper;
	int *pmip, *pmjp, *padjf, *plink,  *pperm, *pinvp, *pxlnz; 
    int *pip, *pjp, *pxlindx, *pxadjf, *plindx, *pxsuper;
	double *panzf, *ppdiag, *pp, *pminvp, *pmperm, *pmxlnz;
	double *pmxsuper, *pmlindx, *plnz, *pmxlindx;

/* CHECK FOR PROPER NUMBER OF ARGUMENTS */

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

/* SPECIFY THE DIMENSIION OF WORKING VECTORS */

/* Computing MAX */
    i__1 = mxGetM(prhs[1]), i__2 = mxGetN(prhs[1]);
    neqns = MAX(i__1,i__2);
/* Computing MAX */
    i__1 = mxGetM(prhs[5]), i__2 = mxGetN(prhs[5]);
    nsuper = MAX(i__1,i__2) - 1;
/* Computing MAX */
    i__1 = mxGetM(prhs[7]), i__2 = mxGetN(prhs[7]);
    maxsub = MAX(i__1,i__2);
    nnzl = (int) mxGetScalar(prhs[8]);
    nzmax = mxGetNzmax(prhs[1]);

/* CREATE OUTPUT PARAMETERS */

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

/* CREATE WORKING PARAMETERS */

    pip = (int *)mxCalloc(nzmax, sizeof(int));
    pjp = (int *)mxCalloc(neqns + 1, sizeof(int));
    pinvp = (int *)mxCalloc(neqns, sizeof(int));
    pperm = (int *)mxCalloc(neqns, sizeof(int));
    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));
    pxadjf = (int *)mxCalloc(neqns + 1, sizeof(int));
    plink = (int *)mxCalloc(neqns, sizeof(int));
    padjf = (int *)mxCalloc(nzmax + neqns, sizeof(int));
    panzf = (double *)mxCalloc(nzmax + neqns, sizeof(double));


/* DEREFERENCE ARGUMENTS TO GET ARRAY POINTERS */

    ppdiag = mxGetPr(prhs[0]);
    pp = mxGetPr(prhs[1]);
    pmip = mxGetIr(prhs[1]);
    pmjp = mxGetJc(prhs[1]);
    pminvp = mxGetPr(prhs[2]);
    pmperm = mxGetPr(prhs[3]);
    pmxlnz = mxGetPr(prhs[4]);
    pmxsuper = mxGetPr(prhs[5]);
    pmxlindx = mxGetPr(prhs[6]);
    pmlindx = mxGetPr(prhs[7]);
    plnz = mxGetPr(plhs[0]);

/* INPUT DATA TRANSFORMATION */
	convertin(neqns, maxsub, nzmax, nsuper, pmip,  pmjp, pminvp, pmperm,
              pip, pjp, pinvp, pperm, pmxlnz, pmxsuper, pmxlindx, pmlindx,
              pxlnz, pxsuper, pxlindx, plindx);
	
/* DO THE ACTUAL COMPUTATIONS IN A FORTRAN SUBROUTINE */
	addiag(neqns, pjp, pip, pp, ppdiag, pxadjf, padjf, panzf);
	
	inpnv_(&neqns, pxadjf, padjf, panzf, pperm, pinvp, &nsuper, pxsuper,
		       pxlindx, plindx, pxlnz, plnz, plink);

/* Release working arrays*/	
	mxFree(pip);
    mxFree(pjp);
    mxFree(pinvp);
    mxFree(pperm);
    mxFree(pxlnz);
    mxFree(pxsuper);
    mxFree(pxlindx);
    mxFree(plindx);
    mxFree(pxadjf);
    mxFree(plink);
    mxFree(padjf);
    mxFree(panzf);

} /* mexFunction */

/* ---------------------------------------------------------------- */
/* Convert from C range [0:n-1] to Fortran range [1:n]; */
/* Convert from type real*8 in Matlab to type integer in Fortran */
/* ---------------------------------------------------------------- */
void convertin(int neqns, int maxsub, int nzmax, int nsuper, int *mip, int *mjp, 
			   double *minvp, double *mperm, int *ip, int *jp, int *invp, int *perm, 
			   double *mxlnz, double *mxsuper, double *mxlindx, double *mlindx, 
			   int *xlnz, int *xsuper, int *xlindx, int *lindx)
{

    /* Local variables */
    int i;

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

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

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

/*   Last modified:  November, 1994 */
/*   Author:         Detong Zhang */
/*                   Department of Mathematics and Statistics */
/*                   University of Maryland Baltimore County */

/* *********************************************************************** */
void addiag(int n, int *jp, int *ip, double *p, double *pdiag, int *xadjf, int *adjf, 
			double *anzf)
{
    /* Local variables */
    int j, ij;

    for (j = 0; j < n; ++j) {
	anzf[jp[j] + j - 1] = pdiag[j];
	adjf[jp[j] + j - 1] = j + 1;
	for (ij = jp[j] - 1; ij < jp[j + 1] - 1; ++ij) {
	    anzf[ij + j + 1] = p[ij];
	    adjf[ij + j + 1] = ip[ij];
/* L50: */
	}
	xadjf[j] = jp[j] + j;
/* L40: */
    }
    xadjf[n] = jp[n] + n;
} /* addiag */

