Index: /trunk/Ohana/src/libohana/src/gaussj.c
===================================================================
--- /trunk/Ohana/src/libohana/src/gaussj.c	(revision 8301)
+++ /trunk/Ohana/src/libohana/src/gaussj.c	(revision 8301)
@@ -0,0 +1,155 @@
+# include <ohana.h>
+
+int dgaussj (double **a, int n, double **b, int m) {
+
+  int *indexCol;
+  int *indexRow;
+  int *pivot;
+  int i, col, row, j, k, l, ll;
+  double big, temp;
+  
+  ALLOCATE (indexCol, int, n);
+  ALLOCATE (indexRow, int, n);
+  ALLOCATE (pivot, int, n);
+  for (j = 0; j < n; j++) pivot[j] = 0;
+
+  row = col = 0;
+  big = fabs(a[0][0]);
+
+  for (i = 0; i < n; i++) {
+    big = 0.0;
+    for (j = 0; j < n; j++) {
+      if (!finite(a[i][j])) goto escape;
+      if (pivot[j] != 1) {
+	for (k = 0; k < n; k++) {
+	  if (pivot[k] == 0) {
+	    if (fabs (a[j][k]) >= big) {
+	      big = fabs (a[j][k]);
+	      row = j;
+	      col = k;
+	    }
+	  } else {
+	    if (pivot[k] > 1) goto escape;
+	  }
+	}
+      }
+    }
+    pivot[col]++;
+    if (row != col) {
+      for (l = 0; l < n; l++) SWAP (a[row][l], a[col][l]);
+      for (l = 0; l < m; l++) SWAP (b[row][l], b[col][l]);
+    }
+    indexRow[i] = row;
+    indexCol[i] = col;
+    if (a[col][col] == 0.0) goto escape;
+
+    /* rescale by pivot reciprocal */
+    temp = 1.0 / a[col][col];
+    a[col][col] = 1.0;
+    for (l = 0; l < n; l++) a[col][l] *= temp;
+    for (l = 0; l < m; l++) b[col][l] *= temp;
+
+    /* adjust the elements above the pivot */
+    for (ll = 0; ll < n; ll++) {
+      if (ll != col) {
+	temp = a[ll][col];
+	a[ll][col] = 0.0;
+	for (l = 0; l < n; l++) a[ll][l] -= a[col][l]*temp;
+	for (l = 0; l < m; l++) b[ll][l] -= b[col][l]*temp;
+      }
+    }
+  }
+
+  for (l = n - 1; l >= 0; l--) {
+    if (indexRow[l] != indexCol[l]) {
+      for (k = 0; k < n; k++) SWAP (a[k][indexRow[l]], a[k][indexCol[l]]);
+    }
+  }
+  free (pivot);
+  free (indexRow);
+  free (indexCol);
+  return (TRUE);
+
+escape:
+  free (pivot);
+  free (indexRow);
+  free (indexCol);
+  return (FALSE);
+}
+
+int fgaussj (float **a, int n, float **b, int m) {
+
+  int *indexCol;
+  int *indexRow;
+  int *pivot;
+  int i, col, row, j, k, l, ll;
+  float big, temp;
+  
+  ALLOCATE (indexCol, int, n);
+  ALLOCATE (indexRow, int, n);
+  ALLOCATE (pivot, int, n);
+  for (j = 0; j < n; j++) pivot[j] = 0;
+
+  row = col = 0;
+  big = fabs(a[0][0]);
+
+  for (i = 0; i < n; i++) {
+    big = 0.0;
+    for (j = 0; j < n; j++) {
+      if (!finite(a[i][j])) goto escape;
+      if (pivot[j] != 1) {
+	for (k = 0; k < n; k++) {
+	  if (pivot[k] == 0) {
+	    if (fabs (a[j][k]) >= big) {
+	      big  = fabs (a[j][k]);
+	      row = j;
+	      col = k;
+	    }
+	  } else {
+	    if (pivot[k] > 1) goto escape;
+	  }
+	}
+      }
+    }
+    pivot[col]++;
+    if (row != col) {
+      for (l = 0; l < n; l++) SWAP (a[row][l], a[col][l]);
+      for (l = 0; l < m; l++) SWAP (b[row][l], b[col][l]);
+    }
+    indexRow[i] = row;
+    indexCol[i] = col;
+    if (a[col][col] == 0.0) goto escape;
+
+    /* rescale by pivot reciprocal */
+    temp = 1.0 / a[col][col];
+    a[col][col] = 1.0;
+    for (l = 0; l < n; l++) a[col][l] *= temp;
+    for (l = 0; l < m; l++) b[col][l] *= temp;
+ 
+    /* adjust the elements above the pivot */
+    for (ll = 0; ll < n; ll++) {
+      if (ll != col) {
+	temp = a[ll][col];
+	a[ll][col] = 0.0;
+	for (l = 0; l < n; l++) a[ll][l] -= a[col][l]*temp;
+	for (l = 0; l < m; l++) b[ll][l] -= b[col][l]*temp;
+      }
+    }
+  }
+
+  for (l = n - 1; l >= 0; l--) {
+    if (indexRow[l] != indexCol[l]) {
+      for (k = 0; k < n; k++) SWAP (a[k][indexRow[l]], a[k][indexCol[l]]);
+    }
+  }
+  free (pivot);
+  free (indexRow);
+  free (indexCol);
+  return (TRUE);
+
+escape:
+  free (pivot);
+  free (indexRow);
+  free (indexCol);
+  return (FALSE);
+}
Index: /trunk/Ohana/src/misc/src/fitdist.c
===================================================================
--- /trunk/Ohana/src/misc/src/fitdist.c	(revision 8300)
+++ /trunk/Ohana/src/misc/src/fitdist.c	(revision 8301)
@@ -626,83 +626,2 @@
   
 }
