Changeset 7102 for trunk/psLib/src/math/psMinimizeLMM.c

Timestamp:

May 9, 2006, 2:49:38 PM (20 years ago)

Author:

Paul Price

Message:

Renaming psGaussJordan to psMatrixGJSolve, moving function into psMatrix.c

File:

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

trunk/psLib/src/math/psMinimizeLMM.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-04-21 21:18:44 $
+ *  @version $Revision: 1.14 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-05-10 00:49:38 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -114 +114 @@
     }
 
-    if (false == psGaussJordan(Alpha, Beta)) {
+    if (false == psMatrixGJSolve(Alpha, Beta)) {
         psTrace (__func__, 4, "singular matrix in Guess ABP\n");
         return(false);
@@ -514 +514 @@
 }
 
-// This is a temporary gauss-jordan solver based on gene's
-// version based on the Numerical Recipes version
-bool psGaussJordan(
-    psImage *a,
-    psVector *b)
-{
-    int *indxc,*indxr,*ipiv;
-    int Nx, icol, irow;
-    int i, j, k, l, ll;
-    float big, dum, pivinv;
-    psF64 *vector;
-    psF64 **matrix;
-
-    Nx = a->numCols;
-    matrix = a->data.F64;
-    vector = b->data.F64;
-
-    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) {
-                            psTrace (__func__, 3, "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) {
-            psTrace (__func__, 3, "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);
-}
-
 static void minimizationFree(psMinimization *min)
 {