Changeset 2671

Timestamp:

Dec 9, 2004, 10:51:35 AM (22 years ago)

Author:

harman

Message:

Added support for psF32 types

Location:

trunk/psLib/src

Files:

• dataManip/psDataManipErrors.dat (modified) (1 diff)
• dataManip/psDataManipErrors.h (modified) (3 diffs)
• dataManip/psMatrix.c (modified) (20 diffs)
• dataManip/psMatrix.h (modified) (10 diffs)
• math/psMatrix.c (modified) (20 diffs)
• math/psMatrix.h (modified) (10 diffs)

trunk/psLib/src/dataManip/psDataManipErrors.dat

@@ -52 +52 @@
 psMatrix_VECTOR_EMPTY                  Input psVector contains no elements.  No data to perform operation with.
 psMatrix_IMAGE_EMPTY                   Input psImage contains no pixels.  No data to perform operation with.
+psMatrix_TRANSPOSE_MISMATCH            Number of rows do not match number of columns.

trunk/psLib/src/dataManip/psDataManipErrors.h

@@ -7 +7 @@
  *  @author Robert DeSonia, MHPCC
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-24 21:59:43 $
+ *  @version $Revision: 1.10 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-09 20:51:35 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -23 +23 @@
  *     $2  The error text (rest of the line in psDataManipErrors.dat)
  *     $n  The order of the source line in psDataManipErrors.dat (comments excluded)
- * 
+ *
  * DO NOT EDIT THE LINES BETWEEN //~Start and //~End!  ANY CHANGES WILL BE OVERWRITTEN.
  */
@@ -67 +67 @@
 #define PS_ERRORTEXT_psMatrix_VECTOR_EMPTY "Input psVector contains no elements.  No data to perform operation with."
 #define PS_ERRORTEXT_psMatrix_IMAGE_EMPTY "Input psImage contains no pixels.  No data to perform operation with."
+#define PS_ERRORTEXT_psMatrix_TRANSPOSE_MISMATCH "Number of rows do not match number of columns."
 //~End
 

trunk/psLib/src/dataManip/psMatrix.c

@@ -1 +1 @@
 /** @file  psMatrix.c
  *
- *  @brief Provides functions for linear algebra operations on psImages and psVectors. 
+ *  @brief Provides functions for linear algebra operations on psImages and psVectors.
  *
  *  Functions are provided to:
@@ -19 +19 @@
  *
  *  @author Ross Harman, MHPCC
- *   
- *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-04 01:04:59 $
+ *
+ *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-09 20:51:16 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -41 +41 @@
 #include "psMatrix.h"
 #include "psConstants.h"
-
-
+#include "psDataManipErrors.h"
+
+
+/*****************************************************************************/
+/* DEFINE STATEMENTS                                                         */
+/*****************************************************************************/
 
 /** Preprocessor macro to generate error for image dimensionality not set to PS_DIMEN_IMAGE */
@@ -50 +54 @@
             "Invalid operation. %s has incorrect dimensionality %d.", #NAME, PS_DIMEN);             \
     return RETURN;                                                                                  \
-} else if(NAME->type.type != PS_TYPE_F64) {                                                         \
+} else if(NAME->type.type!=PS_TYPE_F64 && NAME->type.type!=PS_TYPE_F32) {                           \
     psError(PS_ERR_BAD_PARAMETER_TYPE, true,                                                        \
             "Invalid operation. %s not PS_TYPE_F64.", #NAME);                                       \
@@ -67 +71 @@
 #define PS_CHECK_SQUARE(NAME, RETURN)                                                               \
 if (NAME->numCols != NAME->numRows) {                                                               \
-    psError(PS_ERR_BAD_PARAMETER_SIZE, true,                                                       \
-            "Invalid operation: %s not square array.", #NAME);                             \
+    psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: %s not square array.", #NAME);     \
     return RETURN;                                                                                  \
 }
@@ -79 +82 @@
 LHS_NAME.data  = RHS_NAME;
 
+
+/*****************************************************************************/
+/* FILE STATIC FUNCTIONS                                                     */
+/*****************************************************************************/
+
+static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector);
+static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector);
+static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage);
+static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix);
+
+/** Static function to copy psF32 or psF64 vector data to a GSL vector */
+static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector)
+{
+    psU32 i = 0;
+    psU32 n = 0;
+
+
+    n = inVector->n;
+    for(i=0; i<n; i++) {
+        if(inVector->type.type == PS_TYPE_F32) {
+            outGslVector->data[i] = (psS32)inVector->data.F32[i];
+        } else {
+            outGslVector->data[i] = inVector->data.F64[i];
+        }
+    }
+}
+
+/** Static function to copy GSL vector data to a psF32 or psF64 vector */
+static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector)
+{
+    psU32 i = 0;
+    psU32 n = 0;
+
+
+    n = outVector->n;
+    for(i=0; i<n; i++) {
+        if(outVector->type.type == PS_TYPE_F32) {
+            outVector->data.F32[i] = (psF32)inGslVector->data[i];
+        } else {
+            outVector->data.F64[i] = inGslVector->data[i];
+        }
+    }
+}
+
+/** Static function to copy psF32 or psF64 image data to a GSL matrix */
+static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage)
+{
+    psU32 i = 0;
+    psU32 j = 0;
+    psU32 numRows = 0;
+    psU32 numCols = 0;
+
+
+    numRows = inImage->numRows;
+    numCols = inImage->numCols;
+    for(i=0; i<numRows; i++) {
+        for(j=0; j<numCols; j++) {
+            if(inImage->type.type == PS_TYPE_F32) {
+                outGslMatrix->data[i*numCols+j] = (psS32)inImage->data.F32[i][j];
+            } else {
+                outGslMatrix->data[i*numCols+j] = inImage->data.F64[i][j];
+            }
+        }
+    }
+}
+
+/** Static function to copy GSL matrix data to a psF32 or psF64 image */
+static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix)
+{
+    psU32 i = 0;
+    psU32 j = 0;
+    psU32 numRows = 0;
+    psU32 numCols = 0;
+
+
+    numRows = outImage->numRows;
+    numCols = outImage->numCols;
+    for(i=0; i<numRows; i++) {
+        for(j=0; j<numCols; j++) {
+            if(outImage->type.type == PS_TYPE_F32) {
+                outImage->data.F32[i][j] = (psF32)inGslMatrix->data[i*numCols+j];
+            } else {
+                outImage->data.F64[i][j] = inGslMatrix->data[i*numCols+j];
+            }
+        }
+    }
+}
+
+
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
@@ -86 +178 @@
 {
     psS32 signum = 0;
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix lu;
+    gsl_matrix *lu = NULL;
     gsl_permutation perm;
 
-    #define psMatrixLUD_EXIT { \
-                               psFree(outImage); \
-                               return NULL; \
-                             }
-
-    // Error checks
+
+    #define psMatrixLUD_EXIT {psFree(outImage); return NULL;}
+
+    // Error checks
+    PS_IMAGE_CHECK_NULL_GENERAL(inImage, psMatrixLUD_EXIT);
     PS_CHECK_POINTERS(inImage, outImage, outImage);
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
-
-    PS_IMAGE_CHECK_NULL_GENERAL(inImage, psMatrixLUD_EXIT);
     PS_VECTOR_CHECK_NULL_GENERAL(outPerm, psMatrixLUD_EXIT);
