Changeset 16060

Timestamp:

Jan 14, 2008, 1:50:46 PM (18 years ago)

Author:

eugene

Message:

convert to new gauss-jordan code (non-Press); fixing Wall,Werror warnings

Location:

trunk/Ohana/src

Files:

• libohana/include/ohana.h (modified) (1 diff)
• libohana/src/IOBufferOps.c (modified) (1 diff)
• libohana/src/gaussj.c (modified) (1 diff)
• relastro/include/relastro.h (modified) (1 diff)
• relastro/src/ConfigInit.c (modified) (1 diff)
• relastro/src/FitPM.c (modified) (1 diff)
• relastro/src/FitPMandPar.c (modified) (1 diff)
• relastro/src/FitPar.c (modified) (2 diffs)
• relastro/src/ImageOps.c (modified) (1 diff)
• relastro/src/MosaicOps.c (modified) (1 diff)
• relastro/src/UpdateObjects.c (modified) (1 diff)
• relastro/src/fitpoly.c (modified) (1 diff)
• relastro/src/mkpolyterm.c (modified) (1 diff)
• relastro/src/plotstuff.c (modified) (1 diff)
• relastro/src/relastro.c (modified) (1 diff)

trunk/Ohana/src/libohana/include/ohana.h

@@ -184 +184 @@
 char   *ohana_version          PROTO(());
 
-int     dgaussj                PROTO((double **a, int n, double **b, int m));
-int     fgaussj                PROTO((float **a, int n, float **b, int m));
+int     dgaussjordan           PROTO((double **A, double **B, int N, int M));
+int     fgaussjordan           PROTO((float **A, float **B, int n, int m));
 
 /* in time.c */

trunk/Ohana/src/libohana/src/IOBufferOps.c

@@ -45 +45 @@
 
   Nread = read (fd, &buffer[0].buffer[buffer[0].Nbuffer], buffer[0].Nblock);
-  if (DEBUG) fprintf (stderr, "read IO buffer: (%x) %d from %d\n", buffer, Nread, buffer[0].Nblock); 
+  if (DEBUG) fprintf (stderr, "read IO buffer: (%lx) %d from %d\n", (unsigned long) buffer, Nread, buffer[0].Nblock); 
 
   /* on success, increase the block size for the next read */

trunk/Ohana/src/libohana/src/gaussj.c

@@ -1 +1 @@
 # include <ohana.h>
 
-int dgaussj (double **a, int n, double **b, int m) {
+// Gauss-Jordan elimination using full pivots based on Press et al's description.  Major
+// modifications to conform to C indexing, use a boolean to track the completed pivot rows
+// and catch the singular matrix early on.  Also, much cleaner control loops than their
+// implementation.  XXX this really needs to check on round-off errors (see version by
+// William Kahan
+int dgaussjordan (double **A, double **B, int N, int M) {
 
-  int *indexCol;
-  int *indexRow;
+  int *colIndex;
+  int *rowIndex;
   int *pivot;
-  int i, col, row, j, k, l, ll;
-  double big, temp;
+  int maxcol, maxrow;
+  double maxval, tmpval;
   
-  ALLOCATE (indexCol, int, n);
-  ALLOCATE (indexRow, int, n);
-  ALLOCATE (pivot, int, n);
-  for (j = 0; j < n; j++) pivot[j] = 0;
+  int diag, col, row;
 
-  row = col = 0;
-  big = fabs(a[0][0]);
+  ALLOCATE (colIndex, int, N);
+  ALLOCATE (rowIndex, int, N);
+  ALLOCATE (pivot, int, N);
+  memset (pivot, 0, N*sizeof(int));
 
-  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;
-          }
-        }
+  // init these so gcc does not complain
+  maxval = 0.0;
+  maxrow = maxcol = 0;
+
+  // we loop along the matrix diagonal, but we do not operate on the diagonal elements in
+  // order instead, we are looking for the current max element and operating on that
+  // diagonal element.  this is effectively column pivoting.  row pivoting is perfomed
+  // explicitly 
+
+  for (diag = 0; diag < N; diag++) {
+
+    maxval = 0.0;
+
+    // search for the next pivot
+    for (row = 0; row < N; row++) {
+      if (!finite(A[row][diag])) goto escape;
+
+      // if we have already operated on this row (pivot[row] is true), skip it
+      if (pivot[row]) continue;
+
+      // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row
+      for (col = 0; col < N; col++) {
+        if (pivot[col]) continue;
+        if (fabs (A[row][col]) < maxval) continue;
+        maxval = fabs (A[row][col]);
+        maxrow = row;
+        maxcol = col;
       }
     }
-    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]);
+
+    // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
+    if (pivot[maxcol]) goto escape;
+    pivot[maxcol] = TRUE;
+
+    // if the selected pivot is off the diagonal, do a row swap
+    if (maxrow != maxcol) {
+      for (col = 0; col < N; col++) SWAP (A[maxrow][col], A[maxcol][col]);
+      for (col = 0; col < M; col++) SWAP (B[maxrow][col], B[maxcol][col]);
     }
