Changeset 25842 for branches/pap/psLib/src/math/psMatrix.c

Timestamp:

Oct 14, 2009, 12:50:26 PM (17 years ago)

Author:

Paul Price

Message:

lndet works better.

File:

• branches/pap/psLib/src/math/psMatrix.c (modified) (5 diffs)

branches/pap/psLib/src/math/psMatrix.c

@@ -379 +379 @@
     switch (a->type.type) {
       case PS_TYPE_F32: {
-          psF32 **values = a->data.F32; /* Dereference */
-          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
-          for (int i = 0; i < numRows; i++) {
-              for (int j = 0; j < numCols; j++) {
-                  if (!isfinite(values[i][j])) {
-                      // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
-                      // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
-                      // i, j, values[i][j]);
-                      return false;
-                  }
-              }
-          }
-          break;
+          psF32 **values = a->data.F32; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
+                      // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
+                      // i, j, values[i][j]);
+                      return false;
+                  }
+              }
+          }
+          break;
       }
       case PS_TYPE_F64: {
-          psF64 **values = a->data.F64; /* Dereference */
-          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
-          for (int i = 0; i < numRows; i++) {
-              for (int j = 0; j < numCols; j++) {
-                  if (!isfinite(values[i][j])) {
-                      // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
-                      // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
-                      // i, j, values[i][j]);
-                      return false;
-                  }
-              }
-          }
-          break;
+          psF64 **values = a->data.F64; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
+                      // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
+                      // i, j, values[i][j]);
+                      return false;
+                  }
+              }
+          }
+          break;
       }
-        // MATRIX_CHECK_NONFINITE_CASE(F32, a);
-        // MATRIX_CHECK_NONFINITE_CASE(F64, a);
+        // MATRIX_CHECK_NONFINITE_CASE(F32, a);
+        // MATRIX_CHECK_NONFINITE_CASE(F64, a);
       default:
-        psAbort("Should never get here.");
+        psAbort("Should never get here.");
     }
 
@@ -471 +471 @@
     switch (a->type.type) {
       case PS_TYPE_F32: {
-          psF32 **values = a->data.F32; /* Dereference */
-          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
-          for (int i = 0; i < numRows; i++) {
-              for (int j = 0; j < numCols; j++) {
-                  if (!isfinite(values[i][j])) {
-                      return false;
-                  }
-              }
-          }
-          break;
+          psF32 **values = a->data.F32; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      return false;
+                  }
+              }
+          }
+          break;
       }
       case PS_TYPE_F64: {
-          psF64 **values = a->data.F64; /* Dereference */
-          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
-          for (int i = 0; i < numRows; i++) {
-              for (int j = 0; j < numCols; j++) {
-                  if (!isfinite(values[i][j])) {
-                      return false;
-                  }
-              }
-          }
-          break;
+          psF64 **values = a->data.F64; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      return false;
+                  }
+              }
+          }
+          break;
       }
       default:
-        psAbort("Should never get here.");
-    }
-  
+        psAbort("Should never get here.");
+    }
+
     // Following the algorithm laid out by Press et al., we loop along the matrix diagonal, but
     // we do not operate on the diagonal elements in order.  Instead, we are looking for the