-
     PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage);
-
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     psVectorRecycle(outPerm, inImage->numRows, inImage->type.type);
+    PS_CHECK_SQUARE(inImage, outImage);
+    PS_CHECK_SQUARE(outImage, outImage);
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
 
     // Initialize GSL data
@@ -118 +205 @@
     outPerm->n = numCols;
     perm.data = outPerm->data.V;
-    PS_GSL_MATRIX_INITIALIZE(lu, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    // Copy psImage input data into GSL matrix data to keep input data pristine
-    memcpy(lu.data, inImage->data.V[0], arraySize);
+    lu = gsl_matrix_alloc(numRows, numCols);
+
+    // Copy psImage data into GSL matrix data
+    psImageToGslMatrix(lu, inImage);
 
     // Calculate LU decomposition
-    gsl_linalg_LU_decomp(&lu, &perm, &signum);
+    gsl_linalg_LU_decomp(lu, &perm, &signum);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, lu);
+
+    // Free GSL data
+    gsl_matrix_free(lu);
 
     return outImage;
 }
 
-psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage,
-                          const psVector* inVector, const psVector* inPerm)
-{
-    psS32 arraySize = 0;
+psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage, const psVector* inVector,
+                          const psVector* inPerm)
+{
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix lu;
+    gsl_matrix *lu;
     gsl_permutation perm;
-    gsl_vector b;
-    gsl_vector x;
+    gsl_vector *b = NULL;
+    gsl_vector *x = NULL;
 
     // Error checks
@@ -157 +245 @@
     PS_VECTOR_CHECK_NULL(inPerm, outVector);
     PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector);
-
     psVectorRecycle(outVector, inImage->numRows, inImage->type.type);
-
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
 
     // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(lu, inImage->data.V[0]);
+    lu = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(lu, inImage);
+    b = gsl_vector_alloc(inVector->n);
+    psVectorToGslVector(b, inVector);
+    x = gsl_vector_alloc(inVector->n);
 
     outVector->n = numCols;
-
     perm.size = inPerm->n;
     perm.data = inPerm->data.V;
 
-    b.size = inVector->n;
-    b.stride = 1;
-    b.data = inVector->data.V;
-
-    x.size = numCols;
-    x.stride = 1;
-    x.data = outVector->data.V;
-
     // Solve for {x} in equation: {b} = [A]{x}
-    gsl_linalg_LU_solve(&lu, &perm, &b, &x);
+    gsl_linalg_LU_solve(lu, &perm, b, x);
+
+    // Copy GSL vector data to psVector data
+    gslVectorToPsVector(outVector, x);
+
+    // Free GSL data
+    gsl_vector_free(b);
+    gsl_vector_free(x);
+    gsl_matrix_free(lu);
 
     return outVector;
 }
 
-psImage* psMatrixInvert(psImage* outImage, const psImage* inImage, float *restrict det)
+psImage* psMatrixInvert(psImage* outImage, const psImage* inImage, psF32 *det)
 {
     psS32 signum = 0;
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix inv;
-    gsl_matrix *lu;
+    gsl_matrix *inv = NULL;
+    gsl_matrix *lu = NULL;
     gsl_permutation *perm = NULL;
 
@@ -204 +291 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     PS_IMAGE_CHECK_EMPTY(inImage, outImage);
-
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
+    PS_CHECK_SQUARE(inImage, outImage);
+    PS_CHECK_SQUARE(outImage, outImage);
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
-
-    // Allocate GSL structs
-    perm = gsl_permutation_alloc(inImage->numRows);
+
+    // Initialize GSL data
+    perm = gsl_permutation_alloc(numRows);
     lu = gsl_matrix_alloc(numRows, numCols);
-
-    // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(inv, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    // Copy psImage input data into GSL matrix data to keep input data pristine
-    memcpy(lu->data, inImage->data.V[0], arraySize);
+    inv = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(lu, inImage);
 
     // Invert data and calculate determinant
     gsl_linalg_LU_decomp(lu, perm, &signum);
-    gsl_linalg_LU_invert(lu, perm, &inv);
+    gsl_linalg_LU_invert(lu, perm, inv);
     *det = (float)gsl_linalg_LU_det(lu, signum);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, inv);
 
     // Free GSL structs
     gsl_permutation_free(perm);
     gsl_matrix_free(lu);
+    gsl_matrix_free(inv);
 
     return outImage;
 }
 