-
-
-/***************************************************************************/
-
-int
-gaussj (a, n, b, m)
-double **a, **b;
-int n,m;
-{
-  int *indxc,*indxr,*ipiv;
-  int i, icol, irow, j, k, l, ll, status;
-  double big,dum,pivinv;
-  
-  ALLOCATE (indxc, int, n);
-  ALLOCATE (indxr, int, n);
-  ALLOCATE (ipiv, int, n);
-  for (j = 0; j < n; j++) 
-    ipiv[j] = 0;
- 
-  for (i = 0; i < n; i++) {
-    big = 0.0;
-    for (j = 0; j < n; j++) {
-      if (ipiv[j] != 1) {
-        for (k = 0; 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) {
-              fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n");
-              return (0);
-            }
-        }
-      }
-    }
-    ipiv[icol]++;
-    if (irow != icol) {
-      for (l = 0; l < n; l++) 
-        SWAP (a[irow][l], a[icol][l]);
-      for (l = 0; l < m; l++) 
-        SWAP (b[irow][l], b[icol][l]);
-    }
-    indxr[i] = irow;
-    indxc[i] = icol;
-    if (a[icol][icol] == 0.0) {
-      fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n");
-      return (0);
-    }
-    pivinv = 1.0 / a[icol][icol];
-    a[icol][icol] = 1.0;
-    for (l = 0; l < n; l++) 
-      a[icol][l] *= pivinv;
-    for (l = 0; l < m; l++) 
-      b[icol][l] *= pivinv;
-    for (ll = 0; ll < n; ll++) {
-      if (ll != icol) {
-        dum = a[ll][icol];
-        a[ll][icol] = 0.0;
-        for (l = 0; l < n; l++) 
-          a[ll][l] -= a[icol][l]*dum;
-        for (l = 0; l < m; l++) 
-          b[ll][l] -= b[icol][l]*dum;
-      }
-    }
-  }
- 
-  for (l = n - 1; l >= 0; l--) {
-    if (indxr[l] != indxc[l])
-      for (k = 0; k < n; k++)
-        SWAP (a[k][indxr[l]], a[k][indxc[l]]);
-  }
-  free (ipiv);
-  free (indxr);
-  free (indxc);
-  return (1);
-}
- 
Index: /trunk/Ohana/src/misc/src/list_astro.c
===================================================================
--- /trunk/Ohana/src/misc/src/list_astro.c	(revision 8300)
+++ /trunk/Ohana/src/misc/src/list_astro.c	(revision 8301)
@@ -97,13 +97,4 @@
   RA = R/N - RX*x/N - RY*y/N;
 
