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

 /* b = iatamex(a1,a2) */

# define MEXNAME iatamex

# define B  pout[0]

# define A1 prhs[0] /* double or single, def size=[m1,n] */
# define A2 prhs[1] /* type match A1, def size=[m2,*], size=[m2,n] */

 CHECK_NARGIN_EQ(2);

 int m1 = mxGetM(A1), m2 = mxGetM(A2), n = mxGetN(A1);

 CHECK_REAL(A1);
 CHECK_TYPE_MATCH(A2,A1);
 CHECK_NCOLS(A2,n);

 CHECK_NARGOUT_LE(1);

 mxArray *pout[1];
 
 switch (mxGetClassID(A1)) {
  case mxDOUBLE_CLASS: {
   double *b  = mxGetPr(B = mxCreateDoubleMatrix(n, n, mxREAL));
   double *a1 = mxGetPr(A1);
   double *a2 = mxGetPr(A2);
   switch (n) {
    case 0:  break;
    case 1:  diata1(b, a1, m1, a2, m2); break;
    case 2:  diata2(b, a1, m1, a2, m2); break;
    case 3:  diata3(b, a1, m1, a2, m2); break;
    case 4:  diata4(b, a1, m1, a2, m2); break;
    default: if (diataN(b, a1, m1, a2, m2, n)) mexErrMsgIdAndTxt("iata:lapackerror","LAPACK routines returned with an error");break;
   }
   break;
  }
  case mxSINGLE_CLASS: {
   float *b  = (float *)mxGetData(B = mxCreateNumericMatrix(n, n, mxSINGLE_CLASS, mxREAL));
   float *a1 = (float *)mxGetData(A1);
   float *a2 = (float *)mxGetData(A2);
   switch (n) {
    case 0:  break;
    case 1:  siata1(b, a1, m1, a2, m2); break;
    case 2:  siata2(b, a1, m1, a2, m2); break;
    case 3:  siata3(b, a1, m1, a2, m2); break;
    case 4:  siata4(b, a1, m1, a2, m2); break;
    default: if (siataN(b, a1, m1, a2, m2, n)) mexErrMsgIdAndTxt("iata:lapackerror","LAPACK routines returned with an error");break;
   }
   break;
  }
  /* Other cases impossible due to CHECK_REAL call */
 }

 ASSIGN_PLHS;
}
