#include "mex.h"
#include "math.h"
#include "matrix.h"
#include "myheader.h"
#include "sparse.h"
#include "morph_ibm1.h"

void mexFunction( int nlhs, mxArray *plhs[],
                  int nrhs, const mxArray *prhs[] )
{
  mxArray *lB,*tmpar,*lalfar,*lbetar,*qa0ar,*ql0ar;
  TENSOR2 U,Uniq,Ucnt,lA,lC,alpha,qa0,ql0,qa1,ql1;
  TENSOR2 tmp,lalf,lbet;
  TENSOR3 PWORD;
  auxparamstype AUXPARAMS;

  makeT2ad(U,       IN(3));
  makeT2ad(Uniq,    IN(4));
  makeT2ad(Ucnt,    IN(5));
  makeT2ad(lA,      IN(6));
  makeT2ad(lC,      IN(7));
  makeT3ad(PWORD,   IN(8));
  makeT2ad(alpha,   IN(9));
  setauxparams(&AUXPARAMS, IN(10));

  /* allocate auxilliary variables */
  lB = mxCreateSparse(d1(PWORD),NLOCAL,d2(U)*NLOCAL,mxREAL);
  tmpar = allocTENSOR2(&tmp,NLOCAL,1); /* must free later */
  lalfar = allocTENSOR2(&lalf,NLOCAL,d2(U)); /* need not free -- passed out */
  lbetar = allocTENSOR2(&lbet,NLOCAL,d2(U)); /* need not free -- passed out */
  qa0ar = allocTENSOR2(&qa0,NTOPIC,d2(Uniq)); /* must free later */
  ql0ar = allocTENSOR2(&ql0,NLOCAL,d2(U)); /* must free later */

  OU(1) = allocTENSOR2(&qa1,NTOPIC,d2(Uniq));
  OU(2) = allocTENSOR2(&ql1,NLOCAL,d2(U));
  OU(3) = lB;
  OU(4) = lalfar;
  OU(5) = lbetar;
  copyarr2tens(     IN(1),qa1);
  copyarr2tens(     IN(2),ql1);

  /* [qa,ql,lB,lalf,lbet] = update_qloop(U,Uniq,Ucnt,qa); */
  update_qloop(&qa1,&ql1,&qa0,&ql0,U,Uniq,Ucnt,lA,lC,PWORD,alpha,lB,tmp,lalf,lbet,AUXPARAMS);

  /* technicalities */
  copytens2arr(qa1,OU(1));
  copytens2arr(ql1,OU(2));

  mxDestroyArray(tmpar);
  mxDestroyArray(qa0ar);
  mxDestroyArray(ql0ar);
  freeauxparams(AUXPARAMS);
}