#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include "mex.h" 
/*#include "../../include/mkl_blas.h"*/

extern int omp_get_max_threads();
/* PARDISO prototype. */
#if defined(_WIN32) || defined(_WIN64)
#define pardiso_ PARDISO
#define pardisoinit_ PARDISOINIT
#else
#define PARDISO pardiso_
#define PARDISOINIT pardisoinit_
#endif
extern int PARDISO
	(void *, int *, int *, int *, int *, int *,
	double *, int *, int *, int *, int *, int *,
	int *, double *, double *, int *);

extern int PARDISOINIT
        (void *, int *, int *);


static void ClearAllFactors(void) {
  mexPrintf("Cleanup of PARDISO factors not properly implemented\n");
}

/* Internal solver memory pointer pt, */
/* 32-bit: int pt[64]; 64-bit: long int pt[64] */
/* or void *pt[64] should be OK on both architectures */
/* static void *PardisoPt[64]; */
/*--------------------------------------------------------------------*/
typedef struct _FactorInfo FactorInfo ;
struct _FactorInfo {

  void         *PardisoPt[512];
  int          nRows;
  int          *rows_1;
  int          *columns_1;
  double       *pr;
} ;
/*--------------------------------------------------------------------*/
void ClearMtx(FactorInfo *FI, int cF) {
  mexPrintf("Clear factor\n");
}
/*--------------------------------------------------------------------*/
static FactorInfo     FI[20];
static int            cF=-1, indF[20];

/*--------------------------------------------------------------------*/
/*--------------------------------------------------------------------*/

void PardisoError(int error) {
 if (error==-2)   mexPrintf("error = %i (Out of memory)\n",error);
}

void ClearFactor(int cF, int maxfct, int mnum, int mtype, int* param,
   int msglvl) {

  int             phase, idum, nRhs=0;
  double          ddum;

 phase = -1; /* Release internal memory. */
 cF=0;
 if (msglvl) mexPrintf("Clearing factor %i",cF);
 if (FI[cF].rows_1 == NULL) mexErrMsgTxt("Problem with rows_1 pointers");
 if (FI[cF].columns_1 == NULL) mexErrMsgTxt("Problem with columns_1 pointers");
 if (FI[cF].pr == NULL) mexErrMsgTxt("Problem with pr pointers");
 PARDISO (FI[cF].PardisoPt, &param[100], &param[101], &mtype, &phase,
    &(FI[cF].nRows), FI[cF].pr, FI[cF].rows_1, FI[cF].columns_1, &idum, &nRhs,
    param, &param[102], &ddum, &ddum, &param[103]);

 FI[cF].nRows=0;
 mxFree(FI[cF].rows_1); mxFree(FI[cF].columns_1); mxFree(FI[cF].pr); 
 FI[cF].rows_1=NULL;

}