@@ -511 +511 @@
 
     if (a->type.type == PS_TYPE_F32) {
-        psF32 **A = a->data.F32;
-        psF32  *B = b->data.F32;
-        int *colIndex = colIndexV->data.S32;
-        int *rowIndex = rowIndexV->data.S32;
-        int *pivot    = pivotV->data.S32;
-        psF32 growth = 1.0;
-
-        for (int diag = 0; diag < nSquare; diag++) {
-
-            psF32 maxval = 0.0;
-            int maxrow = 0;
-            int maxcol = 0;
-
-            // search for the next pivot
-            for (int row = 0; row < nSquare; row++) {
-                if (!isfinite(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 (int col = 0; col < nSquare; col++) {
-                    if (pivot[col]) continue;
-                    if (fabs (A[row][col]) < maxval) continue;
-                    maxval = fabs (A[row][col]);
-                    maxrow = row;
-                    maxcol = col;
-                }
-            }
-
-            // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
-            if (pivot[maxcol]) goto escape;
-            pivot[maxcol] = 1;
-
-            // if the selected pivot is off the diagonal, do a row swap
-            if (maxrow != maxcol) {
-                for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
-                PS_SWAP (B[maxrow], B[maxcol]);
-            }
-            rowIndex[diag] = maxrow;
-            colIndex[diag] = maxcol;
-            if (A[maxcol][maxcol] == 0.0) goto escape;
-            // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
-            // Here we are going to raise an error if the dynamic range is too large.
-
-            /* rescale by pivot reciprocal */
-            psF32 tmpval = 1.0 / A[maxcol][maxcol];
-            A[maxcol][maxcol] = 1.0;
-            for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
-            B[maxcol] *= tmpval;
-
-            // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
-            growth *= tmpval;
-            psTrace ("psLib.math", 4, "growth : %e\n", growth);
-            if (fabs(growth) > MAX_RANGE) goto escape;
-
-            /* adjust the elements above the pivot */
-            for (int row = 0; row < nSquare; row++) {
-                if (row == maxcol) continue;
-                tmpval = A[row][maxcol];
-                A[row][maxcol] = 0.0;
-                for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
-                B[row] -= B[maxcol]*tmpval;
-            }
-        }
-
-        // swap back the inverse matrix based on the row swaps above
-        for (int col = nSquare - 1; col >= 0; col--) {
-            if (rowIndex[col] != colIndex[col]) {
-                for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
-            }
-        }
+        psF32 **A = a->data.F32;
+        psF32  *B = b->data.F32;
+        int *colIndex = colIndexV->data.S32;
+        int *rowIndex = rowIndexV->data.S32;
+        int *pivot    = pivotV->data.S32;
+        psF32 growth = 1.0;
+
+        for (int diag = 0; diag < nSquare; diag++) {
+
+            psF32 maxval = 0.0;
+            int maxrow = 0;
+            int maxcol = 0;
+
+            // search for the next pivot
+            for (int row = 0; row < nSquare; row++) {
+                if (!isfinite(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 (int col = 0; col < nSquare; col++) {
+                    if (pivot[col]) continue;
+                    if (fabs (A[row][col]) < maxval) continue;
+                    maxval = fabs (A[row][col]);
+                    maxrow = row;
+                    maxcol = col;
+                }
+            }
+
+            // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
+            if (pivot[maxcol]) goto escape;
+            pivot[maxcol] = 1;
+
+            // if the selected pivot is off the diagonal, do a row swap
+            if (maxrow != maxcol) {
+                for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
+                PS_SWAP (B[maxrow], B[maxcol]);
+            }
+            rowIndex[diag] = maxrow;
+            colIndex[diag] = maxcol;
+            if (A[maxcol][maxcol] == 0.0) goto escape;
+            // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
+            // Here we are going to raise an error if the dynamic range is too large.
+
+            /* rescale by pivot reciprocal */
+            psF32 tmpval = 1.0 / A[maxcol][maxcol];
+            A[maxcol][maxcol] = 1.0;
+            for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
+            B[maxcol] *= tmpval;
+
+            // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
+            growth *= tmpval;
+            psTrace ("psLib.math", 4, "growth : %e\n", growth);
+            if (fabs(growth) > MAX_RANGE) goto escape;
+
+            /* adjust the elements above the pivot */
+            for (int row = 0; row < nSquare; row++) {
+                if (row == maxcol) continue;
+                tmpval = A[row][maxcol];
+                A[row][maxcol] = 0.0;
+                for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
+                B[row] -= B[maxcol]*tmpval;
+            }
+        }
+
+        // swap back the inverse matrix based on the row swaps above
+        for (int col = nSquare - 1; col >= 0; col--) {
+            if (rowIndex[col] != colIndex[col]) {
+                for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
+            }
+        }
     } else {
-        psF64 **A = a->data.F64;
-        psF64  *B = b->data.F64;
-        int *colIndex = colIndexV->data.S32;
-        int *rowIndex = rowIndexV->data.S32;
-        int *pivot    = pivotV->data.S32;
-        psF64 growth = 1.0;
-
-        for (int diag = 0; diag < nSquare; diag++) {
-
-            psF64 maxval = 0.0;
-            int maxrow = 0;
-            int maxcol = 0;
-
-            // search for the next pivot
-            for (int row = 0; row < nSquare; row++) {
-                if (!isfinite(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 (int col = 0; col < nSquare; col++) {
-                    if (pivot[col]) continue;
-                    if (fabs (A[row][col]) < maxval) continue;
-                    maxval = fabs (A[row][col]);
-                    maxrow = row;
-                    maxcol = col;
-                }
-            }
-
-            // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
-            if (pivot[maxcol]) goto escape;
-            pivot[maxcol] = 1;
-
-            // if the selected pivot is off the diagonal, do a row swap
-            if (maxrow != maxcol) {
-                for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
-                PS_SWAP (B[maxrow], B[maxcol]);
-            }
-            rowIndex[diag] = maxrow;
-            colIndex[diag] = maxcol;
-            if (A[maxcol][maxcol] == 0.0) goto escape;
-            // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
-            // Here we are going to raise an error if the dynamic range is too large.
-
-            /* rescale by pivot reciprocal */
-            psF64 tmpval = 1.0 / A[maxcol][maxcol];
-            A[maxcol][maxcol] = 1.0;
-            for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
-            B[maxcol] *= tmpval;
-
-            // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
-            growth *= tmpval;
-            psTrace ("psLib.math", 4, "growth : %e\n", growth);
-            if (fabs(growth) > MAX_RANGE) goto escape;
-
-            /* adjust the elements above the pivot */
-            for (int row = 0; row < nSquare; row++) {
-                if (row == maxcol) continue;
-                tmpval = A[row][maxcol];
-                A[row][maxcol] = 0.0;
-                for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
-                B[row] -= B[maxcol]*tmpval;
-            }
-        }
-
-        // swap back the inverse matrix based on the row swaps above
-        for (int col = nSquare - 1; col >= 0; col--) {
-            if (rowIndex[col] != colIndex[col]) {
-                for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
-            }
-        }
+        psF64 **A = a->data.F64;
+        psF64  *B = b->data.F64;
+        int *colIndex = colIndexV->data.S32;
+        int *rowIndex = rowIndexV->data.S32;
+        int *pivot    = pivotV->data.S32;
+        psF64 growth = 1.0;
+
+        for (int diag = 0; diag < nSquare; diag++) {
+
+            psF64 maxval = 0.0;
+            int maxrow = 0;
+            int maxcol = 0;
+
+            // search for the next pivot
+            for (int row = 0; row < nSquare; row++) {
+                if (!isfinite(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 (int col = 0; col < nSquare; col++) {
+                    if (pivot[col]) continue;
+                    if (fabs (A[row][col]) < maxval) continue;
+                    maxval = fabs (A[row][col]);
+                    maxrow = row;
+                    maxcol = col;
+                }
+            }
+
+            // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
+            if (pivot[maxcol]) goto escape;
+            pivot[maxcol] = 1;
+
+            // if the selected pivot is off the diagonal, do a row swap
+            if (maxrow != maxcol) {
+                for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
+                PS_SWAP (B[maxrow], B[maxcol]);
+            }
+            rowIndex[diag] = maxrow;
+            colIndex[diag] = maxcol;
+            if (A[maxcol][maxcol] == 0.0) goto escape;
+            // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
+            // Here we are going to raise an error if the dynamic range is too large.
+
+            /* rescale by pivot reciprocal */
+            psF64 tmpval = 1.0 / A[maxcol][maxcol];
+            A[maxcol][maxcol] = 1.0;
+            for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
+            B[maxcol] *= tmpval;
+
+            // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
+            growth *= tmpval;
+            psTrace ("psLib.math", 4, "growth : %e\n", growth);
+            if (fabs(growth) > MAX_RANGE) goto escape;
+
+            /* adjust the elements above the pivot */
+            for (int row = 0; row < nSquare; row++) {
+                if (row == maxcol) continue;
+                tmpval = A[row][maxcol];
+                A[row][maxcol] = 0.0;
+                for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
+                B[row] -= B[maxcol]*tmpval;
+            }
+        }
+
+        // swap back the inverse matrix based on the row swaps above
+        for (int col = nSquare - 1; col >= 0; col--) {
+            if (rowIndex[col] != colIndex[col]) {
+                for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
+            }
+        }
     }
 
@@ -708 +708 @@
     if (determinant) {
       // XXX this is getting the wrong value: is it the wrong calculation?
-      // it disagrees with the results of 
+      // it disagrees with the results of
       // det = (psF32)gsl_linalg_LU_det(lu, signum);
       // used in psMatrixDeterminatn
@@ -753 +753 @@
     // Calculate determinant
     gsl_linalg_LU_decomp(lu, perm, &signum);
-    det = (psF32)gsl_linalg_LU_det(lu, signum);
+    det = (psF32)gsl_linalg_LU_lndet(lu);
 
     // Free GSL structs