-float *psMatrixDeterminant(const psImage* restrict inImage)
+psF32 *psMatrixDeterminant(const psImage* inImage)
 {
     psS32 signum = 0;
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    float *det = NULL;
+    psF32 *det = NULL;
     gsl_matrix *lu = NULL;
     gsl_permutation *perm = NULL;
@@ -252 +334 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, NULL);
     PS_IMAGE_CHECK_EMPTY(inImage, NULL);
+    PS_CHECK_SQUARE(inImage, 0);
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
 
     // Allocate GSL structs
     perm = gsl_permutation_alloc(numRows);
     lu = gsl_matrix_alloc(numRows, numCols);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, 0);
-
-    // Copy psImage input data into GSL matrix data to keep input data pristine
-    memcpy(lu->data, inImage->data.V[0], arraySize);
+    psImageToGslMatrix(lu, inImage);
 
     // Calculate determinant
-    det = (float *)psAlloc(sizeof(float));
+    det = (psF32*)psAlloc(sizeof(psF32));
     gsl_linalg_LU_decomp(lu, perm, &signum);
-    *det = (float)gsl_linalg_LU_det(lu, signum);
+    *det = (psF32)gsl_linalg_LU_det(lu, signum);
 
     // Free GSL structs
@@ -282 +359 @@
 psImage* psMatrixMultiply(psImage* outImage, psImage* inImage1, psImage* inImage2)
 {
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix m1;
-    gsl_matrix m2;
-    gsl_matrix m3;
+    gsl_matrix *m1 = NULL;
+    gsl_matrix *m2 = NULL;
+    gsl_matrix *m3 = NULL;
 
     // Error checks
@@ -299 +375 @@
     PS_IMAGE_CHECK_EMPTY(inImage2, outImage);
     PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage);
-
     psImageRecycle(outImage, inImage2->numCols, inImage2->numRows, inImage2->type.type);
+    PS_CHECK_SQUARE(inImage1, outImage);
+    PS_CHECK_SQUARE(inImage2, outImage);
+    PS_CHECK_SQUARE(outImage, outImage);
 
     // Initialize data
     numRows = inImage1->numRows;
     numCols = inImage1->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(outImage->type.type) * numRows * numCols;
 
     // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(m1, inImage1->data.V[0]);
-    PS_GSL_MATRIX_INITIALIZE(m2, inImage2->data.V[0]);
-    PS_GSL_MATRIX_INITIALIZE(m3, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage1, outImage);
-    PS_CHECK_SQUARE(inImage2, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    gsl_linalg_matmult(&m1, &m2, &m3);
+    m1 = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(m1, inImage1);
+    m2 = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(m2, inImage2);
+    m3 = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(m3, outImage);
+
+    // Perform multiplication
+    gsl_linalg_matmult(m1, m2, m3);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, m3);
+
+    // Free GSL structs
+    gsl_matrix_free(m1);
+    gsl_matrix_free(m2);
+    gsl_matrix_free(m3);
 
     return outImage;
@@ -324 +408 @@
 psImage* psMatrixTranspose(psImage* outImage, const psImage* inImage)
 {
-    psS32 arraySize = 0;
-    psS32 numRows = 0;
-    psS32 numCols = 0;
-    gsl_matrix trans;
+    psU32 i = 0;
+    psU32 j = 0;
+    psS32 numRowsIn = 0;
+    psS32 numColsIn = 0;
+    psS32 numRowsOut = 0;
+    psS32 numColsOut = 0;
+
 
     // Error checks
@@ -334 +421 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     PS_IMAGE_CHECK_EMPTY(inImage, outImage);
-
-    psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
+    outImage = psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
 
     // Initialize data
-    numRows = inImage->numRows;
-    numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
-
-    // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(trans, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    // Copy psImage input data into psImage output data to keep input data pristine
-    memcpy(outImage->data.V[0], inImage->data.V[0], arraySize);
-
-    // Transpose data
-    gsl_matrix_transpose(&trans);
+    numRowsIn = inImage->numRows;
+    numColsIn = inImage->numCols;
+    numRowsOut = outImage->numRows;
+    numColsOut = outImage->numCols;
+
+    if(numRowsIn!=numColsOut && numRowsOut!=numColsIn) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psMatrix_TRANSPOSE_MISMATCH);
+    }
+
+    if(outImage->type.type == PS_TYPE_F32) {
+        for(i=0; i<numRowsOut; i++) {
+            for(j=0; j<numColsOut; j++) {
+                outImage->data.F32[i][j] = inImage->data.F32[j][i];
+            }
+        }
+    } else {
+        for(i=0; i<numRowsOut; i++) {
+            for(j=0; j<numColsOut; j++) {
+                outImage->data.F64[i][j] = inImage->data.F64[j][i];
+            }
+        }
+    }
 
     return outImage;
@@ -364 +456 @@
     gsl_vector *eVals = NULL;
     gsl_eigen_symmv_workspace *w = NULL;
-    gsl_matrix out;
-    gsl_matrix in;
+    gsl_matrix *out = NULL;
+    gsl_matrix *in = NULL;
+
 
     // Error checks
@@ -372 +465 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     PS_IMAGE_CHECK_EMPTY(inImage, outImage);
-
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
 
@@ -379 +471 @@
     numCols = inImage->numCols;
 
-    // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(in, inImage->data.V[0]);
-    PS_GSL_MATRIX_INITIALIZE(out, outImage->data.V[0]);
+    in = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(in, inImage);
+    out = gsl_matrix_alloc(numRows, numCols);
 
     // Allocate GSL structs
@@ -392 +484 @@
 
     // Calculate Eigenvalues and Eigenvectors...Eigenvalues not currently used
-    gsl_eigen_symmv(&in, eVals, &out, w);
+    gsl_eigen_symmv(in, eVals, out, w);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, out);
 
     // Free GSL structs
