# include "data.h"

int gaussjordan (int argc, char **argv) {
  
  float *m, *v;
  double **a, **b;
  int i, j, N;
  Vector *B;
  Buffer *A;

  if (argc != 3) goto usage;

  if ((A = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);    
  if ((B = SelectVector (argv[2], OLDVECTOR, TRUE)) == NULL) return (FALSE);    

  N = B[0].Nelements;
  if (A[0].matrix.Naxis[0] != N) goto usage;
  if (A[0].matrix.Naxis[1] != N) goto usage;
  
  ALLOCATE (a, double *, N);
  ALLOCATE (b, double *, N);
  for (i = 0; i < N; i++) {
    ALLOCATE (a[i], double, N);
    ALLOCATE (b[i], double, 1);
  }

  m = (float *) A[0].matrix.buffer;
  v = B[0].elements;
  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
      a[i][j] = m[i+j*N];
    }
    b[i][0] = v[i]; 
  }
  gaussj (a, N, b, 1);

  for (i = 0; i < N; i++) {
    for (j = 0; j < N; j++) {
       m[i+j*N] = a[i][j];
    }
    v[i] = b[i][0]; 
  }

  for (i = 0; i < N; i++) {
    free (a[i]);
    free (b[i]);
  }
  free (a);
  free (b);
  return (TRUE);

 usage:
  fprintf (stderr, "USAGE: gaussj A B\n");
  fprintf (stderr, "  solves Ax = B, returns 1/A in A and x in B\n");
  fprintf (stderr, "  A must be square, B same dimensions\n");
  return (FALSE);
    
}
