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

 /* u = cholmex(n,a)
    x = cholmex(n,u,b) */

# define MEXNAME cholmex

# define U pout[0]
# define X pout[0]

# define N prhs[0]  /* numeric, scalar */
# define A prhs[1]  /* double, def size=[ma,na] with ma*na>=(n*(n+1))/2 */
# define B prhs[2]  /* double, def size=[*,nb], size=[n,nb] */

 CHECK_NARGIN_GE(2);
 CHECK_NARGIN_LE(3);

 int n=mxGetScalar(N);

 CHECK_NUMERIC(N);
 CHECK_SCALAR(N);
 CHECK_DOUBLE(A);
 CHECK_NUMEL(A,n*(n+1)/2);

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];

 if (nrhs == 2) {
  int ma=mxGetM(A), na=mxGetN(A);
  double *u = mxGetPr(U = mxCreateDoubleMatrix(ma,na,mxREAL));
  memcpy((void *)u,(void *)mxGetPr(A),ma*na*sizeof(double));
  LAPACKE_dpptrf(LAPACK_COL_MAJOR,'U',n,u);
 } else {
  int nb=mxGetN(B);
  CHECK_DOUBLE(B);
  CHECK_NROWS(B,n);
  double *x = mxGetPr(X = mxCreateDoubleMatrix(n,nb,mxREAL));
  memcpy((void *)x,(void *)mxGetPr(B),n*nb*sizeof(double));
  LAPACKE_dpptrs(LAPACK_COL_MAJOR,'U',n,nb,mxGetPr(A),x,n);
 }

 ASSIGN_PLHS;
}