/*--------------------------------------------------------------------*/
/*--------------------------------------------------------------------*/
void mexFunction (int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[]) {

char            *buf ;
 int             buflen;
double          *a;
int             j1, j2,j3, nCols, nZmax, nRhs, 
                *rows, *columns, *param;
double          *rhs, *q;
char            statIn[] = "determinant";
int             mtype = -2;   /* Real symmetric matrix */
 int             maxfct, mnum, phase, error, msglvl,nonsym;
int             i;    /* Auxiliary variables. */
double          ddum; /* Double dummy */
int             idum; /* Integer dummy. */

if (nrhs==0) mexErrMsgTxt("no argument");

buflen = (mxGetM (prhs[0]) * mxGetN(prhs[0])) + 1;
buf = mxCalloc (buflen, sizeof(char));
error = mxGetString (prhs[0], buf, buflen);
if (cF==-1) { for (cF=0;cF<20;cF++) {FI[cF].nRows=0;}; cF=0; }

 if ((!strcmp("symbfact",buf)) || (!strcmp("numfact",buf)) || (!strcmp("fact",buf)) || (!strcmp("clear",buf)) ) { 
  if (!mxIsInt32(prhs[2])) mexErrMsgTxt("bad param type");
    param=(int*)mxGetData(prhs[2]); 
} else if (!strcmp("solve",buf)) {
  if (!mxIsInt32(prhs[3])) mexErrMsgTxt("bad param type");
  param=(int*)mxGetData(prhs[3]);
} else if (!strcmp("cvs",buf)) { 
 plhs[0]=mxCreateString("$Revision: 1.2 $  $Date: 2006/05/15 14:06:53 $");
 return;
} else {
  mexErrMsgTxt("bad command\n");
}

maxfct = param[100]; mnum = param[101];  
msglvl = param[102]; error = param[103]; nonsym =param[104]; 


/*----------------------------------------------------------- ismex */
if (!strcmp("ismex",buf))  {

 plhs[0] = mxCreateDoubleMatrix (1, 1, mxREAL);
 a = mxGetPr (plhs[0]);  *a=1.0;


/*--------------------symbolic and numerical fact together ---*/
} else if (!strcmp("fact",buf))  {

cF=0; while (cF<20&&FI[cF].nRows!=0) { cF++; }
if (cF==19) {mexErrMsgTxt("20 factors is the max, use pardiso clear");}
 cF=0; /* xxx can we have more than one pardiso factor ? */

FI[cF].nRows = mxGetM (prhs[1]); nCols = mxGetN (prhs[1]); 
columns  = mxGetIr(prhs[1]);  rows  = mxGetJc(prhs[1]);  


 if (nonsym == 1) {

   nZmax=rows[nCols];
   /* warning pardiso permutes rows/columns */
   FI[cF].columns_1 = mxCalloc(nZmax,sizeof(int)); /* unsymnetric */
   mexMakeMemoryPersistent(FI[cF].columns_1);
   FI[cF].rows_1  = mxCalloc(FI[cF].nRows+1,sizeof(int));
   mexMakeMemoryPersistent(FI[cF].rows_1);
   FI[cF].pr =  mxCalloc(nZmax,sizeof(double)); /* unsymnetric */
   mexMakeMemoryPersistent(FI[cF].pr);
   q=mxGetPr(prhs[1]); 

   j3=0; FI[cF].rows_1[0]=1;for(j1=0; j1<FI[cF].nRows; j1++)  {
     for (j2=rows[j1]; j2<rows[j1+1]; j2++) {  
       FI[cF].columns_1[j3]  = columns[j2]+1;FI[cF].pr[j3]=q[j2]; j3++;  
     }  
     FI[cF].rows_1[j1+1]=j3+1;  
   }

   mtype = 11;
 } else {
   
   nZmax=0; for(j1=0; j1<nCols; j1++)  {
     for(j2=rows[j1]; j2<rows[j1+1]; j2++) {
       if (columns[j2]>=j1) {   
	 nZmax++;
       }
     }
   }

   FI[cF].columns_1 = mxCalloc(nZmax*2,sizeof(int));
   mexMakeMemoryPersistent(FI[cF].columns_1);
   FI[cF].rows_1  = mxCalloc(FI[cF].nRows+1,sizeof(int));
   mexMakeMemoryPersistent(FI[cF].rows_1);
   FI[cF].pr =  mxCalloc(nZmax*2,sizeof(double));
   mexMakeMemoryPersistent(FI[cF].pr);
   q=mxGetPr(prhs[1]); 

   j3=0; FI[cF].rows_1[0]=1;for(j1=0; j1<FI[cF].nRows; j1++)  {
     for (j2=rows[j1]; j2<rows[j1+1]; j2++) {   
       if (columns[j2]>=j1) {   
	 FI[cF].columns_1[j3]  = columns[j2]+1;FI[cF].pr[j3]=q[j2]; j3++;   
       }  
     }   
     FI[cF].rows_1[j1+1]=j3+1;   
   }
   mtype = -2;
 }

plhs[0] = mxCreateDoubleMatrix (1,1, mxREAL);
/*
 for (j1=0;j1<nZmax;j1++) {
 mexPrintf("%i (%g)\n", FI[cF].columns_1[j1],FI[cF].pr[j1]);
 }
*/
if (msglvl>0) {
 mexPrintf("Factor %d :     nNonZeros= %i    nRows=%i \n", 
   cF, nZmax, FI[cF].nRows);
}

/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
for (i = 0; i < 512; i++) {  FI[cF].PardisoPt[i] = 0; }
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
phase = 11; nRhs=0;
PARDISO (FI[cF].PardisoPt, &maxfct, &mnum, &mtype, &phase,
		&(FI[cF].nRows), FI[cF].pr, FI[cF].rows_1, FI[cF].columns_1, &idum, &nRhs,
		param, &msglvl, &ddum, &ddum, &param[103]);
 if (param[103] != 0) return;
 /* {mexErrMsgTxt("\nERROR during symbolic factorization"); }*/

 /* mexPrintf("Reordering completed ... \n");
mexPrintf("Number of nonzeros in factors = %d \n", param[17]);
mexPrintf("Number of factorization MFLOPS = %d \n", param[18]);
 */
 /*
for(j1=0; j1<FI[cF].nRows; j1++)  {
 mexPrintf("\ncol %i : ",j1); 
 for(j2=FI[cF].rows_1[j1]; j2<FI[cF].rows_1[j1+1]; j2++) { 
   mexPrintf("%i (%.1f) ", FI[cF].columns_1[j2-1],FI[cF].pr[j2-1]);
  } 
}
 */

/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
/*phase = 22;mtype = -2; symnetric */
phase = 22;
PARDISO (FI[cF].PardisoPt, &maxfct, &mnum, &mtype, &phase,
   &(FI[cF].nRows), FI[cF].pr, FI[cF].rows_1, FI[cF].columns_1, &idum, &nRhs,
   param, &msglvl, &ddum, &ddum, &param[103]);

if (param[103] != 0) return;
/* { mexErrMsgTxt("\nERROR during numerical factorization"); }*/
if (msglvl>0) mexPrintf("Factorization completed ...\n");
*mxGetPr(plhs[0])=(double)cF;

}

/* --------------------- symbolic factorization------------*/
 else if (!strcmp("symbfact",buf))  {

   cF=0; while (cF<20&&FI[cF].nRows!=0) { cF++; }
   if (cF==20) {mexErrMsgTxt("20 factors is the max, use pardiso clear");}
   cF=0;
   /* xxx can we have more than one pardiso factor ? */
   
   plhs[0] = mxCreateDoubleMatrix (1,1, mxREAL);
   *mxGetPr(plhs[0])=(double)cF;

   PARDISOINIT(FI[cF].PardisoPt, &mtype, param);


   FI[cF].nRows = mxGetM (prhs[1]); nCols = mxGetN (prhs[1]); 
   columns  = mxGetIr(prhs[1]);  rows  = mxGetJc(prhs[1]); 

   if (nonsym == 1) {
     
     nZmax=rows[nCols];
     /* warning pardiso permutes rows/columns */
     FI[cF].columns_1 = mxCalloc(nZmax,sizeof(int)); /* unsymnetric */
     mexMakeMemoryPersistent(FI[cF].columns_1);
     FI[cF].rows_1  = mxCalloc(FI[cF].nRows+1,sizeof(int));
     mexMakeMemoryPersistent(FI[cF].rows_1);
     FI[cF].pr =  mxCalloc(nZmax,sizeof(double)); /* unsymnetric */
     mexMakeMemoryPersistent(FI[cF].pr);
     q=mxGetPr(prhs[1]); 

     j3=0; FI[cF].rows_1[0]=1;for(j1=0; j1<FI[cF].nRows; j1++)  {
       for (j2=rows[j1]; j2<rows[j1+1]; j2++) {  
	 FI[cF].columns_1[j3]  = columns[j2]+1;FI[cF].pr[j3]=q[j2]; j3++;  
       }  
       FI[cF].rows_1[j1+1]=j3+1;  
     }
     
     mtype = 11;
   } else {
     
     nZmax=0; for(j1=0; j1<nCols; j1++)  {
       for(j2=rows[j1]; j2<rows[j1+1]; j2++) {
	 if (columns[j2]>=j1) {   
	   nZmax++;
	 }
       }
     }

     FI[cF].columns_1 = mxCalloc(nZmax*2,sizeof(int));
     mexMakeMemoryPersistent(FI[cF].columns_1);
     FI[cF].rows_1  = mxCalloc(FI[cF].nRows+1,sizeof(int));
     mexMakeMemoryPersistent(FI[cF].rows_1);
     FI[cF].pr =  mxCalloc(nZmax*2,sizeof(double));
     mexMakeMemoryPersistent(FI[cF].pr);
     q=mxGetPr(prhs[1]); 

     j3=0; FI[cF].rows_1[0]=1;for(j1=0; j1<FI[cF].nRows; j1++)  {
       for (j2=rows[j1]; j2<rows[j1+1]; j2++) {   
	 if (columns[j2]>=j1) {   
	   FI[cF].columns_1[j3]  = columns[j2]+1;FI[cF].pr[j3]=q[j2]; j3++;   
	 }  
       }   
       FI[cF].rows_1[j1+1]=j3+1;   
     }
     mtype = -2;
   }

   /*
     for (j1=0;j1<nZmax;j1++) {
     mexPrintf("%i (%g)\n", FI[cF].columns_1[j1],FI[cF].pr[j1]);
     }
   */
   if (msglvl>0) {
     mexPrintf("Factor %d :     nNonZeros= %i    nRows=%i \n", 
	       cF, nZmax, FI[cF].nRows);
   }

/* -------------------------------------------------------------------- */
/* .. Initialize the internal solver memory pointer. This is only */
/* necessary for the FIRST call of the PARDISO solver. */
/* -------------------------------------------------------------------- */
   for (i = 0; i < 512; i++) {  FI[cF].PardisoPt[i] = 0; }
/* -------------------------------------------------------------------- */
/* .. Reordering and Symbolic Factorization. This step also allocates */
/* all memory that is necessary for the factorization. */
/* -------------------------------------------------------------------- */
   phase = 11; nRhs=0;
   PARDISO (FI[cF].PardisoPt, &maxfct, &mnum, &mtype, &phase,
	    &(FI[cF].nRows), FI[cF].pr, FI[cF].rows_1, FI[cF].columns_1, &idum, &nRhs,
	    param, &msglvl, &ddum, &ddum, &param[103]);
   if (param[103] != 0) return;
   if (msglvl>0) mexPrintf("Reordering and symbolic factorization completed ...\n");
   
 }

/* --------------------- nemerical factorization----------------*/
 else if (!strcmp("numfact",buf))  {

/* -------------------------------------------------------------------- */
/* .. Numerical factorization. */
/* -------------------------------------------------------------------- */
    
phase = 22; cF=0;

q=mxGetPr(prhs[1]); 


 if (nonsym == 1) {
   mtype=11;
   j3=0;
   for(j1=0; j1<FI[cF].nRows; j1++)  {
     for (j2= FI[cF].rows_1[j1]-1;j2<(FI[cF].rows_1[j1+1]-1);j2++) {
       FI[cF].pr[j3]=q[j2]; j3++;  
     }
   }
 } else {
   mtype = -2;
   columns  = mxGetIr(prhs[1]);  rows  = mxGetJc(prhs[1]); 
   j3=0;
   for(j1=0; j1<FI[cF].nRows; j1++)  {
     for (j2=rows[j1]; j2<rows[j1+1]; j2++) {   
       if (columns[j2]>=j1) {   
	 FI[cF].pr[j3]=q[j2]; j3++;   
       }  
     }   
   }
 }

 PARDISO (FI[cF].PardisoPt, &maxfct, &mnum, &mtype, &phase,
	  &(FI[cF].nRows), FI[cF].pr, FI[cF].rows_1, FI[cF].columns_1, &idum, &nRhs,
	  param, &msglvl, &ddum, &ddum, &param[103]);

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

 if (param[103] != 0) return;
 /* { mexErrMsgTxt("\nERROR during numerical factorization"); }*/
 if (msglvl>0) {
   mexPrintf("Factor %e :     nNonZeros= %i    nRows=%i \n", FI[cF].pr[1],  FI[cF].rows_1[FI[cF].nRows]-1, FI[cF].nRows);
   mexPrintf("Factorization completed ...\n");
 }
 *mxGetPr(plhs[0])=(double)cF;


 }

/*--------------------------------------------------------- solve */
else if (!strcmp("solve",buf)) {

rhs = mxGetPr(prhs[1]); nRhs=mxGetN(prhs[1]); 

 /* gestion du numero de facteur */
 if (nrhs<3) mexErrMsgTxt("Symbolic factor must be given");
 cF= (int)(*mxGetPr (prhs[2]));
 if (FI[cF].nRows==0) mexErrMsgTxt("Factor is not defined");

 plhs[0] = mxCreateDoubleMatrix (mxGetM(prhs[1]),mxGetN(prhs[1]), mxREAL);
 q=mxGetPr(plhs[0]);

/* -------------------------------------------------------------------- */
/* .. Back substitution and iterative refinement. */
/* -------------------------------------------------------------------- */

 if (nonsym==1) mtype=11;
phase = 33;
PARDISO (FI[cF].PardisoPt, &maxfct, &mnum, &mtype, &phase,
    &(FI[cF].nRows), FI[cF].pr, FI[cF].rows_1, FI[cF].columns_1, &idum, &nRhs,
     param, &msglvl, rhs,q, &param[103]);
if (param[103] != 0) { mexErrMsgTxt("\nERROR during numerical factorization"); }

}/*--------------------------------------------------------- clear */
else if (!strcmp("clear",buf)) {

  if (nrhs<2) {
   cF=-1; if (msglvl>0) mexPrintf("Clearing all spfmex factors");
  } else {
    cF= (int)(*mxGetPr (prhs[1]));
   if (cF<0) if (msglvl>0) mexPrintf("Clearing all spfmex factors");
   else { 
       if (FI[cF].nRows==0) mexErrMsgTxt("Cannot clear undefined factor");
       if (msglvl>0) mexPrintf("Clearing factor %i\n",cF);     
     }
  }
  ClearFactor(cF, maxfct,mnum,mtype,param,msglvl);


}/* end command selection ---------------------------------------- */


 mexAtExit(ClearAllFactors);

 return; 

} /* EOF */

