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

 /* function [a,b] = bslvmex(a1,b1,b0,A,B,ni) */

# define MEXNAME bslvmex

# define a  pout[0]
# define b  pout[1]

# define a1 prhs[0]  /* double, def size=na */
# define b1 prhs[1]  /* double, def size=nb */
# define b0 prhs[2]  /* double, size=nb */
# define A  prhs[3]  /* double, size=na*nb */
# define B  prhs[4]  /* double, size=nb*na */
# define ni prhs[5]  /* numeric, scalar */

 CHECK_NARGIN_EQ(6);

 int na = mxGetNumberOfElements(a1), nb = mxGetNumberOfElements(b1);
 
 CHECK_DOUBLE(a1);
 CHECK_DOUBLE(b1);
 CHECK_DOUBLE(b0);
 CHECK_NUMEL(b0,nb);
 CHECK_DOUBLE(A);
 CHECK_NUMEL(A,na*nb);
 CHECK_DOUBLE(B);
 CHECK_NUMEL(B,nb*na);
 CHECK_NUMERIC(ni);
 CHECK_SCALAR(ni);

 CHECK_NARGOUT_LE(2);

 mxArray *pout[2];

 int niter = (int)mxGetScalar(ni);

 a = mxCreateDoubleMatrix(na, 1, mxREAL);
 b = mxCreateDoubleMatrix(nb, 1, mxREAL);
 dbslv(mxGetPr(a), mxGetPr(b), mxGetPr(a1), mxGetPr(b1), mxGetPr(b0), mxGetPr(A), mxGetPr(B), na, nb, niter);

 ASSIGN_PLHS;
}
