IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Oct 16, 2009, 3:45:34 PM (17 years ago)
Author:
Paul Price
Message:

Cleaning (optimising) the gsl-psLib conversions using memcpy where appropriate. Adding function to solve matrix equation using SVD.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/pap/psLib/src/math/psMatrix.c

    r25842 r25862  
    9494LHS_NAME.data  = RHS_NAME;
    9595
    96 
    97 /*****************************************************************************/
    98 /* FILE STATIC FUNCTIONS                                                     */
    99 /*****************************************************************************/
    100 
    101 static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector);
    102 static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector);
    103 static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage);
    104 static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix);
     96////////////////////////////////////////////////////////////////////////////////
     97// Conversion functions
     98////////////////////////////////////////////////////////////////////////////////
     99
     100// gsl_vector holds *doubles*, so we can directly copy F64, but need to convert F32
    105101
    106102/** Static function to copy psF32 or psF64 vector data to a GSL vector */
    107 static void  psVectorToGslVector(gsl_vector *outGslVector,
    108                                  const psVector *inVector)
    109 {
    110     psU32 i = 0;
    111     psU32 n = 0;
    112 
    113 
    114     n = inVector->n;
    115     for(i=0; i<n; i++) {
    116         if(inVector->type.type == PS_TYPE_F32) {
    117             outGslVector->data[i] = (psF64)inVector->data.F32[i];
     103static void vectorPStoGSL(gsl_vector *out, const psVector *in)
     104{
     105    long n = in->n;                     // Size of input
     106    switch (in->type.type) {
     107      case PS_TYPE_F32:
     108        for (long i = 0; i < n; i++) {
     109            out->data[i] = in->data.F32[i];
     110        }
     111        break;
     112      case PS_TYPE_F64:
     113        memcpy(out->data, in->data.F64, n * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     114        break;
     115      default:
     116        psAbort("Unsupported vector type: %x\n", in->type.type);
     117    }
     118    return;
     119}
     120
     121/** Static function to copy GSL vector data to a psF32 or psF64 vector */
     122static void vectorGSLtoPS(psVector *out, const gsl_vector *in)
     123{
     124    long n = out->n;                    // Size of output
     125    switch (out->type.type) {
     126      case PS_TYPE_F32:
     127        for (long i = 0; i < n; i++) {
     128            out->data.F32[i] = in->data[i];
     129        }
     130        break;
     131      case PS_TYPE_F64:
     132        memcpy(out->data.F64, in->data, n * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     133        break;
     134      default:
     135        psAbort("Unsupported vector type: %x\n", out->type.type);
     136    }
     137    return;
     138}
     139
     140
     141/** Static function to copy psF32 or psF64 image data to a GSL matrix */
     142static void matrixPStoGSL(gsl_matrix *out, const psImage *in)
     143{
     144    int numCols = in->numCols, numRows = in->numRows; // Size of matrix
     145    switch (in->type.type) {
     146      case PS_TYPE_F32:
     147        for (int y = 0, i = 0; y < numRows; y++) {
     148            for (int x = 0; x < numCols; x++, i++) {
     149                out->data[i] = in->data.F32[y][x];
     150            }
     151        }
     152        break;
     153      case PS_TYPE_F64:
     154        if (in->parent) {
     155            for (int y = 0, i = 0; y < numRows; y++, i += numCols) {
     156                memcpy(&out->data[i], in->data.F64[y], numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     157            }
    118158        } else {
    119             outGslVector->data[i] = inVector->data.F64[i];
    120         }
    121     }
    122 }
    123 
    124 /** Static function to copy GSL vector data to a psF32 or psF64 vector */
    125 static void gslVectorToPsVector(psVector *outVector,
    126                                 gsl_vector *inGslVector)
    127 {
    128     psU32 i = 0;
    129     psU32 n = 0;
    130 
    131 
    132     n = outVector->n;
    133     for(i=0; i<n; i++) {
    134         if(outVector->type.type == PS_TYPE_F32) {
    135             outVector->data.F32[i] = (psF32)inGslVector->data[i];
     159            memcpy(out->data, in->data.F64, numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     160        }
     161        break;
     162      default:
     163        psAbort("Unsupported vector type: %x\n", in->type.type);
     164    }
     165    return;
     166}
     167
     168/** Static function to copy GSL matrix data to a psF32 or psF64 image */
     169static void matrixGSLtoPS(psImage *out, gsl_matrix *in)
     170{
     171    int numCols = out->numCols, numRows = out->numRows; // Size of matrix
     172    switch (out->type.type) {
     173      case PS_TYPE_F32:
     174        for (int y = 0, i = 0; y < numRows; y++) {
     175            for (int x = 0; x < numCols; x++, i++) {
     176                out->data.F32[y][x] = in->data[i];
     177            }
     178        }
     179        break;
     180      case PS_TYPE_F64:
     181        if (out->parent) {
     182            for (int y = 0, i = 0; y < numRows; y++, i += numCols) {
     183                memcpy(out->data.F64[y], &in->data[i], numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     184            }
    136185        } else {
    137             outVector->data.F64[i] = inGslVector->data[i];
    138         }
    139     }
    140 }
    141 
    142 /** Static function to copy psF32 or psF64 image data to a GSL matrix */
    143 static void  psImageToGslMatrix(gsl_matrix *outGslMatrix,
    144                                 const psImage *inImage)
    145 {
    146     psU32 i = 0;
    147     psU32 j = 0;
    148     psU32 numRows = 0;
    149     psU32 numCols = 0;
    150 
    151 
    152     numRows = inImage->numRows;
    153     numCols = inImage->numCols;
    154     if(inImage->type.type == PS_TYPE_F32) {
    155         for(i=0; i<numRows; i++) {
    156             for(j=0; j<numCols; j++) {
    157                 outGslMatrix->data[i*numCols+j] = inImage->data.F32[i][j];
    158             }
    159         }
    160     } else {
    161         for(i=0; i<numRows; i++) {
    162             for(j=0; j<numCols; j++) {
    163                 outGslMatrix->data[i*numCols+j] = inImage->data.F64[i][j];
    164             }
    165         }
    166     }
    167 }
    168 
    169 /** Static function to copy GSL matrix data to a psF32 or psF64 image */
    170 static void gslMatrixToPsImage(psImage *outImage,
    171                                gsl_matrix *inGslMatrix)
    172 {
    173     psU32 i = 0;
    174     psU32 j = 0;
    175     psU32 numRows = 0;
    176     psU32 numCols = 0;
    177 
    178 
    179     numRows = outImage->numRows;
    180     numCols = outImage->numCols;
    181     if(outImage->type.type == PS_TYPE_F32) {
    182         for(i=0; i<numRows; i++) {
    183             for(j=0; j<numCols; j++) {
    184                 outImage->data.F32[i][j] = inGslMatrix->data[i*numCols+j];
    185             }
    186         }
    187     } else {
    188         for(i=0; i<numRows; i++) {
    189             for(j=0; j<numCols; j++) {
    190                 outImage->data.F64[i][j] = inGslMatrix->data[i*numCols+j];
    191             }
    192         }
    193     }
     186            memcpy(out->data.F64, in->data, numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
     187        }
     188        break;
     189      default:
     190        psAbort("Unsupported vector type: %x\n", out->type.type);
     191    }
     192    return;
    194193}
    195194
    … …  
    245244
    246245    // Copy psImage data into GSL matrix data
    247     psImageToGslMatrix(lu, in);
     246    matrixPStoGSL(lu, in);
    248247
    249248    // Calculate LU decomposition
    … …  
    251250
    252251    // Copy GSL matrix data to psImage data
    253     gslMatrixToPsImage(out, lu);
     252    matrixGSLtoPS(out, lu);
    254253
    255254    // Free GSL data
    … …  
    293292    // Initialize GSL data
    294293    lu = gsl_matrix_alloc(numRows, numCols);
    295     psImageToGslMatrix(lu, LU);
     294    matrixPStoGSL(lu, LU);
    296295    b = gsl_vector_alloc(RHS->n);
    297     psVectorToGslVector(b, RHS);
     296    vectorPStoGSL(b, RHS);
    298297    x = gsl_vector_alloc(RHS->n);
    299298
    … …  
    306305
    307306    // Copy GSL vector data to psVector data
    308     gslVectorToPsVector(out, x);
     307    vectorGSLtoPS(out, x);
    309308
    310309    // Free GSL data
    … …  
    341340    // Initialize GSL data
    342341    lu = gsl_matrix_alloc(numRows, numCols);
    343     psImageToGslMatrix(lu, LU);
     342    matrixPStoGSL(lu, LU);
    344343
    345344    permGSL.size = perm->n;
    … …  
    352351
    353352    // Copy GSL vector data to psVector data
    354     gslMatrixToPsImage(out, inverse);
     353    matrixGSLtoPS(out, inverse);
    355354
    356355    // Free GSL data
    … …  
    701700    lu = gsl_matrix_alloc(numRows, numCols);
    702701    inv = gsl_matrix_alloc(numRows, numCols);
    703     psImageToGslMatrix(lu, in);
     702    matrixPStoGSL(lu, in);
    704703
    705704    // Invert data and calculate determinant
    … …  
    716715
    717716    // Copy GSL matrix data to psImage data
    718     gslMatrixToPsImage(out, inv);
     717    matrixGSLtoPS(out, inv);
    719718
    720719    // Free GSL structs
    … …  
    749748    perm = gsl_permutation_alloc(numRows);
    750749    lu = gsl_matrix_alloc(numRows, numCols);
    751     psImageToGslMatrix(lu, in);
     750    matrixPStoGSL(lu, in);
    752751
    753752    // Calculate determinant
    … …  
    877876
    878877    inGSL = gsl_matrix_alloc(numRows, numCols);
    879     psImageToGslMatrix(inGSL, in);
     878    matrixPStoGSL(inGSL, in);
    880879    outGSL = gsl_matrix_alloc(numRows, numCols);
    881880
    … …  
    892891
    893892    // Copy GSL matrix data to psImage data
    894     gslMatrixToPsImage(out, outGSL);
     893    matrixGSLtoPS(out, outGSL);
    895894
    896895    // Free GSL structs
    … …  
    10261025}
    10271026
     1027psVector *psMatrixSolveSVD(const psImage *matrix, const psVector *vector)
     1028{
     1029    PS_ASSERT_IMAGE_NON_NULL(matrix, NULL);
     1030    PS_CHECK_DIMEN_AND_TYPE(matrix, PS_DIMEN_IMAGE, {return NULL;});
     1031    PS_ASSERT_VECTOR_NON_NULL(vector, NULL);
     1032    PS_CHECK_DIMEN_AND_TYPE(vector, PS_DIMEN_VECTOR, {return NULL;});
     1033
     1034    int numCols = matrix->numCols, numRows = matrix->numRows; // Size of matrix
     1035
     1036    // Decompose matrix: A = U S V^T
     1037    gsl_matrix *A = gsl_matrix_alloc(numRows, numCols); // Input matrix in GSL-speak; becomes matrix U
     1038    gsl_matrix *V = gsl_matrix_alloc(numCols, numCols); // Untransposed matrix V
     1039    gsl_vector *S = gsl_vector_alloc(numCols);          // Singular values
     1040    gsl_vector *work = gsl_vector_alloc(numCols);       // Work space for GSL
     1041
     1042    matrixPStoGSL(A, matrix);
     1043
     1044    int gslStatus = 0;                  // Status of GSL
     1045    if ((gslStatus = gsl_linalg_SV_decomp(A, V, S, work))) {
     1046        const char *err = gsl_strerror(gslStatus);
     1047        psError(PS_ERR_UNKNOWN, true, "Unable to decompose matrix: %s", err);
     1048        gsl_matrix_free(A);
     1049        gsl_matrix_free(V);
     1050        gsl_vector_free(S);
     1051        gsl_vector_free(work);
     1052        return NULL;
     1053    }
     1054    gsl_vector_free(work);
     1055
     1056    // Solve system (or minimise least-squares if overconstrained): Ax = b
     1057    gsl_vector *b = gsl_vector_alloc(numCols); // Vector b
     1058    gsl_vector *x = gsl_vector_alloc(numCols); // Solution
     1059
     1060    vectorPStoGSL(b, vector);
     1061
     1062    if ((gslStatus = gsl_linalg_SV_solve(A, V, S, b, x))) {
     1063        const char *err = gsl_strerror(gslStatus);
     1064        psError(PS_ERR_UNKNOWN, true, "Unable to solve matrix equation: %s", err);
     1065        gsl_matrix_free(A);
     1066        gsl_matrix_free(V);
     1067        gsl_vector_free(S);
     1068        gsl_vector_free(b);
     1069        gsl_vector_free(x);
     1070        return NULL;
     1071    }
     1072
     1073    gsl_matrix_free(A);
     1074    gsl_matrix_free(V);
     1075    gsl_vector_free(S);
     1076    gsl_vector_free(b);
     1077
     1078    psVector *solution = psVectorAlloc(numCols, PS_TYPE_F64);
     1079
     1080    vectorGSLtoPS(solution, x);
     1081    gsl_vector_free(x);
     1082
     1083    return solution;
     1084}
     1085
     1086
     1087
    10281088// This code supplied by Andy Becker (becker@astro.washington.edu)
    10291089psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in)
    … …  
    10501110
    10511111    // Copy psImage data into GSL matrix data
    1052     psImageToGslMatrix(A, in);
     1112    matrixPStoGSL(A, in);
    10531113
    10541114    // Calculate SVD decomposition
    … …  
    10561116
    10571117    // Copy GSL matrix data to psImage data
    1058     gslMatrixToPsImage(evec, V);
    1059     gslVectorToPsVector(eval, S);
     1118    matrixGSLtoPS(evec, V);
     1119    vectorGSLtoPS(eval, S);
    10601120
    10611121    // Take the square root of eval
Note: See TracChangeset for help on using the changeset viewer.