+    gsl_matrix_free(in);
+    gsl_matrix_free(out);
     gsl_eigen_symmv_free(w);
     gsl_vector_free(eVals);
@@ -405 +502 @@
     psS32 size = 0;
 
-    #define psMatrixToVector_EXIT { \
-                                    psFree(outVector); \
-                                    return NULL; \
-                                  }
+    #define psMatrixToVector_EXIT {psFree(outVector); return NULL;}
 
     // Error checks

trunk/psLib/src/dataManip/psMatrix.h

@@ -2 +2 @@
 /** @file  psMatrix.h
  *
- *  @brief Provides functions for linear algebra operations on psImages and psVectors. 
+ *  @brief Provides functions for linear algebra operations on psImages and psVectors.
  *
  *  Functions are provided to:
@@ -22 +22 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-21 23:44:10 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-09 20:51:22 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -37 +37 @@
  *
  *  Performs a LU decomposition on a psImage matrix and returns the LU matrix. If the user specifies NULL for
- *  the outImage or outPerm arguments, then they will be automatically created. The input image must 
- *  be square. This function operates only with the psF64 data type. Input and output arguments should not be 
- *  the same. GSL indexes the top row as the zero row, not the bottom.  
+ *  the outImage or outPerm arguments, then they will be automatically created. The input image must
+ *  be square. This function operates only with the psF64 data type. Input and output arguments should not be
+ *  the same. GSL indexes the top row as the zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to LU decomposed psImage.
@@ -51 +51 @@
 /** LU Solution of psImage matrix.
  *
- *  Solves for and returns the psVector, {x} in the equation [A]{x} = {b}. If the user specifies NULL as the 
- *  outVector argument, then it will automatically be created. The input image must be square. This function 
- *  operates only with the psF64 data type. Input and output arguments should not be the same. GSL indexes 
- *  the top row as the zero row, not the bottom.  
+ *  Solves for and returns the psVector, {x} in the equation [A]{x} = {b}. If the user specifies NULL as the
+ *  outVector argument, then it will automatically be created. The input image must be square. This function
+ *  operates only with the psF64 data type. Input and output arguments should not be the same. GSL indexes
+ *  the top row as the zero row, not the bottom.
  *
  *  @return  psVector* : Pointer to psVector solution of matrix equation.
@@ -68 +68 @@
  *
  *  Inverts a psImage matrix and returns the determinant as an option through the argument list. If the user
- *  specifies NULL as the outImage argument, then it will automatically be created. The input image must be 
- *  square. This function operates only with the psF64 data type. Input and output arguments should not be 
- *  the same. GSL indexes the top row as the zero row, not the bottom.  
+ *  specifies NULL as the outImage argument, then it will automatically be created. The input image must be
+ *  square. This function operates only with the psF64 data type. Input and output arguments should not be
+ *  the same. GSL indexes the top row as the zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to inverted psImage.
@@ -77 +77 @@
     psImage* outImage,                 ///< Image to return, or NULL for in-place substitution.
     const psImage* inImage,            ///< Image to be inverted
-    float *restrict det                ///< Determinant to return, or NULL
+    psF32 *det                         ///< Determinant to return, or NULL
 );
 
@@ -83 +83 @@
  *
  *  Calculates the determinant of a psImage matrix and returns the single precision floating point result. The
- *  input image must be square. This function operates only with the psF64 data type. GSL indexes the top row 
- *  as the zero row, not the bottom.  
+ *  input image must be square. This function operates only with the psF64 data type. GSL indexes the top row
+ *  as the zero row, not the bottom.
  *
  *  @return  float: Determinant from psImage.
  */
