Index: /trunk/psLib/src/dataManip/Makefile =================================================================== --- /trunk/psLib/src/dataManip/Makefile (revision 760) +++ /trunk/psLib/src/dataManip/Makefile (revision 761) @@ -6,5 +6,9 @@ include ../Makefile.Globals CFLAGS := $(CFLAGS_RELOC) -I. -I../sysUtils -I../collections -I.. -SRC_OBJS = psStats.o psFunctions.o psMinimize.o + +SRC_OBJS = psStats.o \ + psFunctions.o \ + psMinimize.o \ + psMatrix.o all: $(TARGET_STATIC) Index: /trunk/psLib/src/dataManip/psMatrix.c =================================================================== --- /trunk/psLib/src/dataManip/psMatrix.c (revision 761) +++ /trunk/psLib/src/dataManip/psMatrix.c (revision 761) @@ -0,0 +1,374 @@ +/** @file psMatrix.c + * + * @brief Provides functions for linear algebra operations on psImages and psVectors. + * + * Functions are provided to: + * Transpose a psImage + * Compute LUD + * Solve LUD + * Matrix inversion + * Calculate determinant + * Matrix addition + * Matrix subtraction + * Matrix multiplication + * Calculate Eigenvectors + * Convert matrix to vector + * Convert vector to matrix + * + * These functions treat psImages as if they were matrices, therefore there is no psMatrix. + * + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-05-24 21:09:39 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include +#include +#include + +#include "psError.h" +#include "psImage.h" +#include "psVector.h" +#include "psMatrix.h" + +/******************************************************************************/ +/* DEFINE STATEMENTS */ +/******************************************************************************/ + +// None + +/******************************************************************************/ +/* TYPE DEFINITIONS */ +/******************************************************************************/ + +// None + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +#define PS_MATRIX_TRANSPOSE(PS_TYPE) \ +{ \ + int arraySize = 0; \ + int numRows = 0; \ + int numCols = 0; \ + gsl_matrix##PS_TYPE trans; \ + \ + /* Initialize data */ \ + numRows = inImage->numRows; \ + numCols = inImage->numCols; \ + arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; \ + \ + /* Copy psImage input data into psImage output data to keep input data pristine */ \ + memcpy(outImage->data.v[0], inImage->data.v[0], arraySize); \ + \ + /* Manually fill inverted output matrix so it will be aligned with output image */ \ + trans.size1 = numRows; \ + trans.size2 = numCols; \ + trans.tda = numCols; \ + trans.data = outImage->data.v[0]; \ + \ + /* Transpose data */ \ + gsl_matrix##PS_TYPE##_transpose(&trans); \ +} + +psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage) +{ + psElemType elemType = 0; + + elemType = inImage->type.type; + switch(elemType) { + case PS_TYPE_FLOAT: + PS_MATRIX_TRANSPOSE(_float); + break; + case PS_TYPE_DOUBLE: + PS_MATRIX_TRANSPOSE(); + break; + default: + psError(__func__, " : Line %d - Invalid psElemType: %d\n", __LINE__, elemType); + } + + return outImage; +} + +psImage* psMatrixOp(psImage *outImage, psImage *inImage1, const char op, psImage *inImage2) +{ + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix m2; + gsl_matrix m3; + + // Initialize data + numRows = inImage1->numRows; + numCols = inImage1->numCols; + arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols; + + m2.size1 = numRows; + m2.size2 = numCols; + m2.tda = numCols; + m2.data = inImage2->data.v[0]; + + m3.size1 = numRows; + m3.size2 = numCols; + m3.tda = numCols; + m3.data = outImage->data.v[0]; + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(m3.data, inImage1->data.v[0], arraySize); + + + switch(op) { + case '+': + gsl_matrix_add(&m3, &m2); + break; + case '-': + gsl_matrix_sub(&m3, &m2); + break; + case '*': + gsl_linalg_matmult(&m3, &m2, &m3); + break; + default: + psError(__func__, " : Line %d - Invalid psMatrixOp operation: %c\n", __LINE__, op); + \ + } + + return NULL; +} + +psImage *psMatrixLUD(psImage *outImage, psVector *outPerm, psImage *inImage) +{ + int signum = 0; + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix lu; + gsl_permutation perm; + + // Initialize data + numRows = inImage->numRows; + numCols = inImage->numCols; + arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; + + // Manually fill GSL structs so they will be aligned with output data + perm.size = numCols; + outPerm->n = numCols; + perm.data = outPerm->vec.ui32; + lu.size1 = numRows; + lu.size2 = numCols; + lu.tda = numCols; + lu.data = outImage->data.v[0]; + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(lu.data, inImage->data.v[0], arraySize); + + // Calculate LU decomposition + gsl_linalg_LU_decomp(&lu, &perm, &signum); + + return outImage; +} + +psVector *psMatrixLUSolve(psVector *outVector, const psImage *luImage, const psVector *inVector, const psVector *inPerm) +{ + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix lu; + gsl_permutation perm; + gsl_vector b; + gsl_vector x; + + // Initialize data + numRows = luImage->numRows; + numCols = luImage->numCols; + arraySize = PSELEMTYPE_SIZEOF(luImage->type.type)*numRows*numCols; + + // Manually fill GSL structs so they will be aligned with output data + lu.size1 = luImage->numRows; + lu.size2 = luImage->numCols; + lu.tda = luImage->numCols; + lu.data = luImage->data.v[0]; + + outVector->n = numCols; + + perm.size = inPerm->n; + perm.data = inPerm->vec.ui32; + + b.size = inVector->n; + b.stride = 1; + b.data = inVector->vec.v; + + x.size = numCols; + x.stride = 1; + x.data = outVector->vec.v; + + // Solve for {x} in equation: {b} = [A]{x} + gsl_linalg_LU_solve(&lu, &perm, &b, &x); + + return outVector; +} + +psImage *psMatrixInvert(psImage *outImage, const psImage *inImage, float *restrict det) +{ + int signum = 0; + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix inv; + gsl_matrix *lu; + gsl_permutation *perm = NULL; + + // 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); + lu = gsl_matrix_alloc(numRows, numCols); + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(lu->data, inImage->data.v[0], arraySize); + + // Manually fill inverted output matrix so it will be aligned with output image + inv.size1 = outImage->numRows; + inv.size2 = outImage->numCols; + inv.tda = outImage->numCols; + inv.data = outImage->data.v[0]; + + // Invert data and calculate determinant + gsl_linalg_LU_decomp(lu, perm, &signum); + gsl_linalg_LU_invert(lu, perm, &inv); + *det = (float)gsl_linalg_LU_det(lu, signum); + + // Free GSL structs + gsl_permutation_free(perm); + gsl_matrix_free(lu); + + return outImage; +} + +float psMatrixDeterminant(const psImage *restrict inImage) +{ + int signum = 0; + int arraySize = 0; + int numRows = 0; + int numCols = 0; + float det = 0.0f; + gsl_matrix *lu = NULL; + gsl_permutation *perm = NULL; + + // 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); + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(lu->data, inImage->data.v[0], arraySize); + + // Calculate determinant + gsl_linalg_LU_decomp(lu, perm, &signum); + det = (float)gsl_linalg_LU_det(lu, signum); + + // Free GSL structs + gsl_permutation_free(perm); + gsl_matrix_free(lu); + + return det; +} + +psImage *psMatrixEigenvectors(psImage *outImage, psImage *inImage) +{ + int numRows = 0; + int numCols = 0; + gsl_vector *eVals = NULL; + gsl_eigen_symmv_workspace *w = NULL; + gsl_matrix out; + gsl_matrix in; + + // Initialize data + numRows = inImage->numRows; + numCols = inImage->numCols; + out.data = outImage->data.v[0]; + + // Manually fill GSL structs so they will be aligned with output data + in.size1 = numRows; + in.size2 = numCols; + in.tda = numCols; + in.data = inImage->data.v[0]; + + out.size1 = numRows; + out.size2 = numCols; + out.tda = numCols; + out.data = outImage->data.v[0]; + + // Allocate GSL structs + eVals = gsl_vector_alloc(numRows); + w = gsl_eigen_symmv_alloc(numRows); + + // Calculate Eigenvalues and Eigenvectors. Eigenvalues not currently used + gsl_eigen_symmv(&in, eVals, &out, w); + + // Free GSL structs + gsl_eigen_symmv_free(w); + gsl_vector_free(eVals); + + return outImage; +} + +psVector *psMatrixToVector(psVector *outVector, psImage *inImage) +{ + // if inimage isn't 1d + // if inimage nrows != vector len + + int colSize = 0; + + colSize = PSELEMTYPE_SIZEOF(inImage->type.type)*inImage->numRows; + memcpy(outVector->vec.v, inImage->data.v[0], colSize); + + return outVector; +} + +psImage *psVectorToMatrix(psImage *outImage, psVector *inVector) +{ + // if inimage isn't 1d + // if inimage nrows != vector len + + int colSize = 0; + + colSize = PSELEMTYPE_SIZEOF(outImage->type.type)*outImage->numRows; + memcpy(outImage->data.v[0], inVector->vec.v, colSize); + + return outImage; +} + Index: /trunk/psLib/src/math/psMatrix.c =================================================================== --- /trunk/psLib/src/math/psMatrix.c (revision 761) +++ /trunk/psLib/src/math/psMatrix.c (revision 761) @@ -0,0 +1,374 @@ +/** @file psMatrix.c + * + * @brief Provides functions for linear algebra operations on psImages and psVectors. + * + * Functions are provided to: + * Transpose a psImage + * Compute LUD + * Solve LUD + * Matrix inversion + * Calculate determinant + * Matrix addition + * Matrix subtraction + * Matrix multiplication + * Calculate Eigenvectors + * Convert matrix to vector + * Convert vector to matrix + * + * These functions treat psImages as if they were matrices, therefore there is no psMatrix. + * + * @author Ross Harman, MHPCC + * + * @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-05-24 21:09:39 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +/******************************************************************************/ +/* INCLUDE FILES */ +/******************************************************************************/ +#include +#include +#include +#include +#include + +#include "psError.h" +#include "psImage.h" +#include "psVector.h" +#include "psMatrix.h" + +/******************************************************************************/ +/* DEFINE STATEMENTS */ +/******************************************************************************/ + +// None + +/******************************************************************************/ +/* TYPE DEFINITIONS */ +/******************************************************************************/ + +// None + +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + +#define PS_MATRIX_TRANSPOSE(PS_TYPE) \ +{ \ + int arraySize = 0; \ + int numRows = 0; \ + int numCols = 0; \ + gsl_matrix##PS_TYPE trans; \ + \ + /* Initialize data */ \ + numRows = inImage->numRows; \ + numCols = inImage->numCols; \ + arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; \ + \ + /* Copy psImage input data into psImage output data to keep input data pristine */ \ + memcpy(outImage->data.v[0], inImage->data.v[0], arraySize); \ + \ + /* Manually fill inverted output matrix so it will be aligned with output image */ \ + trans.size1 = numRows; \ + trans.size2 = numCols; \ + trans.tda = numCols; \ + trans.data = outImage->data.v[0]; \ + \ + /* Transpose data */ \ + gsl_matrix##PS_TYPE##_transpose(&trans); \ +} + +psImage* psMatrixTranspose(psImage *outImage, const psImage *inImage) +{ + psElemType elemType = 0; + + elemType = inImage->type.type; + switch(elemType) { + case PS_TYPE_FLOAT: + PS_MATRIX_TRANSPOSE(_float); + break; + case PS_TYPE_DOUBLE: + PS_MATRIX_TRANSPOSE(); + break; + default: + psError(__func__, " : Line %d - Invalid psElemType: %d\n", __LINE__, elemType); + } + + return outImage; +} + +psImage* psMatrixOp(psImage *outImage, psImage *inImage1, const char op, psImage *inImage2) +{ + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix m2; + gsl_matrix m3; + + // Initialize data + numRows = inImage1->numRows; + numCols = inImage1->numCols; + arraySize = PSELEMTYPE_SIZEOF(outImage->type.type)*numRows*numCols; + + m2.size1 = numRows; + m2.size2 = numCols; + m2.tda = numCols; + m2.data = inImage2->data.v[0]; + + m3.size1 = numRows; + m3.size2 = numCols; + m3.tda = numCols; + m3.data = outImage->data.v[0]; + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(m3.data, inImage1->data.v[0], arraySize); + + + switch(op) { + case '+': + gsl_matrix_add(&m3, &m2); + break; + case '-': + gsl_matrix_sub(&m3, &m2); + break; + case '*': + gsl_linalg_matmult(&m3, &m2, &m3); + break; + default: + psError(__func__, " : Line %d - Invalid psMatrixOp operation: %c\n", __LINE__, op); + \ + } + + return NULL; +} + +psImage *psMatrixLUD(psImage *outImage, psVector *outPerm, psImage *inImage) +{ + int signum = 0; + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix lu; + gsl_permutation perm; + + // Initialize data + numRows = inImage->numRows; + numCols = inImage->numCols; + arraySize = PSELEMTYPE_SIZEOF(inImage->type.type)*numRows*numCols; + + // Manually fill GSL structs so they will be aligned with output data + perm.size = numCols; + outPerm->n = numCols; + perm.data = outPerm->vec.ui32; + lu.size1 = numRows; + lu.size2 = numCols; + lu.tda = numCols; + lu.data = outImage->data.v[0]; + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(lu.data, inImage->data.v[0], arraySize); + + // Calculate LU decomposition + gsl_linalg_LU_decomp(&lu, &perm, &signum); + + return outImage; +} + +psVector *psMatrixLUSolve(psVector *outVector, const psImage *luImage, const psVector *inVector, const psVector *inPerm) +{ + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix lu; + gsl_permutation perm; + gsl_vector b; + gsl_vector x; + + // Initialize data + numRows = luImage->numRows; + numCols = luImage->numCols; + arraySize = PSELEMTYPE_SIZEOF(luImage->type.type)*numRows*numCols; + + // Manually fill GSL structs so they will be aligned with output data + lu.size1 = luImage->numRows; + lu.size2 = luImage->numCols; + lu.tda = luImage->numCols; + lu.data = luImage->data.v[0]; + + outVector->n = numCols; + + perm.size = inPerm->n; + perm.data = inPerm->vec.ui32; + + b.size = inVector->n; + b.stride = 1; + b.data = inVector->vec.v; + + x.size = numCols; + x.stride = 1; + x.data = outVector->vec.v; + + // Solve for {x} in equation: {b} = [A]{x} + gsl_linalg_LU_solve(&lu, &perm, &b, &x); + + return outVector; +} + +psImage *psMatrixInvert(psImage *outImage, const psImage *inImage, float *restrict det) +{ + int signum = 0; + int arraySize = 0; + int numRows = 0; + int numCols = 0; + gsl_matrix inv; + gsl_matrix *lu; + gsl_permutation *perm = NULL; + + // 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); + lu = gsl_matrix_alloc(numRows, numCols); + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(lu->data, inImage->data.v[0], arraySize); + + // Manually fill inverted output matrix so it will be aligned with output image + inv.size1 = outImage->numRows; + inv.size2 = outImage->numCols; + inv.tda = outImage->numCols; + inv.data = outImage->data.v[0]; + + // Invert data and calculate determinant + gsl_linalg_LU_decomp(lu, perm, &signum); + gsl_linalg_LU_invert(lu, perm, &inv); + *det = (float)gsl_linalg_LU_det(lu, signum); + + // Free GSL structs + gsl_permutation_free(perm); + gsl_matrix_free(lu); + + return outImage; +} + +float psMatrixDeterminant(const psImage *restrict inImage) +{ + int signum = 0; + int arraySize = 0; + int numRows = 0; + int numCols = 0; + float det = 0.0f; + gsl_matrix *lu = NULL; + gsl_permutation *perm = NULL; + + // 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); + + // Copy psImage input data into GSL matrix data to keep input data pristine + memcpy(lu->data, inImage->data.v[0], arraySize); + + // Calculate determinant + gsl_linalg_LU_decomp(lu, perm, &signum); + det = (float)gsl_linalg_LU_det(lu, signum); + + // Free GSL structs + gsl_permutation_free(perm); + gsl_matrix_free(lu); + + return det; +} + +psImage *psMatrixEigenvectors(psImage *outImage, psImage *inImage) +{ + int numRows = 0; + int numCols = 0; + gsl_vector *eVals = NULL; + gsl_eigen_symmv_workspace *w = NULL; + gsl_matrix out; + gsl_matrix in; + + // Initialize data + numRows = inImage->numRows; + numCols = inImage->numCols; + out.data = outImage->data.v[0]; + + // Manually fill GSL structs so they will be aligned with output data + in.size1 = numRows; + in.size2 = numCols; + in.tda = numCols; + in.data = inImage->data.v[0]; + + out.size1 = numRows; + out.size2 = numCols; + out.tda = numCols; + out.data = outImage->data.v[0]; + + // Allocate GSL structs + eVals = gsl_vector_alloc(numRows); + w = gsl_eigen_symmv_alloc(numRows); + + // Calculate Eigenvalues and Eigenvectors. Eigenvalues not currently used + gsl_eigen_symmv(&in, eVals, &out, w); + + // Free GSL structs + gsl_eigen_symmv_free(w); + gsl_vector_free(eVals); + + return outImage; +} + +psVector *psMatrixToVector(psVector *outVector, psImage *inImage) +{ + // if inimage isn't 1d + // if inimage nrows != vector len + + int colSize = 0; + + colSize = PSELEMTYPE_SIZEOF(inImage->type.type)*inImage->numRows; + memcpy(outVector->vec.v, inImage->data.v[0], colSize); + + return outVector; +} + +psImage *psVectorToMatrix(psImage *outImage, psVector *inVector) +{ + // if inimage isn't 1d + // if inimage nrows != vector len + + int colSize = 0; + + colSize = PSELEMTYPE_SIZEOF(outImage->type.type)*outImage->numRows; + memcpy(outImage->data.v[0], inVector->vec.v, colSize); + + return outImage; +} +