-/*
-  gaussj (a, 4, b, 1);
-
-  RA = b[0][0];
-  RX = b[1][0];
-  RY = b[2][0];
-  RX2 = b[3][0];
-*/
-  
   DX  = (Sdx*Sy2 - Sdy*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
   DY  = (Sdy*Sx2 - Sdx*Sxy) / (Sx2*Sy2 - Sxy*Sxy);
@@ -125,78 +116,2 @@
   exit (0);
 }
-
-
-int gaussj (double **a, int n, double **b, int m) {
-
-  int *indxc,*indxr,*ipiv;
-  int i, icol, irow, j, k, l, ll;
-  double big,dum,pivinv;
-  
-  ALLOCATE (indxc, int, n);
-  ALLOCATE (indxr, int, n);
-  ALLOCATE (ipiv, int, n);
-  for (j = 0; j < n; j++) 
-    ipiv[j] = 0;
-
-  for (i = 0; i < n; i++) {
-    big = 0.0;
-    for (j = 0; j < n; j++) {
-      if (ipiv[j] != 1) {
-	for (k = 0; 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) {
-	      fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n");
-	      return (0);
-	    }
-	}
-      }
-    }
-    ipiv[icol]++;
-    if (irow != icol) {
-      for (l = 0; l < n; l++) 
-	SWAP (a[irow][l], a[icol][l]);
-      for (l = 0; l < m; l++) 
-	SWAP (b[irow][l], b[icol][l]);
-    }
-    indxr[i] = irow;
-    indxc[i] = icol;
-    if (a[icol][icol] == 0.0) {
-      fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n");
-      return (0);
-    }
-    pivinv = 1.0 / a[icol][icol];
-    a[icol][icol] = 1.0;
-    for (l = 0; l < n; l++) 
-      a[icol][l] *= pivinv;
-    for (l = 0; l < m; l++) 
-      b[icol][l] *= pivinv;
-    for (ll = 0; ll < n; ll++) {
-      if (ll != icol) {
-	dum = a[ll][icol];
-	a[ll][icol] = 0.0;
-	for (l = 0; l < n; l++) 
-	  a[ll][l] -= a[icol][l]*dum;
-	for (l = 0; l < m; l++) 
-	  b[ll][l] -= b[icol][l]*dum;
-      }
-    }
-  }
-
-  for (l = n - 1; l >= 0; l--) {
-    if (indxr[l] != indxc[l])
-      for (k = 0; k < n; k++)
-	SWAP (a[k][indxr[l]], a[k][indxc[l]]);
-  }
-  free (ipiv);
-  free (indxr);
-  free (indxc);
-  return (1);
-}
-
Index: /trunk/Ohana/src/misc/src/mosastro.c
===================================================================
--- /trunk/Ohana/src/misc/src/mosastro.c	(revision 8300)
+++ /trunk/Ohana/src/misc/src/mosastro.c	(revision 8301)
@@ -18,5 +18,4 @@
 int mrqinit (int *index, int Npts, double *par, int Npar,                             double (funcs)(int, double *, int, double *));
 int mrqfree (int Npar);
-int gaussj (double **a, int n, double **b, int m);
 int fit (double *X, double *Y, int N, double *A, double *B);
 