-float *psMatrixDeterminant(
-    const psImage* restrict inMatrix   ///< Image used to calculate determinant.
+psF32 *psMatrixDeterminant(
+    const psImage* inMatrix        ///< Image used to calculate determinant.
 );
 
 /** Performs psImage matrix multiplication.
  *
- *  Performs a classical matrix multiplication involving row and column operations. Input images must be square 
- *  and the same size. If the user specifies NULL as the outImage argument, then it will automatically be 
+ *  Performs a classical matrix multiplication involving row and column operations. Input images must be square
+ *  and the same size. If the user specifies NULL as the outImage argument, then it will automatically be
  *  created. This function operates only with the psF64 data type. GSL indexes the top row as the
- *  zero row, not the bottom.  
+ *  zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to resulting psImage.
@@ -109 +109 @@
 /** Transpose matrix.
  *
- *  Performs psImage matrix transpose by substituting existing rows for columns. The input image must be 
+ *  Performs psImage matrix transpose by substituting existing rows for columns. The input image must be
  *  square. If the user specifies NULL as the outImage argument, then it will automaticallty be created.
- *  This function operates only with the psF64 data type. GSL indexes the top row as the zero 
- *  row, not the bottom.  
+ *  This function operates only with the psF64 data type. GSL indexes the top row as the zero
+ *  row, not the bottom.
  *
  *  @return  psImage* : Pointer to transposed psImage.
@@ -123 +123 @@
 /** Calculate matrix eigenvectors.
  *
- *  Calculates the eigenvectors for a matrix. The input image must be symmetric and square. If the user 
- *  specifies NULL as the outImage argument, then it will automatically be created. This function operates 
- *  only with the psF64 data type. GSL indexes the top row as the zero row, not the bottom.  
+ *  Calculates the eigenvectors for a matrix. The input image must be symmetric and square. If the user
+ *  specifies NULL as the outImage argument, then it will automatically be created. This function operates
+ *  only with the psF64 data type. GSL indexes the top row as the zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to matrix of Eigenvectors.
@@ -152 +152 @@
  *  Converts a vector into a psImage matrix. If the dimensionality of the vector is PS_DIMEN_VECTOR, then the
  *  resulting psImage is a 1d column. If the dimensionality of the vector is PS_DIMEN_TRANSV, then the
- *  resulting psImage is a 1d row. If the user specifies NULL as the outImage argument,  then it will 
- *  automatically be created. This function operates only with the psF64 data type.  
+ *  resulting psImage is a 1d row. If the user specifies NULL as the outImage argument,  then it will
+ *  automatically be created. This function operates only with the psF64 data type.
  *
  *  @return  psVector* : Pointer to psIamge.

trunk/psLib/src/math/psMatrix.c

@@ -1 +1 @@
 /** @file  psMatrix.c
  *
- *  @brief Provides functions for linear algebra operations on psImages and psVectors. 
+ *  @brief Provides functions for linear algebra operations on psImages and psVectors.
  *
  *  Functions are provided to:
@@ -19 +19 @@
  *
  *  @author Ross Harman, MHPCC
- *   
- *  @version $Revision: 1.19 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-11-04 01:04:59 $
+ *
+ *  @version $Revision: 1.20 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-09 20:51:16 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -41 +41 @@
 #include "psMatrix.h"
 #include "psConstants.h"
-
-
+#include "psDataManipErrors.h"
+
+
+/*****************************************************************************/
+/* DEFINE STATEMENTS                                                         */
+/*****************************************************************************/
 
 /** Preprocessor macro to generate error for image dimensionality not set to PS_DIMEN_IMAGE */
@@ -50 +54 @@
             "Invalid operation. %s has incorrect dimensionality %d.", #NAME, PS_DIMEN);             \
     return RETURN;                                                                                  \
-} else if(NAME->type.type != PS_TYPE_F64) {                                                         \
+} else if(NAME->type.type!=PS_TYPE_F64 && NAME->type.type!=PS_TYPE_F32) {                           \
     psError(PS_ERR_BAD_PARAMETER_TYPE, true,                                                        \
             "Invalid operation. %s not PS_TYPE_F64.", #NAME);                                       \
@@ -67 +71 @@
 #define PS_CHECK_SQUARE(NAME, RETURN)                                                               \
 if (NAME->numCols != NAME->numRows) {                                                               \
-    psError(PS_ERR_BAD_PARAMETER_SIZE, true,                                                       \
-            "Invalid operation: %s not square array.", #NAME);                             \
+    psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: %s not square array.", #NAME);     \
     return RETURN;                                                                                  \
 }
@@ -79 +82 @@
 LHS_NAME.data  = RHS_NAME;
 
+
+/*****************************************************************************/
+/* FILE STATIC FUNCTIONS                                                     */
+/*****************************************************************************/
+
+static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector);
+static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector);
+static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage);
+static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix);
+
+/** Static function to copy psF32 or psF64 vector data to a GSL vector */
+static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector)
+{
+    psU32 i = 0;
+    psU32 n = 0;
+
+
+    n = inVector->n;
+    for(i=0; i<n; i++) {
+        if(inVector->type.type == PS_TYPE_F32) {
+            outGslVector->data[i] = (psS32)inVector->data.F32[i];
+        } else {
+            outGslVector->data[i] = inVector->data.F64[i];
+        }
+    }
+}
+
+/** Static function to copy GSL vector data to a psF32 or psF64 vector */
+static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector)
+{
+    psU32 i = 0;
+    psU32 n = 0;
+
+
+    n = outVector->n;
+    for(i=0; i<n; i++) {
+        if(outVector->type.type == PS_TYPE_F32) {
+            outVector->data.F32[i] = (psF32)inGslVector->data[i];
+        } else {
+            outVector->data.F64[i] = inGslVector->data[i];
+        }
+    }
+}
+
+/** Static function to copy psF32 or psF64 image data to a GSL matrix */
+static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage)
+{
+    psU32 i = 0;
+    psU32 j = 0;
+    psU32 numRows = 0;
+    psU32 numCols = 0;
+
+
+    numRows = inImage->numRows;
+    numCols = inImage->numCols;
+    for(i=0; i<numRows; i++) {
+        for(j=0; j<numCols; j++) {
+            if(inImage->type.type == PS_TYPE_F32) {
+                outGslMatrix->data[i*numCols+j] = (psS32)inImage->data.F32[i][j];
+            } else {
+                outGslMatrix->data[i*numCols+j] = inImage->data.F64[i][j];
+            }
+        }
+    }
+}
+
+/** Static function to copy GSL matrix data to a psF32 or psF64 image */
+static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix)
+{
+    psU32 i = 0;
+    psU32 j = 0;
+    psU32 numRows = 0;
+    psU32 numCols = 0;
+
+
+    numRows = outImage->numRows;
+    numCols = outImage->numCols;
+    for(i=0; i<numRows; i++) {
+        for(j=0; j<numCols; j++) {
+            if(outImage->type.type == PS_TYPE_F32) {
+                outImage->data.F32[i][j] = (psF32)inGslMatrix->data[i*numCols+j];
+            } else {
+                outImage->data.F64[i][j] = inGslMatrix->data[i*numCols+j];
+            }
+        }
+    }
+}
+
+
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
@@ -86 +178 @@
 {
     psS32 signum = 0;
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix lu;
+    gsl_matrix *lu = NULL;
     gsl_permutation perm;
 
-    #define psMatrixLUD_EXIT { \
-                               psFree(outImage); \
-                               return NULL; \
-                             }
-
-    // Error checks
+
+    #define psMatrixLUD_EXIT {psFree(outImage); return NULL;}
+
+    // Error checks
+    PS_IMAGE_CHECK_NULL_GENERAL(inImage, psMatrixLUD_EXIT);
     PS_CHECK_POINTERS(inImage, outImage, outImage);
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
-
-    PS_IMAGE_CHECK_NULL_GENERAL(inImage, psMatrixLUD_EXIT);
     PS_VECTOR_CHECK_NULL_GENERAL(outPerm, psMatrixLUD_EXIT);
-
     PS_CHECK_DIMEN_AND_TYPE(outPerm, PS_DIMEN_VECTOR, outImage);
