Changeset 10381

Timestamp:

Dec 1, 2006, 11:48:47 AM (20 years ago)

Author:

Paul Price

Message:

Adding SVD provided by Andy Becker.

Location:

trunk/psLib/src/math

Files:

• psMatrix.c (modified) (2 diffs)
• psMatrix.h (modified) (2 diffs)

trunk/psLib/src/math/psMatrix.c

@@ -20 +20 @@
  *  @author Ross Harman, MHPCC
  *  @author Robert DeSonia, MHPCC
+ *  @author Andy Becker, University of Washington (SVD).
  *
- *  @version $Revision: 1.45 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-11-29 02:14:49 $
+ *  @version $Revision: 1.46 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-01 21:48:47 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -871 +872 @@
     return outImage;
 }
+
+// This code supplied by Andy Becker (becker@astro.washington.edu)
+psImage *psMatrixSVDSolve(psImage* evec, psVector* eval, const psImage* in)
+{
+    #define psMatrixSVD_EXIT {psFree(evec); psFree(eval); return NULL;}
+
+    // Error checks  Missing one for eval
+    PS_ASSERT_GENERAL_IMAGE_NON_NULL(in, psMatrixSVD_EXIT);
+    PS_CHECK_POINTERS(in, evec, psMatrixSVD_EXIT);
+    PS_CHECK_DIMEN_AND_TYPE(in, PS_DIMEN_IMAGE, psMatrixSVD_EXIT);
+
+    // evec ends up : numCols x numCols
+    evec = psImageRecycle(evec,  in->numCols, in->numCols, in->type.type);
+    eval = psVectorRecycle(eval, in->numCols, in->type.type);
+
+    // Initialize data
+    int numRows = in->numRows;
+    int numCols = in->numCols;
+
+    gsl_matrix *A = gsl_matrix_alloc(numRows, numCols);
+    gsl_matrix *V = gsl_matrix_alloc(numCols, numCols);
+    gsl_vector *S = gsl_vector_alloc(numCols);
+    gsl_vector *work = gsl_vector_alloc(numCols);
+
+    // Copy psImage data into GSL matrix data
+    psImageToGslMatrix(A, in);
+
+    // Calculate SVD decomposition
+    gsl_linalg_SV_decomp(A, V, S, work);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(evec, V);
+    gslVectorToPsVector(eval, S);
+
+    // Take the square root of eval
+    for (int i = 0; i < eval->n; i++) {
+        eval->data.F64[i] = sqrt(eval->data.F64[i]);
+        /* make sure that these things are sorted! */
+        if (i > 0) {
+            assert(eval->data.F64[i] > eval->data.F64[i-1]);
+        }
+    }
+
+    // Free GSL data
+    gsl_matrix_free(A);
+    gsl_matrix_free(V);
+    gsl_vector_free(S);
+    gsl_vector_free(work);
+
+    return evec;
+}

trunk/psLib/src/math/psMatrix.h

@@ -21 +21 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-11-29 02:14:49 $
+ *  @version $Revision: 1.25 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-12-01 21:48:47 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -179 +179 @@
 );
 
+/// Single value decomposition, provided by Andy Becker
+psImage *psMatrixSVDSolve(psImage* evec, psVector* eval, const psImage* in);
+
 /// @}