Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gaussj.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gaussj.c	(revision 37631)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/gaussj.c	(revision 37632)
@@ -1,4 +1,4 @@
 # include <ohana.h>
-# define GROWTHTEST 0
+# define GROWTHTEST 1
 # define MAX_RANGE 1.0e9
 
@@ -48,5 +48,8 @@
     // search for the next pivot
     for (row = 0; row < N; row++) {
-      if (!finite(A[row][diag])) goto escape;
+      if (!finite(A[row][diag])) { 
+	fprintf (stderr, "infinity\n"); 
+	goto escape; 
+      }
 
       // if we have already operated on this row (pivot[row] is true), skip it
@@ -65,16 +68,19 @@
 # if (GROWTHTEST)
     fprintf (stderr, "maxcol: %d\n", maxcol);
-    fprintf (stderr, "full A matrix:\n");
-    for (row = 0; row < N; row++) {
-	for (col = 0; col < N; col++) {
-	    fprintf (stderr, "%10.3e ", A[row][col]);
-	}
-	fprintf (stderr, "\n");
-    }
-    fprintf (stderr, "\n");
+    // fprintf (stderr, "full A matrix:\n");
+    // for (row = 0; row < N; row++) {
+    //   for (col = 0; col < N; col++) {
+    // 	fprintf (stderr, "%10.3e ", A[row][col]);
+    //   }
+    //   fprintf (stderr, "\n");
+    // }
+    // fprintf (stderr, "\n");
 # endif
 
     // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
-    if (pivot[maxcol]) goto escape;
+    if (pivot[maxcol]) { 
+      fprintf (stderr, "singular\n"); 
+      goto escape; 
+    }
     pivot[maxcol] = TRUE;
 
@@ -86,5 +92,9 @@
     rowIndex[diag] = maxrow;
     colIndex[diag] = maxcol;
-    if (A[maxcol][maxcol] == 0.0) goto escape;
+    if (A[maxcol][maxcol] == 0.0) { 
+      fprintf (stderr, "zero pivot\n"); 
+      goto escape;
+    }
+
     // XXX Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow
 
@@ -102,10 +112,13 @@
 #   if (GROWTHTEST)
     fprintf (stderr, "column: %d, maxval : %f, growth: %e, epsilon: %e\n", maxcol, tmpval, growth, epsilon);
-    fprintf (stderr, "A diagonal: ");
-    for (col = 0; col < N; col++) fprintf (stderr, "%f ", A[col][col]);
-    fprintf (stderr, "\n");
-# endif
-
-    if (fabs(growth) > MAX_RANGE) goto escape;
+    // fprintf (stderr, "A diagonal: ");
+    // for (col = 0; col < N; col++) fprintf (stderr, "%f ", A[col][col]);
+    // fprintf (stderr, "\n");
+# endif
+
+    if (fabs(growth) > MAX_RANGE) { 
+      fprintf (stderr, "max range\n"); 
+      goto escape;
+    }
 
     /* adjust the elements above the pivot */
@@ -120,12 +133,12 @@
 
 # if (GROWTHTEST)
-    fprintf (stderr, "final A matrix:\n");
-    for (row = 0; row < N; row++) {
-	for (col = 0; col < N; col++) {
-	    fprintf (stderr, "%10.3e ", A[row][col]);
-	}
-	fprintf (stderr, "\n");
-    }
-    fprintf (stderr, "\n");
+  fprintf (stderr, "final A matrix:\n");
+  // for (row = 0; row < N; row++) {
+  //   for (col = 0; col < N; col++) {
+  //     fprintf (stderr, "%10.3e ", A[row][col]);
+  //   }
+  //   fprintf (stderr, "\n");
+  // }
+  // fprintf (stderr, "\n");
 # endif
 
Index: /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/ohana_allocate.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/ohana_allocate.c	(revision 37631)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/ohana_allocate.c	(revision 37632)
@@ -181,4 +181,7 @@
 
   ref = (Memblock *) in - 1;
+
+  // fprintf (stderr, "  file: %s, line: %d, func: %s, size: %zd, addr: %zx\n", 
+  // ref->file, ref->line, ref->func, ref->size, (size_t) ref);
 
   pthread_mutex_lock(&memBlockListMutex);
Index: /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c
===================================================================
--- /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c	(revision 37631)
+++ /branches/eam_branches/ipp-20140904/Ohana/src/opihi/cmd.data/shfit.c	(revision 37632)
@@ -34,14 +34,40 @@
   ResetVector (Mvec,  OPIHI_FLT, terms->Nterms);
 
-  double **Ar, **br, **Ai, **bi;
-  ALLOCATE (Ar, double *, terms->Nterms);
-  ALLOCATE (Ai, double *, terms->Nterms);
-  ALLOCATE (br, double *, terms->Nterms);
-  ALLOCATE (bi, double *, terms->Nterms);
+  // we only fit the linearlly independent terms: Re(m >= 0), Im(m > 0)
+  int Nre = 0;
+  int Nim = 0;
   for (i = 0; i < terms->Nterms; i++) {
-    ALLOCATE_ZERO (Ar[i], double, terms->Nterms);
-    ALLOCATE_ZERO (Ai[i], double, terms->Nterms);
-    ALLOCATE_ZERO (br[i], double, 1);
-    ALLOCATE_ZERO (bi[i], double, 1);
+    if (terms->m[i] >= 0) Nre ++;
+    if (terms->m[i] >  0) Nim ++;
+  }
+  int *Jre = NULL;
+  int *Jim = NULL;
+  ALLOCATE (Jre, int, Nre);
+  ALLOCATE (Jim, int, Nim);
+
+  Nre = Nim = 0;
+  for (i = 0; i < terms->Nterms; i++) {
+    if (terms->m[i] >= 0) {
+      Jre[Nre] = i;
+      Nre ++;
+    }
+    if (terms->m[i] >  0) {
+      Jim[Nim] = i;
+      Nim ++;
+    }
+  }
+
+  double **Are, **bre, **Aim, **bim;
+  ALLOCATE (Are, double *, Nre);
+  ALLOCATE (bre, double *, Nre);
+  ALLOCATE (Aim, double *, Nim);
+  ALLOCATE (bim, double *, Nim);
+  for (i = 0; i < Nre; i++) {
+    ALLOCATE_ZERO (Are[i], double, Nre);
+    ALLOCATE_ZERO (bre[i], double, 1);
+  }
+  for (i = 0; i < Nim; i++) {
+    ALLOCATE_ZERO (Aim[i], double, Nim);
+    ALLOCATE_ZERO (bim[i], double, 1);
   }
 
@@ -53,87 +79,74 @@
     double Fv = Fvec->elements.Flt[i];
 
-    for (j = 0; j < terms->Nterms; j++) {
-      br[j][0] += Fv * terms->Fr[j];
-      bi[j][0] += Fv * terms->Fi[j];
+    for (j = 0; j < Nre; j++) {
+      int jre = Jre[j];
+      bre[j][0] += Fv * terms->Fr[jre];
+    }
+    for (j = 0; j < Nim; j++) {
+      int jim = Jim[j];
+      bim[j][0] += Fv * terms->Fi[jim];
     }
 
-    for (j = 0; j < terms->Nterms; j++) {
-      for (k = j; k < terms->Nterms; k++) {
-	Ar[j][k] += terms->Fr[j] * terms->Fr[k];
-	Ai[j][k] += terms->Fi[j] * terms->Fi[k];
-      }	
+    for (j = 0; j < Nre; j++) {
+      int jre = Jre[j];
+      for (k = j; k < Nre; k++) {
+	int kre = Jre[k];
+	Are[j][k] += terms->Fr[jre] * terms->Fr[kre];
+      }
     }
-  }
 
-  for (j = 1; j < terms->Nterms; j++) {
-    for (k = 0; k < j; k++) {
-      Ar[j][k] = Ar[k][j];
-      Ai[j][k] = Ai[k][j];
-    }	
-  }
-
-  for (j = 0; j < terms->Nterms; j++) {
-    if (terms->m[j] == 0) {
-      Ai[j][j] = 1.0; // Im(Y_l0) == 0.0
-      bi[j][0] = 0.0;
-    } else {
-      if (terms->m[j] < 0) {
-	for (k = 0; k < terms->Nterms; k++) {
-	  Ar[j][k] = 0.0;
-	  Ar[k][j] = 0.0;
-	  Ai[j][k] = 0.0;
-	  Ai[k][j] = 0.0;
-	}
-	Ar[j][j] = 1.0;
-	Ai[j][j] = 1.0;
-	br[j][0] = 0.0;
-	bi[j][0] = 0.0;
-	// Re(Y_lm) = Re(Y_l(-m)) [m even], Re(Y_lm) = Re(-Y_l(-m)) [m odd]
+    for (j = 0; j < Nim; j++) {
+      int jim = Jim[j];
+      for (k = j; k < Nim; k++) {
+	int kim = Jim[k];
+	Aim[j][k] += terms->Fi[jim] * terms->Fi[kim];
       }
     }
   }
 
-  if (0) {
-    fprintf (stderr, "--- Ar --- : br \n");
-    for (i = 0; i < terms->Nterms; i++) {
-      for (j = 0; j < terms->Nterms; j++) {
-	fprintf (stderr, "%10.6f ", Ar[i][j]);
-      }
-      fprintf (stderr, " : %10.6f\n", br[i][0]);
-    }
-
-    fprintf (stderr, "--- Ai --- : bi \n");
-    for (i = 0; i < terms->Nterms; i++) {
-      for (j = 0; j < terms->Nterms; j++) {
-	fprintf (stderr, "%10.6f ", Ai[i][j]);
-      }
-      fprintf (stderr, " : %10.6f\n", bi[i][0]);
-    }
+  for (j = 1; j < Nre; j++) {
+    for (k = 0; k < j; k++) {
+      Are[j][k] = Are[k][j];
+    }	
+  }
+  for (j = 1; j < Nim; j++) {
+    for (k = 0; k < j; k++) {
+      Aim[j][k] = Aim[k][j];
+    }	
   }
 
-  fprintf (stderr, "--- Ar --- : br \n");
-  for (j = 0; j < terms->Nterms; j++) {
-    fprintf (stderr, "%10.6f ", Ar[j][j]);
-    fprintf (stderr, " : %10.6f\n", br[j][0]);
+  fprintf (stderr, "--- Are --- : bre \n");
+  for (j = 0; j < Nre; j++) {
+    fprintf (stderr, "%10.6f ", Are[j][j]);
+    fprintf (stderr, " : %10.6f\n", bre[j][0]);
   }
 
-  fprintf (stderr, "--- Ai --- : bi \n");
-  for (i = 0; i < terms->Nterms; i++) {
-    fprintf (stderr, "%10.6f : ", Ai[i][i]);
-    fprintf (stderr, " : %10.6f\n", bi[i][0]);
+  fprintf (stderr, "--- Aim --- : bim \n");
+  for (i = 0; i < Nim; i++) {
+    fprintf (stderr, "%10.6f : ", Aim[i][i]);
+    fprintf (stderr, " : %10.6f\n", bim[i][0]);
   }
 
-  if (!dgaussjordan (Ar, br, terms->Nterms, 1)) {
+  if (!dgaussjordan (Are, bre, Nre, 1)) {
     gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
     return FALSE;
   }
-  if (!dgaussjordan (Ai, bi, terms->Nterms, 1)) {
+  if (!dgaussjordan (Aim, bim, Nim, 1)) {
     gprint (GP_ERR, "failed to fit data : ill-conditioned matrix\n");
     return FALSE;
   }
 
-  for (j = 0; j < terms->Nterms; j++) {
-    Fr[j] = br[j][0] / (float) Frvec[0].Nelements;
-    Fi[j] = bi[j][0] / (float) Fivec[0].Nelements;
+  // XXX EAM : why is this factor needed?
+  for (j = 0; j < terms->Nterms; j++) { 
+    Fr[j] = 0.0;
+    Fi[i] = 0.0;
+  }
+  for (j = 0; j < Nre; j++) {
+    int jre = Jre[j];
+    Fr[jre] = bre[j][0] * terms->Nterms / (float) Frvec[0].Nelements;
+  }
+  for (j = 0; j < Nim; j++) {
+    int jim = Jim[j];
+    Fi[jim] = bim[j][0] * terms->Nterms / (float) Frvec[0].Nelements;
   }
 
@@ -142,6 +155,6 @@
     Lvec[0].elements.Flt[j]  = terms->l[j];
     Mvec[0].elements.Flt[j]  = terms->m[j];
-    Frvec[0].elements.Flt[j] = 4*M_PI*Fr[j];
-    Fivec[0].elements.Flt[j] = 4*M_PI*Fi[j];
+    Frvec[0].elements.Flt[j] = Fr[j];
+    Fivec[0].elements.Flt[j] = Fi[j];
   }