Changeset 10251

Timestamp:

Nov 28, 2006, 4:14:49 PM (19 years ago)

Author:

magnier

Message:

added F32 version of psMatrixGJSolve

Location:

trunk/psLib/src/math

Files:

• psMatrix.c (modified) (2 diffs)
• psMatrix.h (modified) (2 diffs)

trunk/psLib/src/math/psMatrix.c

@@ -21 +21 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.44 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-13 21:13:48 $
+ *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-11-29 02:14:49 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -409 +409 @@
 }
 
+// This is a temporary gauss-jordan solver based on gene's
+// version based on the Numerical Recipes version
+bool psMatrixGJSolveF32(psImage *a,
+                        psVector *b
+                       )
+{
+    PS_ASSERT_IMAGE_NON_NULL(a, false);
+    PS_ASSERT_VECTOR_NON_NULL(b, false);
+    PS_ASSERT_IMAGE_TYPE(a, PS_TYPE_F32, false);
+    PS_ASSERT_VECTOR_TYPE(b, PS_TYPE_F32, false);
+    PS_ASSERT_INT_EQUAL(a->numCols, a->numRows, false);
+    int Nx = a->numCols;
+    PS_ASSERT_VECTOR_SIZE(b, (long int)Nx, false);
+
+    psF32 *vector = b->data.F32;
+    psF32 **matrix = a->data.F32;
+    int *indxc = psAlloc(Nx*sizeof(int));
+    int *indxr = psAlloc(Nx*sizeof(int));
+    int *ipiv  = psAlloc(Nx*sizeof(int));
+    memset(ipiv, 0, Nx*sizeof(int));
+
+    int irow = 0;
+    int icol = 0;
+
+    for (int i = 0; i < Nx; i++) {
+        double big = 0.0;
+        for (int j = 0; j < Nx; j++) {
+            if (!isfinite(matrix[i][j])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
+                        "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
+                        i, j, matrix[i][j]);
+                goto fescape;
+            }
+            if (ipiv[j] != 1) {
+                for (int k = 0; k < Nx; k++) {
+                    if (ipiv[k] == 0) {
+                        if (fabs(matrix[j][k]) >= big) {
+                            big  = fabs(matrix[j][k]);
+                            irow = j;
+                            icol = k;
+                        }
+                    } else {
+                        if (ipiv[k] > 1) {
+                            // psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Singular Matrix (1).\n");
+                            psTrace("psLib.math", 4, "Singular Matrix (1).\n");
+                            goto fescape;
+                        }
+                    }
+                }
+            }
+        }
+        ipiv[icol]++;
+        if (irow != icol) {
+            for (int l = 0; l < Nx; l++) {
+                PS_SWAP(matrix[irow][l], matrix[icol][l]);
+            }
+            PS_SWAP(vector[irow], vector[icol]);
+        }
+        indxr[i] = irow;
+        indxc[i] = icol;
+        if (matrix[icol][icol] == 0.0) {
+            // psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Singular Matrix (2).\n");
+            psTrace (__func__, 4, "Singular Matrix (2).\n");
+            goto fescape;
+        }
+        double pivinv = 1.0 / matrix[icol][icol];
+        matrix[icol][icol] = 1.0;
+        for (int l = 0; l < Nx; l++) {
+            matrix[icol][l] *= pivinv;
+        }
+        vector[icol] *= pivinv;
+
+        for (int ll = 0; ll < Nx; ll++) {
+            if (ll != icol) {
+                double dum = matrix[ll][icol];
+                matrix[ll][icol] = 0.0;
+                for (int l = 0; l < Nx; l++) {
+                    matrix[ll][l] -= matrix[icol][l]*dum;
+                }
+                vector[ll] -= vector[icol]*dum;
+            }
+        }
+    }
+
+    for (int l = Nx - 1; l >= 0; l--) {
+        if (indxr[l] != indxc[l]) {
+            for (int k = 0; k < Nx; k++) {
+                PS_SWAP(matrix[k][indxr[l]], matrix[k][indxc[l]]);
+            }
+        }
+    }
+    psFree(ipiv);
+    psFree(indxr);
+    psFree(indxc);
+    return true;
+
+fescape:
+    psFree(ipiv);
+    psFree(indxr);
+    psFree(indxc);
+    return false;
+}
+
 psImage* psMatrixInvert(psImage* out,
                         const psImage* in,

trunk/psLib/src/math/psMatrix.h

@@ -21 +21 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-13 22:27:21 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-11-29 02:14:49 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -69 +69 @@
  */
 bool psMatrixGJSolve(
+    psImage *A,                   ///< Matrix to be solved
+    psVector *b                   ///< Vector of values
+);
+
+/** Gauss-Jordan numerical solver for F32 input data
+ *
+ *  @return bool:   True if successful.
+ */
+bool psMatrixGJSolveF32(
     psImage *A,                   ///< Matrix to be solved
     psVector *b                   ///< Vector of values