-
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     psVectorRecycle(outPerm, inImage->numRows, inImage->type.type);
+    PS_CHECK_SQUARE(inImage, outImage);
+    PS_CHECK_SQUARE(outImage, outImage);
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
 
     // Initialize GSL data
@@ -118 +205 @@
     outPerm->n = numCols;
     perm.data = outPerm->data.V;
-    PS_GSL_MATRIX_INITIALIZE(lu, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    // Copy psImage input data into GSL matrix data to keep input data pristine
-    memcpy(lu.data, inImage->data.V[0], arraySize);
+    lu = gsl_matrix_alloc(numRows, numCols);
+
+    // Copy psImage data into GSL matrix data
+    psImageToGslMatrix(lu, inImage);
 
     // Calculate LU decomposition
-    gsl_linalg_LU_decomp(&lu, &perm, &signum);
+    gsl_linalg_LU_decomp(lu, &perm, &signum);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, lu);
+
+    // Free GSL data
+    gsl_matrix_free(lu);
 
     return outImage;
 }
 
-psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage,
-                          const psVector* inVector, const psVector* inPerm)
-{
-    psS32 arraySize = 0;
+psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage, const psVector* inVector,
+                          const psVector* inPerm)
+{
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix lu;
+    gsl_matrix *lu;
     gsl_permutation perm;
-    gsl_vector b;
-    gsl_vector x;
+    gsl_vector *b = NULL;
+    gsl_vector *x = NULL;
 
     // Error checks
@@ -157 +245 @@
     PS_VECTOR_CHECK_NULL(inPerm, outVector);
     PS_CHECK_DIMEN_AND_TYPE(inPerm, PS_DIMEN_VECTOR, outVector);
-
     psVectorRecycle(outVector, inImage->numRows, inImage->type.type);
-
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
 
     // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(lu, inImage->data.V[0]);
+    lu = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(lu, inImage);
+    b = gsl_vector_alloc(inVector->n);
+    psVectorToGslVector(b, inVector);
+    x = gsl_vector_alloc(inVector->n);
 
     outVector->n = numCols;
-
     perm.size = inPerm->n;
     perm.data = inPerm->data.V;
 
-    b.size = inVector->n;
-    b.stride = 1;
-    b.data = inVector->data.V;
-
-    x.size = numCols;
-    x.stride = 1;
-    x.data = outVector->data.V;
-
     // Solve for {x} in equation: {b} = [A]{x}
-    gsl_linalg_LU_solve(&lu, &perm, &b, &x);
+    gsl_linalg_LU_solve(lu, &perm, b, x);
+
+    // Copy GSL vector data to psVector data
+    gslVectorToPsVector(outVector, x);
+
+    // Free GSL data
+    gsl_vector_free(b);
+    gsl_vector_free(x);
+    gsl_matrix_free(lu);
 
     return outVector;
 }
 
-psImage* psMatrixInvert(psImage* outImage, const psImage* inImage, float *restrict det)
+psImage* psMatrixInvert(psImage* outImage, const psImage* inImage, psF32 *det)
 {
     psS32 signum = 0;
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix inv;
-    gsl_matrix *lu;
+    gsl_matrix *inv = NULL;
+    gsl_matrix *lu = NULL;
     gsl_permutation *perm = NULL;
 
@@ -204 +291 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     PS_IMAGE_CHECK_EMPTY(inImage, outImage);
-
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
+    PS_CHECK_SQUARE(inImage, outImage);
+    PS_CHECK_SQUARE(outImage, outImage);
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
-
-    // Allocate GSL structs
-    perm = gsl_permutation_alloc(inImage->numRows);
+
+    // Initialize GSL data
+    perm = gsl_permutation_alloc(numRows);
     lu = gsl_matrix_alloc(numRows, numCols);
-
-    // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(inv, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    // Copy psImage input data into GSL matrix data to keep input data pristine
-    memcpy(lu->data, inImage->data.V[0], arraySize);
+    inv = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(lu, inImage);
 
     // Invert data and calculate determinant
     gsl_linalg_LU_decomp(lu, perm, &signum);
-    gsl_linalg_LU_invert(lu, perm, &inv);
+    gsl_linalg_LU_invert(lu, perm, inv);
     *det = (float)gsl_linalg_LU_det(lu, signum);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, inv);
 
     // Free GSL structs
     gsl_permutation_free(perm);
     gsl_matrix_free(lu);
+    gsl_matrix_free(inv);
 
     return outImage;
 }
 
-float *psMatrixDeterminant(const psImage* restrict inImage)
+psF32 *psMatrixDeterminant(const psImage* inImage)
 {
     psS32 signum = 0;
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    float *det = NULL;
+    psF32 *det = NULL;
     gsl_matrix *lu = NULL;
     gsl_permutation *perm = NULL;
@@ -252 +334 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, NULL);
     PS_IMAGE_CHECK_EMPTY(inImage, NULL);
+    PS_CHECK_SQUARE(inImage, 0);
 
     // Initialize data
     numRows = inImage->numRows;
     numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
 
     // Allocate GSL structs
     perm = gsl_permutation_alloc(numRows);
     lu = gsl_matrix_alloc(numRows, numCols);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, 0);
-
-    // Copy psImage input data into GSL matrix data to keep input data pristine
-    memcpy(lu->data, inImage->data.V[0], arraySize);
+    psImageToGslMatrix(lu, inImage);
 
     // Calculate determinant
-    det = (float *)psAlloc(sizeof(float));
+    det = (psF32*)psAlloc(sizeof(psF32));
     gsl_linalg_LU_decomp(lu, perm, &signum);
-    *det = (float)gsl_linalg_LU_det(lu, signum);
+    *det = (psF32)gsl_linalg_LU_det(lu, signum);
 
     // Free GSL structs
