Changeset 7101

Timestamp:

May 9, 2006, 2:48:50 PM (20 years ago)

Author:

Paul Price

Message:

Moving psGaussJordan into psMatrix.c, renaming to psMatrixGJSolve.

File:

• trunk/psLib/src/math/psMatrix.c (modified) (3 diffs)

trunk/psLib/src/math/psMatrix.c

@@ -21 +21 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-07-21 01:40:10 $
+ *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-05-10 00:48:50 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -44 +44 @@
 #include "psConstants.h"
 #include "psErrorText.h"
-
+#include "psTrace.h"
 
 /*****************************************************************************/
@@ -304 +304 @@
 
     return out;
+}
+
+// This is a temporary gauss-jordan solver based on gene's
+// version based on the Numerical Recipes version
+bool psMatrixGJSolve(
+    psImage *a,
+    psVector *b)
+{
+    int *indxc,*indxr,*ipiv;
+    int icol, irow;
+    int i, j, k, l, ll;
+    float big, dum, pivinv;
+    psF64 *vector = b->data.F64;
+    psF64 **matrix = a->data.F64;
+
+    int Nx = a->numCols;
+
+    indxc = psAlloc(Nx*sizeof(int));
+    indxr = psAlloc(Nx*sizeof(int));
+    ipiv  = psAlloc(Nx*sizeof(int));
+    for (j = 0; j < Nx; j++) {
+        ipiv[j] = 0;
+    }
+
+    irow = icol = 0;
+    big = fabs(matrix[0][0]);
+
+    for (i = 0; i < Nx; i++) {
+        big = 0.0;
+        for (j = 0; j < Nx; j++) {
+            if (!isfinite(matrix[i][j])) {
+                psTrace (__func__, 3, "Input matrix contains NaNs: matrix[%d][%d] is %.2f\n", i, j, matrix[i][j]);
+                goto fescape;
+            }
+            if (ipiv[j] != 1) {
+                for (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");
+                            goto fescape;
+                        }
+                    }
+                }
+            }
+        }
+        ipiv[icol]++;
+        if (irow != icol) {
+            for (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");
+            goto fescape;
+        }
+        pivinv = 1.0 / matrix[icol][icol];
+        matrix[icol][icol] = 1.0;
+        for (l = 0; l < Nx; l++) {
+            matrix[icol][l] *= pivinv;
+        }
+        vector[icol] *= pivinv;
+
+        for (ll = 0; ll < Nx; ll++) {
+            if (ll != icol) {
+                dum = matrix[ll][icol];
+                matrix[ll][icol] = 0.0;
+                for (l = 0; l < Nx; l++) {
+                    matrix[ll][l] -= matrix[icol][l]*dum;
+                }
+                vector[ll] -= vector[icol]*dum;
+            }
+        }
+    }
+
+    for (l = Nx - 1; l >= 0; l--) {
+        if (indxr[l] != indxc[l]) {
+            for (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;
 }