Index: /trunk/Ohana/src/libohana/include/ohana.h =================================================================== --- /trunk/Ohana/src/libohana/include/ohana.h (revision 16059) +++ /trunk/Ohana/src/libohana/include/ohana.h (revision 16060) @@ -184,6 +184,6 @@ 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 */ Index: /trunk/Ohana/src/libohana/src/IOBufferOps.c =================================================================== --- /trunk/Ohana/src/libohana/src/IOBufferOps.c (revision 16059) +++ /trunk/Ohana/src/libohana/src/IOBufferOps.c (revision 16060) @@ -45,5 +45,5 @@ 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 */ Index: /trunk/Ohana/src/libohana/src/gaussj.c =================================================================== --- /trunk/Ohana/src/libohana/src/gaussj.c (revision 16059) +++ /trunk/Ohana/src/libohana/src/gaussj.c (revision 16060) @@ -1,155 +1,192 @@ # include -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); } Index: /trunk/Ohana/src/relastro/include/relastro.h =================================================================== --- /trunk/Ohana/src/relastro/include/relastro.h (revision 16059) +++ /trunk/Ohana/src/relastro/include/relastro.h (revision 16060) @@ -207,4 +207,5 @@ 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)); Index: /trunk/Ohana/src/relastro/src/ConfigInit.c =================================================================== --- /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/ConfigInit.c (revision 16060) @@ -3,5 +3,4 @@ void ConfigInit (int *argc, char **argv) { - double ZERO_POINT; char *config, *file; char CatdirPhotcodeFile[256]; Index: /trunk/Ohana/src/relastro/src/FitPM.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPM.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/FitPM.c (revision 16060) @@ -53,5 +53,5 @@ B[3][0] = YT; - dgaussj (A, 4, B, 1); + dgaussjordan (A, B, 4, 1); fit[0].Ro = B[0][0]; Index: /trunk/Ohana/src/relastro/src/FitPMandPar.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/FitPMandPar.c (revision 16060) @@ -76,5 +76,5 @@ B[4][0] = PRX + PDY; - dgaussj (A, 5, B, 1); + dgaussjordan (A, B, 5, 1); fit[0].Ro = B[0][0]; Index: /trunk/Ohana/src/relastro/src/FitPar.c =================================================================== --- /trunk/Ohana/src/relastro/src/FitPar.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/FitPar.c (revision 16060) @@ -7,6 +7,6 @@ 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,5 +50,5 @@ B[2][0] = PRX + PDY; - dgaussj (A, 3, B, 1); + dgaussjordan (A, B, 3, 1); fit[0].Ro = B[0][0]; Index: /trunk/Ohana/src/relastro/src/ImageOps.c =================================================================== --- /trunk/Ohana/src/relastro/src/ImageOps.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/ImageOps.c (revision 16060) @@ -289,4 +289,7 @@ 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(); } } Index: /trunk/Ohana/src/relastro/src/MosaicOps.c =================================================================== --- /trunk/Ohana/src/relastro/src/MosaicOps.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/MosaicOps.c (revision 16060) @@ -22,7 +22,6 @@ 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; Index: /trunk/Ohana/src/relastro/src/UpdateObjects.c =================================================================== --- /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/UpdateObjects.c (revision 16060) @@ -179,5 +179,5 @@ default: - fprintf (stderr, "programming error at %s, %s", __FILE__, __LINE__); + fprintf (stderr, "programming error at %s, %d", __FILE__, __LINE__); exit (2); } Index: /trunk/Ohana/src/relastro/src/fitpoly.c =================================================================== --- /trunk/Ohana/src/relastro/src/fitpoly.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/fitpoly.c (revision 16060) @@ -137,5 +137,5 @@ } - dgaussj (matrix, fit[0].Nelems, vector, 2); + dgaussjordan (matrix, vector, fit[0].Nelems, 2); for (i = 0; i < fit[0].Nelems; i++) { Index: /trunk/Ohana/src/relastro/src/mkpolyterm.c =================================================================== --- /trunk/Ohana/src/relastro/src/mkpolyterm.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/mkpolyterm.c (revision 16060) @@ -61,5 +61,5 @@ beta[1][0] = poly2d_eval (yfit, Nx, Ny, Xo, Yo); - dgaussj (alpha, 2, beta, 1); + dgaussjordan (alpha, beta, 2, 1); Xo -= beta[0][0]; Index: /trunk/Ohana/src/relastro/src/plotstuff.c =================================================================== --- /trunk/Ohana/src/relastro/src/plotstuff.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/plotstuff.c (revision 16060) @@ -79,5 +79,5 @@ float *values; - int i, Nbytes; + int i; if (Npts < 1) return; Index: /trunk/Ohana/src/relastro/src/relastro.c =================================================================== --- /trunk/Ohana/src/relastro/src/relastro.c (revision 16059) +++ /trunk/Ohana/src/relastro/src/relastro.c (revision 16060) @@ -55,5 +55,5 @@ default: - fprintf (stderr, "programming error at %s:%s", __FILE__, __LINE__); + fprintf (stderr, "programming error at %s:%d", __FILE__, __LINE__); exit (2); }