/* [+MEQ MatlabEQuilibrium Toolbox+] Swiss Plasma Center EPFL Lausanne 2022. All rights reserved. */
# include "meq.h"

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

 /* [Ie,Jh,st] = ipmhmex(Ahd,Aed,Ahe,Aeh,Ahh,uAhh,Yd,Ie0,Jh0,sIp,nit,tol) */

# define MEXNAME ipmhmex

# define IE   pout[ 0]
# define JH   pout[ 1]
# define ST   pout[ 2]

# define AHD  prhs[ 0]  /* double or single, def size=[nh,nd] */
# define AED  prhs[ 1]  /* type match AHD, def size=[ne,*], size=[ne,nd] */
# define AHE  prhs[ 2]  /* type match AHD, size=nh*ne */
# define AEH  prhs[ 3]  /* type match AHD, size=ne*nh */
# define AHH  prhs[ 4]  /* type match AHD, size=nh*nh */
# define UAHH prhs[ 5]  /* type match AHD, size=nh*(nh+1)/2 */
# define YD   prhs[ 6]  /* type match AHD, size=nd */
# define IE0  prhs[ 7]  /* type match AHD, size=ne */
# define JH0  prhs[ 8]  /* type match AHD, size=nh */
# define SIP  prhs[ 9]  /* numeric, scalar */
# define NIT  prhs[10]  /* numeric, scalar */
# define TOL  prhs[11]  /* numeric, scalar */

 CHECK_NARGIN_EQ(12);

 int nh = mxGetM(AHD),
     ne = mxGetM(AED),
     nd = mxGetN(AHD);

 CHECK_REAL(AHD);
 CHECK_TYPE_MATCH(AED,AHD);
 CHECK_NCOLS(AED,nd);
 CHECK_TYPE_MATCH(AHE,AHD);
 CHECK_NUMEL(AHE,nh*ne);
 CHECK_TYPE_MATCH(AEH,AHD);
 CHECK_NUMEL(AEH,ne*nh);
 CHECK_TYPE_MATCH(AHH,AHD);
 CHECK_NUMEL(AHH,nh*nh);
 CHECK_TYPE_MATCH(UAHH,AHD);
 CHECK_NUMEL(UAHH,nh*(nh+1)/2);
 CHECK_TYPE_MATCH(YD,AHD);
 CHECK_NUMEL(YD,nd);
 CHECK_TYPE_MATCH(IE0,AHD);
 CHECK_NUMEL(IE0,ne);
 CHECK_TYPE_MATCH(JH0,AHD);
 CHECK_NUMEL(JH0,nh);
 CHECK_NUMERIC(SIP);
 CHECK_SCALAR(SIP);
 CHECK_NUMERIC(NIT);
 CHECK_SCALAR(NIT);
 CHECK_NUMERIC(TOL);
 CHECK_SCALAR(TOL);

 CHECK_NARGOUT_LE(3);

 mxArray *pout[3];

 int nit = mxGetScalar(NIT);
 bool stat;
 
 switch (mxGetClassID(AHD)) {
  case mxDOUBLE_CLASS: {
   double *w  = (double *)mxCalloc(nh*(nh+1)/2+8*nh+ne,sizeof(double));
   IE = mxCreateDoubleMatrix(ne, 1, mxREAL);
   JH = mxCreateDoubleMatrix(nh, 1, mxREAL);
   dipmh(mxGetPr(IE),mxGetPr(JH),&stat,mxGetPr(AHD),mxGetPr(AED),mxGetPr(AHE),mxGetPr(AEH),mxGetPr(AHH),mxGetPr(UAHH),
         mxGetPr(YD),mxGetPr(IE0),mxGetPr(JH0),mxGetScalar(SIP),nit,mxGetScalar(TOL),w,nh,ne,nd);
   mxFree(w);
   break;
  }
  case mxSINGLE_CLASS: {
   float *w  = (float *)mxCalloc(nh*(nh+1)/2+8*nh+ne,sizeof(float)),
         *sip = (float *)mxGetData(SIP),
         *tol = (float *)mxGetData(TOL);
   IE = mxCreateNumericMatrix(ne, 1, mxSINGLE_CLASS, mxREAL);
   JH = mxCreateNumericMatrix(nh, 1, mxSINGLE_CLASS, mxREAL);
   sipmh((float *)mxGetData(IE),(float *)mxGetData(JH),&stat,(float *)mxGetData(AHD),(float *)mxGetData(AED),(float *)mxGetData(AHE),(float *)mxGetData(AEH),(float *)mxGetData(AHH),(float *)mxGetData(UAHH),
         (float *)mxGetData(YD),(float *)mxGetData(IE0),(float *)mxGetData(JH0),*sip,nit,*tol,w,nh,ne,nd);
   mxFree(w);
   break;
  }
  /* Other cases impossible due to CHECK_REAL call */ 
 }
 ST = mxCreateLogicalScalar(stat);

 ASSIGN_PLHS;
}
