Index: /branches/eam_branches/20091201/psLib/src/math/psMatrix.c
===================================================================
--- /branches/eam_branches/20091201/psLib/src/math/psMatrix.c	(revision 26343)
+++ /branches/eam_branches/20091201/psLib/src/math/psMatrix.c	(revision 26344)
@@ -1093,8 +1093,6 @@
 }
 
-
-
 // This code supplied by Andy Becker (becker@astro.washington.edu)
-psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in)
+psImage *psMatrixSVD_old(psImage* evec, psVector* eval, const psImage* in)
 {
     #define psMatrixSVD_EXIT {psFree(evec); psFree(eval); return NULL;}
@@ -1145,2 +1143,54 @@
     return evec;
 }
+
+// this is basically a wrapper for the gsl function: gsl_linalg_SV_decomp()
+// SVD decomposes matrix A based on the following equation:  A = U w V^T .
+// This function (as usual for SVD implementations) returns V not V^T.
+// U and V are returned to images; w is returned to a vector.  The input image is not modified.
+// U, V, and w must be supplied as allocated structures, but their lengths are set here to match the 
+// dimensionality of A. 
+// XXX there is no error handling for the gsl functions (anywhere in psMatrix.c)
+bool psMatrixSVD(psImage *U, psVector *w, psImage *V, const psImage *A)
+{
+    // Error checks  Missing one for eval
+    PS_ASSERT_GENERAL_IMAGE_NON_NULL(A, false);
+    PS_CHECK_POINTERS(A, evec, false);
+    PS_CHECK_DIMEN_AND_TYPE(A, PS_DIMEN_IMAGE, false);
+
+    // A is provided with size Nx,Ny = numCols,numRows
+    // U has size Nx,Ny
+    // V has size Nx,Nx
+    // w has size Nx
+
+    // Initialize data
+    int numRows = A->numRows;
+    int numCols = A->numCols;
+
+    U = psImageRecycle(U,  numCols, numRows, A->type.type);
+    V = psImageRecycle(V,  numCols, numCols, A->type.type);
+    w = psVectorRecycle(w, numCols, A->type.type);
+
+    gsl_matrix *Agsl = gsl_matrix_alloc(numRows, numCols);
+    gsl_matrix *Vgsl = gsl_matrix_alloc(numCols, numCols);
+    gsl_vector *Sgsl = gsl_vector_alloc(numCols);
+    gsl_vector *work = gsl_vector_alloc(numCols);
+
+    // Copy psImage data into GSL matrix data
+    matrixPStoGSL(Agsl, A);
+
+    // Calculate SVD decomposition
+    gsl_linalg_SV_decomp(Agsl, Vgsl, Sgsl, work);
+
+    // Copy GSL matrix data to psImage data
+    matrixGSLtoPS(V, Vgsl);
+    matrixGSLtoPS(U, Agsl);  // gsl_linalg_SV_decomp replaces A with U
+    vectorGSLtoPS(S, Sgsl);
+
+    // Free GSL data
+    gsl_matrix_free(A);
+    gsl_matrix_free(V);
+    gsl_vector_free(S);
+    gsl_vector_free(work);
+
+    return true;
+}