@@ -282 +359 @@
 psImage* psMatrixMultiply(psImage* outImage, psImage* inImage1, psImage* inImage2)
 {
-    psS32 arraySize = 0;
     psS32 numRows = 0;
     psS32 numCols = 0;
-    gsl_matrix m1;
-    gsl_matrix m2;
-    gsl_matrix m3;
+    gsl_matrix *m1 = NULL;
+    gsl_matrix *m2 = NULL;
+    gsl_matrix *m3 = NULL;
 
     // Error checks
@@ -299 +375 @@
     PS_IMAGE_CHECK_EMPTY(inImage2, outImage);
     PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, outImage);
-
     psImageRecycle(outImage, inImage2->numCols, inImage2->numRows, inImage2->type.type);
+    PS_CHECK_SQUARE(inImage1, outImage);
+    PS_CHECK_SQUARE(inImage2, outImage);
+    PS_CHECK_SQUARE(outImage, outImage);
 
     // Initialize data
     numRows = inImage1->numRows;
     numCols = inImage1->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(outImage->type.type) * numRows * numCols;
 
     // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(m1, inImage1->data.V[0]);
-    PS_GSL_MATRIX_INITIALIZE(m2, inImage2->data.V[0]);
-    PS_GSL_MATRIX_INITIALIZE(m3, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage1, outImage);
-    PS_CHECK_SQUARE(inImage2, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    gsl_linalg_matmult(&m1, &m2, &m3);
+    m1 = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(m1, inImage1);
+    m2 = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(m2, inImage2);
+    m3 = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(m3, outImage);
+
+    // Perform multiplication
+    gsl_linalg_matmult(m1, m2, m3);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, m3);
+
+    // Free GSL structs
+    gsl_matrix_free(m1);
+    gsl_matrix_free(m2);
+    gsl_matrix_free(m3);
 
     return outImage;
@@ -324 +408 @@
 psImage* psMatrixTranspose(psImage* outImage, const psImage* inImage)
 {
-    psS32 arraySize = 0;
-    psS32 numRows = 0;
-    psS32 numCols = 0;
-    gsl_matrix trans;
+    psU32 i = 0;
+    psU32 j = 0;
+    psS32 numRowsIn = 0;
+    psS32 numColsIn = 0;
+    psS32 numRowsOut = 0;
+    psS32 numColsOut = 0;
+
 
     // Error checks
@@ -334 +421 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     PS_IMAGE_CHECK_EMPTY(inImage, outImage);
-
-    psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
+    outImage = psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
 
     // Initialize data
-    numRows = inImage->numRows;
-    numCols = inImage->numCols;
-    arraySize = PSELEMTYPE_SIZEOF(inImage->type.type) * numRows * numCols;
-
-    // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(trans, outImage->data.V[0]);
-
-    // Non-square matrices not allowed
-    PS_CHECK_SQUARE(inImage, outImage);
-    PS_CHECK_SQUARE(outImage, outImage);
-
-    // Copy psImage input data into psImage output data to keep input data pristine
-    memcpy(outImage->data.V[0], inImage->data.V[0], arraySize);
-
-    // Transpose data
-    gsl_matrix_transpose(&trans);
+    numRowsIn = inImage->numRows;
+    numColsIn = inImage->numCols;
+    numRowsOut = outImage->numRows;
+    numColsOut = outImage->numCols;
+
+    if(numRowsIn!=numColsOut && numRowsOut!=numColsIn) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, PS_ERRORTEXT_psMatrix_TRANSPOSE_MISMATCH);
+    }
+
+    if(outImage->type.type == PS_TYPE_F32) {
+        for(i=0; i<numRowsOut; i++) {
+            for(j=0; j<numColsOut; j++) {
+                outImage->data.F32[i][j] = inImage->data.F32[j][i];
+            }
+        }
+    } else {
+        for(i=0; i<numRowsOut; i++) {
+            for(j=0; j<numColsOut; j++) {
+                outImage->data.F64[i][j] = inImage->data.F64[j][i];
+            }
+        }
+    }
 
     return outImage;
@@ -364 +456 @@
     gsl_vector *eVals = NULL;
     gsl_eigen_symmv_workspace *w = NULL;
-    gsl_matrix out;
-    gsl_matrix in;
+    gsl_matrix *out = NULL;
+    gsl_matrix *in = NULL;
+
 
     // Error checks
@@ -372 +465 @@
     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, outImage);
     PS_IMAGE_CHECK_EMPTY(inImage, outImage);
-
     psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
 
@@ -379 +471 @@
     numCols = inImage->numCols;
 
-    // Initialize GSL data
-    PS_GSL_MATRIX_INITIALIZE(in, inImage->data.V[0]);
-    PS_GSL_MATRIX_INITIALIZE(out, outImage->data.V[0]);
+    in = gsl_matrix_alloc(numRows, numCols);
+    psImageToGslMatrix(in, inImage);
+    out = gsl_matrix_alloc(numRows, numCols);
 
     // Allocate GSL structs
@@ -392 +484 @@
 
     // Calculate Eigenvalues and Eigenvectors...Eigenvalues not currently used
-    gsl_eigen_symmv(&in, eVals, &out, w);
+    gsl_eigen_symmv(in, eVals, out, w);
+
+    // Copy GSL matrix data to psImage data
+    gslMatrixToPsImage(outImage, out);
 
     // Free GSL structs
+    gsl_matrix_free(in);
+    gsl_matrix_free(out);
     gsl_eigen_symmv_free(w);
     gsl_vector_free(eVals);
@@ -405 +502 @@
     psS32 size = 0;
 
-    #define psMatrixToVector_EXIT { \
-                                    psFree(outVector); \
-                                    return NULL; \
-                                  }
+    #define psMatrixToVector_EXIT {psFree(outVector); return NULL;}
 
     // Error checks

trunk/psLib/src/math/psMatrix.h

