/* [+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[])
{

 /* [x,s,z,stat] = ipm2mex(uH,c,A,b,x0,s0,z0,tol,nit) */

# define MEXNAME ipm2mex

# define X   pout[0]
# define S   pout[1]
# define Z   pout[2]
# define ST  pout[3]

# define H   prhs[0]  /* double, size=n*(n+1)/2 */
# define C   prhs[1]  /* double, def size=n */
# define A   prhs[2]  /* double, def size=[ni,*], size=[ni,n] */
# define B   prhs[3]  /* double, size=ni */
# define X0  prhs[4]  /* double, size=n */
# define S0  prhs[5]  /* double, size=ni */
# define Z0  prhs[6]  /* double, size=ni */
# define TOL prhs[7]  /* numeric, scalar */
# define NIT prhs[8]  /* numeric, scalar */

 CHECK_NARGIN_EQ(9);

 int n = mxGetNumberOfElements(C), ni = mxGetM(A);

 CHECK_DOUBLE(H);
 CHECK_NUMEL(H,n*(n+1)/2);
 CHECK_DOUBLE(C);
 CHECK_DOUBLE(A);
 CHECK_NCOLS(A,n);
 CHECK_DOUBLE(B);
 CHECK_NUMEL(B,ni);
 CHECK_DOUBLE(X0);
 CHECK_NUMEL(X0,n);
 CHECK_DOUBLE(S0);
 CHECK_NUMEL(S0,ni);
 CHECK_DOUBLE(Z0);
 CHECK_NUMEL(Z0,ni);
 CHECK_NUMERIC(TOL);
 CHECK_SCALAR(TOL);
 CHECK_NUMERIC(NIT);
 CHECK_SCALAR(NIT);

 CHECK_NARGOUT_LE(4);

 mxArray *pout[4];
 
 bool stat;
 double *w  = (double *)mxCalloc(n*(n+1)/2+2*n+6*ni,sizeof(double));
 X = mxCreateDoubleMatrix(n ,1,mxREAL);
 S = mxCreateDoubleMatrix(ni,1,mxREAL);
 Z = mxCreateDoubleMatrix(ni,1,mxREAL);
 
 dipm2(mxGetPr(X), mxGetPr(S), mxGetPr(Z), &stat, 
       mxGetPr(H), mxGetPr(C), mxGetPr(A), mxGetPr(B),
       mxGetPr(X0), mxGetPr(S0), mxGetPr(Z0),
       mxGetScalar(TOL), (int)mxGetScalar(NIT), w, n, ni);
 ST = mxCreateLogicalScalar(stat);
 mxFree(w);

 ASSIGN_PLHS;
}
