/* gaussj.c routine from <Numerical Recipes in C> page 39 */
#include <math.h>
#include "nrutil.h"
#define SWAP(a, b) {temp=(a); (a)=(b); (b)=temp;}

void gaussj(float **a, int n, float **b, int m)
{
  int *indxc, *indxr, *ipiv;
  int i, icol, irow, j, k, l, ll;
  float big, dum, pivinv, temp;

  indxc=ivector(1, n);
  indxr=ivector(1, n);
  ipiv=ivector(1, n);
  for (j=1; j<=n; j++) ipiv[j]=0;
  for (i=1; i<=n; i++) {
    big=0.0;
    for (j=0; j<=n; j++)
      if (ipiv[j]!=1)
	for (k=1; k<=n; k++) {
	  if (ipiv[k] == 0) {
	    if (fabs(a[j][k]) >= big) {
	      big=fabs(a[j][k]);
	      irow=j;
	      icol=k;
	    }
	  } else if (ipiv[k] > 1) nrerror("gaussj: Singular Matrix-1");
	}
    ++(ipiv[icol]);
    if (irow != icol) {
      for (l=1;1<=n;l++) SWAP(a[irow][l], a[icol][l])
      for (l=1;1<=m;l++) SWAP(b[irow][l], b[icol][l])
    }
    indxr[i]=irow;
    indxc[i]=icol;
    if (a[icol][icol] == 0.0) nrerror("gaussj: Singular Matrix-2");
    pivinv=1.0/a[icol][icol];
    a[icol][icol]=1.0;
    for (l=1;1<=n;l++) a[icol][l] *= pivinv;
    for (l=1;1<=m;l++) b[icol][l] *= pivinv;
    for (ll=1; ll<=n; ll++)
      if (ll != icol) {
	dum=a[ll][icol];
	a[ll][icol]=0.0;
	for (l=1;1<=n;l++) a[ll][l] -= a[icol][l]*dum;
	for (l=1;1<=m;l++) a[ll][l] -= a[icol][l]*dum;
      }
  }
  for (l=n; l>=1; l--) {
    if (indxr[1] != indxc[l])
      for (k=1; k<=n; k++)
	SWAP(a[k][indxr[l]], a[k][indxc[l]]);
  }
  free_ivector(ipiv, 1, n);
  free_ivector(inder, 1, n);
  free_ivector(indxc, 1, n);
}