#include <stdio.h>
#include "types.h"

int
udu( Real *a,        /* coefficient matrix in udu format */
     Real *y,        /* RHS on input, solution on output */
     const int neq,  /* No. of equations to solve        */
     const int ib,   /* 1/2 bandwidth of the system      */
     const int lu    /* lu=1 factorize, lu=0 resolve     */ );

int
main()
{

  Real *a, *y;
  int neq,ib,lu;
  int i,j;
  lu=1;
  neq=9;
  ib=3;
  a=NULL;
  y=NULL;

  a=(Real *)malloc(neq*ib*sizeof(Real));
  y=(Real *)malloc(neq*sizeof(Real));

  for (j=0;j<ib;j++) {
    for (i=0;i<neq;i++) {
      a[i+j]=0.0;
    }
  }
  for (i=0;i<neq;i++) {
    a[i]=-1.0;
    a[i+neq]=2.0;
    a[i+2*neq]=-1.0;
    y[i]=0.0;
  }
  a[0]=0.0;
  a[3*neq-1]=0.0;
  y[0]=1.0;
  y[neq-1]=0.0;

  for (i=0; i<neq; i++) {
    fprintf(stdout,"a[i,j]= %g %g %g\n",a[i],a[i+neq],a[i+2*neq]);
  }
   
  if (ldu(a,y,neq,ib,lu)<0) {
    fprintf(stdout,"\n\tError termindation from linear solver ...\n");
    exit(0);
  }

  for (i=0;i<neq;i++) {
    fprintf(stdout,"i,y[i]=%4d %10.4e\n",i,y[i]);
  }
  
}