Changeset 7103

Timestamp:

May 9, 2006, 3:02:55 PM (20 years ago)

Author:

Paul Price

Message:

Reorganising psMatrixGJSolve.

File:

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

trunk/psLib/src/math/psMatrix.c

@@ -21 +21 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.37 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-05-10 00:48:50 $
+ *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-05-10 01:02:55 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -312 +312 @@
     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++) {
+    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])) {
-                psTrace (__func__, 3, "Input matrix contains NaNs: matrix[%d][%d] is %.2f\n", i, j, 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 (k = 0; k < Nx; k++) {
+                for (int k = 0; k < Nx; k++) {
                     if (ipiv[k] == 0) {
-                        if (fabs (matrix[j][k]) >= big) {
+                        if (fabs(matrix[j][k]) >= big) {
                             big  = fabs(matrix[j][k]);
                             irow = j;
@@ -357 +351 @@
         ipiv[icol]++;
         if (irow != icol) {
-            for (l = 0; l < Nx; l++) {
+            for (int l = 0; l < Nx; l++) {
                 PS_SWAP(matrix[irow][l], matrix[icol][l]);
             }
@@ -368 +362 @@
             goto fescape;
         }
-        pivinv = 1.0 / matrix[icol][icol];
+        double pivinv = 1.0 / matrix[icol][icol];
         matrix[icol][icol] = 1.0;
-        for (l = 0; l < Nx; l++) {
+        for (int l = 0; l < Nx; l++) {
             matrix[icol][l] *= pivinv;
         }
         vector[icol] *= pivinv;
 
-        for (ll = 0; ll < Nx; ll++) {
+        for (int ll = 0; ll < Nx; ll++) {
             if (ll != icol) {
-                dum = matrix[ll][icol];
+                double dum = matrix[ll][icol];
                 matrix[ll][icol] = 0.0;
-                for (l = 0; l < Nx; l++) {
+                for (int l = 0; l < Nx; l++) {
                     matrix[ll][l] -= matrix[icol][l]*dum;
                 }
@@ -387 +381 @@
     }
 
-    for (l = Nx - 1; l >= 0; l--) {
+    for (int l = Nx - 1; l >= 0; l--) {
         if (indxr[l] != indxc[l]) {
-            for (k = 0; k < Nx; k++) {
+            for (int k = 0; k < Nx; k++) {
                 PS_SWAP(matrix[k][indxr[l]], matrix[k][indxc[l]]);
             }