/*     M*u=b     */
#include <stdio.h>
#include <stdlib.h>
#include "lupack.h"

#define N 5

float M[N*N]={2,  3,  0,  0,  0,
              3,  0,  4,  0,  6,
              0, -1, -3,  2,  0,
              0,  0,  1,  0,  0,
              0,  4,  2,  0,  1};

float b[N]={8, 45, -3, 3, 19};

int main()
{
  int i;
  int *indx;
  float d;
 
  indx=(int *)malloc(N*sizeof(int));
  ludcmp(M, N, indx, &d);
  lubksb(M, N, indx, b);
  free(indx);

  for(i=0;i<N;i++)
  {
    printf("u[%d]= %g\n", i, b [i]);
  }
  return 0;
}