Index: /trunk/Ohana/src/mosastro/Makefile
===================================================================
--- /trunk/Ohana/src/mosastro/Makefile	(revision 8300)
+++ /trunk/Ohana/src/mosastro/Makefile	(revision 8301)
@@ -52,5 +52,4 @@
 $(SRC)/output.$(ARCH).o \
 $(SRC)/wstars.$(ARCH).o \
-$(SRC)/gaussj.$(ARCH).o \
 $(SRC)/mkpolyterm.$(ARCH).o \
 $(SRC)/mkmosaic.$(ARCH).o \
Index: /trunk/Ohana/src/mosastro/src/GetGradients.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/GetGradients.c	(revision 8300)
+++ /trunk/Ohana/src/mosastro/src/GetGradients.c	(revision 8301)
@@ -90,5 +90,5 @@
 
 	if (Npts < 5) continue;
-	if (!gaussj (a, 3, b, 2)) continue;
+	if (!dgaussj (a, 3, b, 2)) continue;
 
 	/* we only care about the slopes, not the offsets */
Index: /trunk/Ohana/src/mosastro/src/fitpoly.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/fitpoly.c	(revision 8300)
+++ /trunk/Ohana/src/mosastro/src/fitpoly.c	(revision 8301)
@@ -109,5 +109,5 @@
 # endif
 
-  gaussj (matrix, NPARS, vector, 2); 
+  dgaussj (matrix, NPARS, vector, 2); 
 
 # if (0)
Index: unk/Ohana/src/mosastro/src/gaussj.c
===================================================================
--- /trunk/Ohana/src/mosastro/src/gaussj.c	(revision 8300)
+++ 	(revision )
@@ -1,76 +1,0 @@
-# include "mosastro.h"
-
-int gaussj (double **a, int n, double **b, int m) {
-
-  int *indxc,*indxr,*ipiv;
-  int i, icol, irow, j, k, l, ll;
-  double big,dum,pivinv;
-  
-  ALLOCATE (indxc, int, n);
-  ALLOCATE (indxr, int, n);
-  ALLOCATE (ipiv, int, n);
-  for (j = 0; j < n; j++) 
-    ipiv[j] = 0;
-
-  icol = irow = 0;
-  for (i = 0; i < n; i++) {
-    big = 0.0;
-    for (j = 0; j < n; j++) {
-      if (ipiv[j] != 1) {
-	for (k = 0; 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) {
-	      /* fprintf (stderr, "GAUSSJ: Singular Matrix! (1)\n"); */
-	      return (0);
-	    }
-	}
-      }
-    }
-    ipiv[icol]++;
-    if (irow != icol) {
-      for (l = 0; l < n; l++) 
-	SWAP (a[irow][l], a[icol][l]);
-      for (l = 0; l < m; l++) 
-	SWAP (b[irow][l], b[icol][l]);
-    }
-    indxr[i] = irow;
-    indxc[i] = icol;
-    if (a[icol][icol] == 0.0) {
-      /* fprintf (stderr, "GAUSSJ: Singular Matrix! (2)\n"); */
-      return (0);
-    }
-    pivinv = 1.0 / a[icol][icol];
-    a[icol][icol] = 1.0;
-    for (l = 0; l < n; l++) 
-      a[icol][l] *= pivinv;
-    for (l = 0; l < m; l++) 
-      b[icol][l] *= pivinv;
-    for (ll = 0; ll < n; ll++) {
-      if (ll != icol) {
-	dum = a[ll][icol];
-	a[ll][icol] = 0.0;
-	for (l = 0; l < n; l++) 
-	  a[ll][l] -= a[icol][l]*dum;
-	for (l = 0; l < m; l++) 
-	  b[ll][l] -= b[icol][l]*dum;
-      }
-    }
-  }
-
-  for (l = n - 1; l >= 0; l--) {
-    if (indxr[l] != indxc[l])
-      for (k = 0; k < n; k++)
-	SWAP (a[k][indxr[l]], a[k][indxc[l]]);
-  }
-  free (ipiv);
-  free (indxr);
-  free (indxc);
-  return (1);
-}