/* [+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,z,st] = ipm1mex(H,c,a,x0,z0) */

# define MEXNAME ipm1mex

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

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

 CHECK_NARGIN_EQ(7);

 int n = mxGetNumberOfElements(C);

 CHECK_DOUBLE(H);
 CHECK_NUMEL(H,n*(n+1)/2);
 CHECK_DOUBLE(C);
 CHECK_NUMERIC(A);
 CHECK_SCALAR(A);
 CHECK_DOUBLE(X0);
 CHECK_NUMEL(X0,n);
 CHECK_DOUBLE(Z0);
 CHECK_NUMEL(Z0,n);
 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;
 double *w  = (double *)mxCalloc(n*(n+1)/2+4*n,sizeof(double));
 X = mxCreateDoubleMatrix(n,1,mxREAL);
 Z = mxCreateDoubleMatrix(n,1,mxREAL);
 
 dipm1(mxGetPr(X), mxGetPr(Z), &stat, 
       mxGetPr(H), mxGetPr(C), mxGetScalar(A), mxGetPr(X0), mxGetPr(Z0),
       mxGetScalar(TOL), nit, w, n);
 ST = mxCreateLogicalScalar(stat);
 mxFree(w);

 ASSIGN_PLHS;
}
