/*
 * C/Matlab/mex replacement for get_A.m
 *
 * R. Perry, 8 June 2000
 *
 * Note: assumes that zstate = 1, and only uses states(1,1)
 * since all elements of states(:,1) should be the same
 * (unlike get_A.m which is more general and uses flipud).
 *
 * A test using get_ML.m with N=13 shows that use of get_A_C.c
 * decreases the runtime by about 12.5% compared to get_A.m
 */

#include "mex.h"		/* for Matlab/mex stuff */

static int N, Memory;

static double z;		/* states(1,1) */
static double *A;		/* result */
static const double *b;		/* argument */

static void get_A( void)
{
  double *pA;
  const double *pb;
  int i, j;

  for( i = 0; i <= Memory; ++i) {	/* for each column of A */

    pA = A + i*N;
    pb = b;

    for( j = 0; j < i; ++j)
      *pA++ = z;

    for( ; j < N; ++j)
      *pA++ = *pb++;
  }
}

void mexFunction( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
  const mxArray *mx_tmp;
 
/* get z, N, Memory */

  if( (mx_tmp = mexGetArrayPtr( "states", "global")) == NULL)
    mexErrMsgTxt("global variable ``states'' not found");

  z = mxGetScalar( mx_tmp);
  Memory = mxGetM( mx_tmp);

  if( (mx_tmp = mexGetArrayPtr( "N", "global")) == NULL)
    mexErrMsgTxt("global variable ``N'' not found");

  N = mxGetScalar( mx_tmp);

#ifdef DEBUG
  mexPrintf( "z = %g, Memory = %d, N = %d\n", z, Memory, N);
  mexPrintf( "nlhs = %d\n", nlhs);
#endif

/* get rhs */

  if( nrhs != 1)
    mexErrMsgTxt("one argument required");

  b = mxGetPr( prhs[0]);

/* get lhs */

  plhs[0] = mxCreateDoubleMatrix( N, Memory+1, mxREAL);
  A = mxGetPr( plhs[0]);

/* doit */

  get_A();
}