-    indexRow[i] = row;
-    indexCol[i] = col;
-    if (a[col][col] == 0.0) goto escape;
+    rowIndex[diag] = maxrow;
+    colIndex[diag] = maxcol;
+    if (A[maxcol][maxcol] == 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;
+    tmpval = 1.0 / A[maxcol][maxcol];
+    A[maxcol][maxcol] = 1.0;
+    for (col = 0; col < N; col++) A[maxcol][col] *= tmpval;
+    for (col = 0; col < M; col++) B[maxcol][col] *= tmpval;
 
     /* 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 (row = 0; row < N; row++) {
+      if (row == maxcol) continue;
+      tmpval = A[row][maxcol];
+      A[row][maxcol] = 0.0;
+      for (col = 0; col < N; col++) A[row][col] -= A[maxcol][col]*tmpval;
+      for (col = 0; col < M; col++) B[row][col] -= B[maxcol][col]*tmpval;
     }
   }
 
-  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]]);
+  // swap back the inverse matrix based on the row swaps above
+  for (col = N - 1; col >= 0; col--) {
+    if (rowIndex[col] != colIndex[col]) {
+      for (row = 0; row < N; row++) SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
     }
   }
+
   free (pivot);
-  free (indexRow);
-  free (indexCol);
+  free (rowIndex);
+  free (colIndex);
   return (TRUE);
 
 escape:
   free (pivot);
-  free (indexRow);
-  free (indexCol);
+  free (rowIndex);
+  free (colIndex);
   return (FALSE);
 }
 
-int fgaussj (float **a, int n, float **b, int m) {
+int fgaussjordan (float **A, float **B, int N, int M) {
 
-  int *indexCol;
-  int *indexRow;
+  int *colIndex;
+  int *rowIndex;
   int *pivot;
-  int i, col, row, j, k, l, ll;
-  float big, temp;
+  int maxcol, maxrow;
+  float maxval, tmpval;
   
-  ALLOCATE (indexCol, int, n);
-  ALLOCATE (indexRow, int, n);
-  ALLOCATE (pivot, int, n);
-  for (j = 0; j < n; j++) pivot[j] = 0;
+  int diag, col, row;
 
-  row = col = 0;
-  big = fabs(a[0][0]);
+  ALLOCATE (colIndex, int, N);
+  ALLOCATE (rowIndex, int, N);
+  ALLOCATE (pivot, int, N);
+  memset (pivot, 0, N*sizeof(int));
 
-  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;
-          }
-        }
+  // init these so gcc does not complain
+  maxval = 0.0;
+  maxrow = maxcol = 0;
+
+  // we loop along the matrix diagonal, but we do not operate on the diagonal elements in
+  // order instead, we are looking for the current max element and operating on that
+  // diagonal element.  this is effectively column pivoting.  row pivoting is perfomed
+  // explicitly 
+
+  for (diag = 0; diag < N; diag++) {
+
+    maxval = 0.0;
+
+    // search for the next pivot
+    for (row = 0; row < N; row++) {
+      if (!finite(A[row][diag])) goto escape;
+
+      // if we have already operated on this row (pivot[row] is true), skip it
+      if (pivot[row]) continue;
+
+      // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row
+      for (col = 0; col < N; col++) {
+        if (pivot[col]) continue;
+        if (fabs (A[row][col]) < maxval) continue;
+        maxval = fabs (A[row][col]);
+        maxrow = row;
+        maxcol = col;
       }
     }
-    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]);
+
+    // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
+    if (pivot[maxcol]) goto escape;
+    pivot[maxcol] = TRUE;
+
+    // if the selected pivot is off the diagonal, do a row swap
+    if (maxrow != maxcol) {
+      for (col = 0; col < N; col++) SWAP (A[maxrow][col], A[maxcol][col]);
+      for (col = 0; col < M; col++) SWAP (B[maxrow][col], B[maxcol][col]);
     }
-    indexRow[i] = row;
-    indexCol[i] = col;
-    if (a[col][col] == 0.0) goto escape;
+    rowIndex[diag] = maxrow;
+    colIndex[diag] = maxcol;
+    if (A[maxcol][maxcol] == 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;
- 
+    tmpval = 1.0 / A[maxcol][maxcol];
+    A[maxcol][maxcol] = 1.0;
+    for (col = 0; col < N; col++) A[maxcol][col] *= tmpval;
+    for (col = 0; col < M; col++) B[maxcol][col] *= tmpval;
+
     /* 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 (row = 0; row < N; row++) {
+      if (row == maxcol) continue;
+      tmpval = A[row][maxcol];
+      A[row][maxcol] = 0.0;
+      for (col = 0; col < N; col++) A[row][col] -= A[maxcol][col]*tmpval;
+      for (col = 0; col < M; col++) B[row][col] -= B[maxcol][col]*tmpval;
     }
   }
 
-  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]]);
+  // swap back the inverse matrix based on the row swaps above
+  for (col = N - 1; col >= 0; col--) {
+    if (rowIndex[col] != colIndex[col]) {
+      for (row = 0; row < N; row++) SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
     }
   }
+
   free (pivot);
-  free (indexRow);
-  free (indexCol);
+  free (rowIndex);
+  free (colIndex);
   return (TRUE);
 
 escape:
   free (pivot);
-  free (indexRow);
-  free (indexCol);
+  free (rowIndex);
+  free (colIndex);
   return (FALSE);
 }

trunk/Ohana/src/relastro/include/relastro.h

@@ -207 +207 @@
 void unlock_image_db (FITS_DB *db);
 void create_image_db (FITS_DB *db);
+void save_catalogs (Catalog *catalog, int Ncatalog);
 
 int           main                PROTO((int argc, char **argv));

trunk/Ohana/src/relastro/src/ConfigInit.c

@@ -3 +3 @@
 void ConfigInit (int *argc, char **argv) {
 
-  double ZERO_POINT;
   char  *config, *file;
   char CatdirPhotcodeFile[256];

trunk/Ohana/src/relastro/src/FitPM.c

@@ -53 +53 @@
   B[3][0] = YT;
 
-  dgaussj (A, 4, B, 1);
+  dgaussjordan (A, B, 4, 1);
 
   fit[0].Ro = B[0][0];

trunk/Ohana/src/relastro/src/FitPMandPar.c

@@ -76 +76 @@
   B[4][0] = PRX + PDY;
 
-  dgaussj (A, 5, B, 1);
+  dgaussjordan (A, B, 5, 1);
 
   fit[0].Ro = B[0][0];

trunk/Ohana/src/relastro/src/FitPar.c

@@ -7 +7 @@
 
   double **A, **B;
-  double wx, wy, Wx, Wy, Tx, Ty, Tx2, Ty2, Xs, Ys, XT, YT;
-  double PR, PD, PRT, PDT, PRX, PDY, PR2, PD2;
+  double wx, wy, Wx, Wy, Xs, Ys;
+  double PR, PD, PRX, PDY, PR2, PD2;
 
   A = array_init (3, 3);
@@ -50 +50 @@
   B[2][0] = PRX + PDY;
 
-  dgaussj (A, 3, B, 1);
+  dgaussjordan (A, B, 3, 1);
 
   fit[0].Ro = B[0][0];

trunk/Ohana/src/relastro/src/ImageOps.c

@@ -289 +289 @@
       LM_to_XY (&ref[i].X, &ref[i].Y, ref[i].L, ref[i].M, &image[im].coords);
       break;
+      default:
+        fprintf (stderr, "invalid case");
+        abort();
     }
   }

trunk/Ohana/src/relastro/src/MosaicOps.c

@@ -22 +22 @@
 void initMosaics (Image *image, int Nimage) {
 
-  int i, j, status, found, NMOSAIC, *NIMLIST;
+  int i, j, found, NMOSAIC;
   unsigned int start, stop;
-  char *pname;
 
   Nmosaic = 0;

trunk/Ohana/src/relastro/src/UpdateObjects.c

@@ -179 +179 @@
 
         default:
-          fprintf (stderr, "programming error at %s, %s", __FILE__, __LINE__);
+          fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__);
           exit (2);
       }   

trunk/Ohana/src/relastro/src/fitpoly.c

@@ -137 +137 @@
   }     
 
-  dgaussj (matrix, fit[0].Nelems, vector, 2); 
+  dgaussjordan (matrix, vector, fit[0].Nelems, 2); 
 
   for (i = 0; i < fit[0].Nelems; i++) {

trunk/Ohana/src/relastro/src/mkpolyterm.c

@@ -61 +61 @@
       beta[1][0] = poly2d_eval (yfit, Nx, Ny, Xo, Yo);
 
-      dgaussj (alpha, 2, beta, 1);
+      dgaussjordan (alpha, beta, 2, 1);
 
       Xo -= beta[0][0];

trunk/Ohana/src/relastro/src/plotstuff.c

@@ -79 +79 @@
 
   float *values;
-  int i, Nbytes;
+  int i;
 
   if (Npts < 1) return;

trunk/Ohana/src/relastro/src/relastro.c

@@ -55 +55 @@
 
     default:
-      fprintf (stderr, "programming error at %s:%s", __FILE__, __LINE__);
+      fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__);
       exit (2);
   }