@@ -2 +2 @@
 /** @file  psMatrix.h
  *
- *  @brief Provides functions for linear algebra operations on psImages and psVectors. 
+ *  @brief Provides functions for linear algebra operations on psImages and psVectors.
  *
  *  Functions are provided to:
@@ -22 +22 @@
  *  @author Ross Harman, MHPCC
  *
- *  @version $Revision: 1.12 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-21 23:44:10 $
+ *  @version $Revision: 1.13 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-12-09 20:51:22 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -37 +37 @@
  *
  *  Performs a LU decomposition on a psImage matrix and returns the LU matrix. If the user specifies NULL for
- *  the outImage or outPerm arguments, then they will be automatically created. The input image must 
- *  be square. This function operates only with the psF64 data type. Input and output arguments should not be 
- *  the same. GSL indexes the top row as the zero row, not the bottom.  
+ *  the outImage or outPerm arguments, then they will be automatically created. The input image must
+ *  be square. This function operates only with the psF64 data type. Input and output arguments should not be
+ *  the same. GSL indexes the top row as the zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to LU decomposed psImage.
@@ -51 +51 @@
 /** LU Solution of psImage matrix.
  *
- *  Solves for and returns the psVector, {x} in the equation [A]{x} = {b}. If the user specifies NULL as the 
- *  outVector argument, then it will automatically be created. The input image must be square. This function 
- *  operates only with the psF64 data type. Input and output arguments should not be the same. GSL indexes 
- *  the top row as the zero row, not the bottom.  
+ *  Solves for and returns the psVector, {x} in the equation [A]{x} = {b}. If the user specifies NULL as the
+ *  outVector argument, then it will automatically be created. The input image must be square. This function
+ *  operates only with the psF64 data type. Input and output arguments should not be the same. GSL indexes
+ *  the top row as the zero row, not the bottom.
  *
  *  @return  psVector* : Pointer to psVector solution of matrix equation.
@@ -68 +68 @@
  *
  *  Inverts a psImage matrix and returns the determinant as an option through the argument list. If the user
- *  specifies NULL as the outImage argument, then it will automatically be created. The input image must be 
- *  square. This function operates only with the psF64 data type. Input and output arguments should not be 
- *  the same. GSL indexes the top row as the zero row, not the bottom.  
+ *  specifies NULL as the outImage argument, then it will automatically be created. The input image must be
+ *  square. This function operates only with the psF64 data type. Input and output arguments should not be
+ *  the same. GSL indexes the top row as the zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to inverted psImage.
@@ -77 +77 @@
     psImage* outImage,                 ///< Image to return, or NULL for in-place substitution.
     const psImage* inImage,            ///< Image to be inverted
-    float *restrict det                ///< Determinant to return, or NULL
+    psF32 *det                         ///< Determinant to return, or NULL
 );
 
@@ -83 +83 @@
  *
  *  Calculates the determinant of a psImage matrix and returns the single precision floating point result. The
- *  input image must be square. This function operates only with the psF64 data type. GSL indexes the top row 
- *  as the zero row, not the bottom.  
+ *  input image must be square. This function operates only with the psF64 data type. GSL indexes the top row
+ *  as the zero row, not the bottom.
  *
  *  @return  float: Determinant from psImage.
  */
-float *psMatrixDeterminant(
-    const psImage* restrict inMatrix   ///< Image used to calculate determinant.
+psF32 *psMatrixDeterminant(
+    const psImage* inMatrix        ///< Image used to calculate determinant.
 );
 
 /** Performs psImage matrix multiplication.
  *
- *  Performs a classical matrix multiplication involving row and column operations. Input images must be square 
- *  and the same size. If the user specifies NULL as the outImage argument, then it will automatically be 
+ *  Performs a classical matrix multiplication involving row and column operations. Input images must be square
+ *  and the same size. If the user specifies NULL as the outImage argument, then it will automatically be
  *  created. This function operates only with the psF64 data type. GSL indexes the top row as the
- *  zero row, not the bottom.  
+ *  zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to resulting psImage.
@@ -109 +109 @@
 /** Transpose matrix.
  *
- *  Performs psImage matrix transpose by substituting existing rows for columns. The input image must be 
+ *  Performs psImage matrix transpose by substituting existing rows for columns. The input image must be
  *  square. If the user specifies NULL as the outImage argument, then it will automaticallty be created.
- *  This function operates only with the psF64 data type. GSL indexes the top row as the zero 
- *  row, not the bottom.  
+ *  This function operates only with the psF64 data type. GSL indexes the top row as the zero
+ *  row, not the bottom.
  *
  *  @return  psImage* : Pointer to transposed psImage.
@@ -123 +123 @@
 /** Calculate matrix eigenvectors.
  *
- *  Calculates the eigenvectors for a matrix. The input image must be symmetric and square. If the user 
- *  specifies NULL as the outImage argument, then it will automatically be created. This function operates 
- *  only with the psF64 data type. GSL indexes the top row as the zero row, not the bottom.  
+ *  Calculates the eigenvectors for a matrix. The input image must be symmetric and square. If the user
+ *  specifies NULL as the outImage argument, then it will automatically be created. This function operates
+ *  only with the psF64 data type. GSL indexes the top row as the zero row, not the bottom.
  *
  *  @return  psImage* : Pointer to matrix of Eigenvectors.
@@ -152 +152 @@
  *  Converts a vector into a psImage matrix. If the dimensionality of the vector is PS_DIMEN_VECTOR, then the
  *  resulting psImage is a 1d column. If the dimensionality of the vector is PS_DIMEN_TRANSV, then the
- *  resulting psImage is a 1d row. If the user specifies NULL as the outImage argument,  then it will 
- *  automatically be created. This function operates only with the psF64 data type.  
+ *  resulting psImage is a 1d row. If the user specifies NULL as the outImage argument,  then it will
+ *  automatically be created. This function operates only with the psF64 data type.
  *
  *  @return  psVector* : Pointer to psIamge.