Index: trunk/psLib/src/math/psMatrix.c
===================================================================
--- trunk/psLib/src/math/psMatrix.c	(revision 15785)
+++ trunk/psLib/src/math/psMatrix.c	(revision 15820)
@@ -22,9 +22,13 @@
  *  @author Andy Becker, University of Washington (SVD).
  *
- *  @version $Revision: 1.52 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-12-11 22:10:13 $
+ *  @version $Revision: 1.53 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-12-14 00:41:17 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  */
+
+// XXX Future optimisation: use GSL type appropriate to the psLib type, and use memcpy instead of copying each
+// value individually.
+
 
 #ifdef HAVE_CONFIG_H
@@ -312,6 +316,6 @@
 }
 
-// This is a temporary gauss-jordan solver based on gene's
-// version based on the Numerical Recipes version
+// This used to be "a temporary gauss-jordan solver based on gene's version based on the Numerical Recipes
+// version".  However, it's been removed due to copyright, and replaced with LU Decomposition solving.
 bool psMatrixGJSolve(psImage *a,
                      psVector *b
@@ -320,201 +324,53 @@
     PS_ASSERT_IMAGE_NON_NULL(a, false);
     PS_ASSERT_VECTOR_NON_NULL(b, false);
-    PS_ASSERT_IMAGE_TYPE(a, PS_TYPE_F64, false);
-    PS_ASSERT_VECTOR_TYPE(b, PS_TYPE_F64, false);
+    PS_ASSERT_IMAGE_TYPE_F32_OR_F64(a, false);
+    PS_ASSERT_VECTOR_TYPE_F32_OR_F64(b, false);
     PS_ASSERT_INT_EQUAL(a->numCols, a->numRows, false);
-    int Nx = a->numCols;
-    PS_ASSERT_VECTOR_SIZE(b, (long int)Nx, false);
-
-    psF64 *vector = b->data.F64;
-    psF64 **matrix = a->data.F64;
-    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);
+    PS_ASSERT_VECTOR_SIZE(b, (long int)a->numCols, false);
+
+    // Check for non-finite entries
+#define MATRIX_CHECK_NONFINITE_CASE(TYPE,MATRIX) \
+  case PS_TYPE_##TYPE: { \
+      ps##TYPE **values = MATRIX->data.TYPE; /* Dereference */ \
+      int numCols = (MATRIX)->numCols, numRows = (MATRIX)->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; \
+              } \
+          } \
+      } \
+  }
+
+    // Check for non-finite entries in matrix
+    switch (a->type.type) {
+        MATRIX_CHECK_NONFINITE_CASE(F32, a);
+        MATRIX_CHECK_NONFINITE_CASE(F64, a);
+      default:
+        psAbort("Should never get here.");
+    }
+
+    // Decompose the matrix and solve
+    psVector *perm = NULL;              // Permutation vector
+    psImage *lu = psMatrixLUD(NULL, &perm, a); // LU decomposed matrix
+    if (!lu) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to generate LU decomposed matrix");
+        psFree(perm);
+        return false;
+    }
+    b = psMatrixLUSolve(b, lu, b, perm);
+    psFree(lu);
+    psFree(perm);
+    if (!b) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to solve matrix equation.");
+        return false;
+    }
+
     return true;
-
-fescape:
-    psFree(ipiv);
-    psFree(indxr);
-    psFree(indxc);
-    return false;
-}
-
-// 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,