IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 24, 2005, 3:15:01 PM (21 years ago)
Author:
drobbin
Message:

changed arguments as requested in apidelta-report-cycle6

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psMatrix.c

    r4385 r4388  
    2121 *  @author Robert DeSonia, MHPCC
    2222 *
    23  *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
    24  *  @date $Date: 2005-06-25 00:51:28 $
     23 *  @version $Revision: 1.34 $ $Name: not supported by cvs2svn $
     24 *  @date $Date: 2005-06-25 01:15:01 $
    2525 *
    2626 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    185185/*****************************************************************************/
    186186
    187 psImage* psMatrixLUD(psImage* outImage, psVector** outPerm, psImage* inImage)
     187psImage* psMatrixLUD(psImage* out, psVector** perm, psImage* in)
    188188{
    189189    psS32 signum = 0;
    … …  
    194194
    195195
    196     #define psMatrixLUD_EXIT {psFree(outImage); return NULL;}
    197 
    198     // Error checks
    199     PS_ASSERT_GENERAL_IMAGE_NON_NULL(inImage, psMatrixLUD_EXIT);
    200     PS_CHECK_POINTERS(inImage, outImage, psMatrixLUD_EXIT);
    201     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, psMatrixLUD_EXIT);
    202     PS_ASSERT_GENERAL_PTR_NON_NULL(outPerm, psMatrixLUD_EXIT);
    203 
    204     outImage = psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
    205 
    206     PS_CHECK_SQUARE(inImage, psMatrixLUD_EXIT); // gsl_linalg_LU_decomp would fail on non-square input.
    207     PS_CHECK_SQUARE(outImage, psMatrixLUD_EXIT);
     196    #define psMatrixLUD_EXIT {psFree(out); return NULL;}
     197
     198    // Error checks
     199    PS_ASSERT_GENERAL_IMAGE_NON_NULL(in, psMatrixLUD_EXIT);
     200    PS_CHECK_POINTERS(in, out, psMatrixLUD_EXIT);
     201    PS_CHECK_DIMEN_AND_TYPE(in, PS_DIMEN_IMAGE, psMatrixLUD_EXIT);
     202    PS_ASSERT_GENERAL_PTR_NON_NULL(perm, psMatrixLUD_EXIT);
     203
     204    out = psImageRecycle(out, in->numCols, in->numRows, in->type.type);
     205
     206    PS_CHECK_SQUARE(in, psMatrixLUD_EXIT); // gsl_linalg_LU_decomp would fail on non-square input.
     207    PS_CHECK_SQUARE(out, psMatrixLUD_EXIT);
    208208
    209209    // Initialize data
    210     numRows = inImage->numRows;
    211     numCols = inImage->numCols;
     210    numRows = in->numRows;
     211    numCols = in->numCols;
    212212
    213213    // Initialize GSL data
    214214    permGSL.size = numCols;
    215215    if (sizeof(size_t) == 4) {
    216         *outPerm = psVectorRecycle(*outPerm, numCols, PS_TYPE_S32);
     216        *perm = psVectorRecycle(*perm, numCols, PS_TYPE_S32);
    217217    } else if (sizeof(size_t) == 8) {
    218         *outPerm = psVectorRecycle(*outPerm, numCols, PS_TYPE_S64);
     218        *perm = psVectorRecycle(*perm, numCols, PS_TYPE_S64);
    219219    } else {
    220220        psError(PS_ERR_UNKNOWN, true,
    … …  
    224224    }
    225225
    226     (*outPerm)->n = numCols;
    227     permGSL.data = (psPtr)((*outPerm)->data.U8);
     226    (*perm)->n = numCols;
     227    permGSL.data = (psPtr)((*perm)->data.U8);
    228228    lu = gsl_matrix_alloc(numRows, numCols);
    229229
    230230    // Copy psImage data into GSL matrix data
    231     psImageToGslMatrix(lu, inImage);
     231    psImageToGslMatrix(lu, in);
    232232
    233233    // Calculate LU decomposition
    … …  
    235235
    236236    // Copy GSL matrix data to psImage data
    237     gslMatrixToPsImage(outImage, lu);
     237    gslMatrixToPsImage(out, lu);
    238238
    239239    // Free GSL data
    240240    gsl_matrix_free(lu);
    241241
    242     return outImage;
    243 }
    244 
    245 psVector* psMatrixLUSolve(psVector* outVector, const psImage* inImage, const psVector* inVector,
    246                           const psVector* inPerm)
     242    return out;
     243}
     244
     245psVector* psMatrixLUSolve(psVector* out, const psImage* LU, const psVector* RHS,
     246                          const psVector* perm)
    247247{
    248248    psS32 numRows = 0;
    249249    psS32 numCols = 0;
    250250    gsl_matrix *lu;
    251     gsl_permutation perm;
     251    gsl_permutation permGSL;
    252252    gsl_vector *b = NULL;
    253253    gsl_vector *x = NULL;
    254254
    255     #define LUSOLVE_CLEANUP {psFree(outVector); return NULL;}
    256 
    257     // Error checks
    258     PS_ASSERT_GENERAL_IMAGE_NON_NULL(inImage, LUSOLVE_CLEANUP);
    259     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, LUSOLVE_CLEANUP);
    260     PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(inImage, LUSOLVE_CLEANUP);
    261     PS_ASSERT_GENERAL_VECTOR_NON_NULL(inVector, LUSOLVE_CLEANUP);
    262     PS_CHECK_DIMEN_AND_TYPE(inVector, PS_DIMEN_VECTOR, LUSOLVE_CLEANUP);
    263     PS_ASSERT_GENERAL_VECTOR_NON_NULL(inPerm, LUSOLVE_CLEANUP);
    264 
    265     outVector = psVectorRecycle(outVector, inImage->numRows, inImage->type.type);
    266 
    267     PS_CHECK_POINTERS(outVector, inVector, LUSOLVE_CLEANUP);
    268     PS_CHECK_POINTERS(inVector, inPerm, LUSOLVE_CLEANUP);
    269     PS_CHECK_POINTERS(outVector, inPerm, LUSOLVE_CLEANUP);
     255    #define LUSOLVE_CLEANUP {psFree(out); return NULL;}
     256
     257    // Error checks
     258    PS_ASSERT_GENERAL_IMAGE_NON_NULL(LU, LUSOLVE_CLEANUP);
     259    PS_CHECK_DIMEN_AND_TYPE(LU, PS_DIMEN_IMAGE, LUSOLVE_CLEANUP);
     260    PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(LU, LUSOLVE_CLEANUP);
     261    PS_ASSERT_GENERAL_VECTOR_NON_NULL(RHS, LUSOLVE_CLEANUP);
     262    PS_CHECK_DIMEN_AND_TYPE(RHS, PS_DIMEN_VECTOR, LUSOLVE_CLEANUP);
     263    PS_ASSERT_GENERAL_VECTOR_NON_NULL(perm, LUSOLVE_CLEANUP);
     264
     265    out = psVectorRecycle(out, LU->numRows, LU->type.type);
     266
     267    PS_CHECK_POINTERS(out, RHS, LUSOLVE_CLEANUP);
     268    PS_CHECK_POINTERS(RHS, perm, LUSOLVE_CLEANUP);
     269    PS_CHECK_POINTERS(out, perm, LUSOLVE_CLEANUP);
    270270
    271271    // Initialize data
    272     numRows = inImage->numRows;
    273     numCols = inImage->numCols;
     272    numRows = LU->numRows;
     273    numCols = LU->numCols;
    274274
    275275    // Initialize GSL data
    276276    lu = gsl_matrix_alloc(numRows, numCols);
    277     psImageToGslMatrix(lu, inImage);
    278     b = gsl_vector_alloc(inVector->n);
    279     psVectorToGslVector(b, inVector);
    280     x = gsl_vector_alloc(inVector->n);
    281 
    282     outVector->n = numCols;
    283     perm.size = inPerm->n;
    284     perm.data = (psPtr)(inPerm->data.U8);
     277    psImageToGslMatrix(lu, LU);
     278    b = gsl_vector_alloc(RHS->n);
     279    psVectorToGslVector(b, RHS);
     280    x = gsl_vector_alloc(RHS->n);
     281
     282    out->n = numCols;
     283    permGSL.size = perm->n;
     284    permGSL.data = (psPtr)(perm->data.U8);
    285285
    286286    // Solve for {x} in equation: {b} = [A]{x}
    287     gsl_linalg_LU_solve(lu, &perm, b, x);
     287    gsl_linalg_LU_solve(lu, &permGSL, b, x);
    288288
    289289    // Copy GSL vector data to psVector data
    290     gslVectorToPsVector(outVector, x);
     290    gslVectorToPsVector(out, x);
    291291
    292292    // Free GSL data
    … …  
    295295    gsl_matrix_free(lu);
    296296
    297     return outVector;
     297    return out;
    298298}
    299299
    … …  
    383383}
    384384
    385 psImage* psMatrixMultiply(psImage* outImage, psImage* inImage1, psImage* inImage2)
     385psImage* psMatrixMultiply(psImage* out, psImage* in1, psImage* in2)
    386386{
    387387    gsl_matrix *m1 = NULL;
    … …  
    389389    gsl_matrix *m3 = NULL;
    390390
    391     #define MULTIPLY_CLEANUP { psFree(outImage); return NULL; }
    392 
    393     // Error checks
    394     PS_ASSERT_GENERAL_IMAGE_NON_NULL(inImage1, MULTIPLY_CLEANUP);
    395     PS_ASSERT_GENERAL_IMAGE_NON_NULL(inImage2, MULTIPLY_CLEANUP);
    396     PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(inImage1, MULTIPLY_CLEANUP);
    397     PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(inImage2, MULTIPLY_CLEANUP);
    398     PS_CHECK_DIMEN_AND_TYPE(inImage1, PS_DIMEN_IMAGE, MULTIPLY_CLEANUP);
    399     PS_CHECK_DIMEN_AND_TYPE(inImage2, PS_DIMEN_IMAGE, MULTIPLY_CLEANUP);
    400     PS_CHECK_POINTERS(inImage1, outImage, MULTIPLY_CLEANUP);
    401     PS_CHECK_POINTERS(inImage1, inImage2, MULTIPLY_CLEANUP);
    402 
    403     if (inImage1->numRows != inImage2->numCols) {
    404         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: number of rows of inImage1 != number of cols of inImage2.");
     391    #define MULTIPLY_CLEANUP { psFree(out); return NULL; }
     392
     393    // Error checks
     394    PS_ASSERT_GENERAL_IMAGE_NON_NULL(in1, MULTIPLY_CLEANUP);
     395    PS_ASSERT_GENERAL_IMAGE_NON_NULL(in2, MULTIPLY_CLEANUP);
     396    PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(in1, MULTIPLY_CLEANUP);
     397    PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(in2, MULTIPLY_CLEANUP);
     398    PS_CHECK_DIMEN_AND_TYPE(in1, PS_DIMEN_IMAGE, MULTIPLY_CLEANUP);
     399    PS_CHECK_DIMEN_AND_TYPE(in2, PS_DIMEN_IMAGE, MULTIPLY_CLEANUP);
     400    PS_CHECK_POINTERS(in1, out, MULTIPLY_CLEANUP);
     401    PS_CHECK_POINTERS(in1, in2, MULTIPLY_CLEANUP);
     402
     403    if (in1->numRows != in2->numCols) {
     404        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: number of rows of in1 != number of cols of in2.");
    405405        MULTIPLY_CLEANUP;
    406406    }
    407     if (inImage1->type.type != inImage2->type.type) {
    408         psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: data types of inImage1 and inImage2 must match.");
     407    if (in1->type.type != in2->type.type) {
     408        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Invalid operation: data types of in1 and in2 must match.");
    409409        MULTIPLY_CLEANUP;
    410410    }
    411411
    412     outImage = psImageRecycle(outImage, inImage2->numCols, inImage1->numRows, inImage2->type.type);
     412    out = psImageRecycle(out, in2->numCols, in1->numRows, in2->type.type);
    413413
    414414    // Initialize GSL data
    415     m1 = gsl_matrix_alloc(inImage1->numRows, inImage1->numCols);
    416     psImageToGslMatrix(m1, inImage1);
    417     m2 = gsl_matrix_alloc(inImage2->numRows, inImage2->numCols);
    418     psImageToGslMatrix(m2, inImage2);
    419     m3 = gsl_matrix_alloc(outImage->numRows, outImage->numCols);
     415    m1 = gsl_matrix_alloc(in1->numRows, in1->numCols);
     416    psImageToGslMatrix(m1, in1);
     417    m2 = gsl_matrix_alloc(in2->numRows, in2->numCols);
     418    psImageToGslMatrix(m2, in2);
     419    m3 = gsl_matrix_alloc(out->numRows, out->numCols);
    420420
    421421    // Perform multiplication
    … …  
    424424
    425425    // Copy GSL matrix data to psImage data
    426     gslMatrixToPsImage(outImage, m3);
     426    gslMatrixToPsImage(out, m3);
    427427
    428428    // Free GSL structs
    … …  
    431431    gsl_matrix_free(m3);
    432432
    433     return outImage;
    434 }
    435 
    436 psImage* psMatrixTranspose(psImage* outImage, const psImage* inImage)
     433    return out;
     434}
     435
     436psImage* psMatrixTranspose(psImage* out, const psImage* in)
    437437{
    438438    psU32 i = 0;
    … …  
    443443    psS32 numColsOut = 0;
    444444
    445     #define TRANSPOSE_CLEANUP { psFree(outImage); return NULL; }
    446     // Error checks
    447     PS_ASSERT_GENERAL_IMAGE_NON_NULL(inImage, TRANSPOSE_CLEANUP);
    448     PS_CHECK_DIMEN_AND_TYPE(inImage, PS_DIMEN_IMAGE, TRANSPOSE_CLEANUP);
    449     PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(inImage, TRANSPOSE_CLEANUP);
    450     PS_CHECK_POINTERS(inImage, outImage, TRANSPOSE_CLEANUP);
    451 
    452     outImage = psImageRecycle(outImage, inImage->numCols, inImage->numRows, inImage->type.type);
     445    #define TRANSPOSE_CLEANUP { psFree(out); return NULL; }
     446    // Error checks
     447    PS_ASSERT_GENERAL_IMAGE_NON_NULL(in, TRANSPOSE_CLEANUP);
     448    PS_CHECK_DIMEN_AND_TYPE(in, PS_DIMEN_IMAGE, TRANSPOSE_CLEANUP);
     449    PS_ASSERT_GENERAL_IMAGE_NON_EMPTY(in, TRANSPOSE_CLEANUP);
     450    PS_CHECK_POINTERS(in, out, TRANSPOSE_CLEANUP);
     451
     452    out = psImageRecycle(out, in->numCols, in->numRows, in->type.type);
    453453
    454454    // Initialize data
    455     numRowsIn = inImage->numRows;
    456     numColsIn = inImage->numCols;
    457     numRowsOut = outImage->numRows;
    458     numColsOut = outImage->numCols;
     455    numRowsIn = in->numRows;
     456    numColsIn = in->numCols;
     457    numRowsOut = out->numRows;
     458    numColsOut = out->numCols;
    459459
    460460    if(numRowsIn!=numColsOut && numRowsOut!=numColsIn) {
    … …  
    463463    }
    464464
    465     if(outImage->type.type == PS_TYPE_F32) {
     465    if(out->type.type == PS_TYPE_F32) {
    466466        for(i=0; i<numRowsOut; i++) {
    467467            for(j=0; j<numColsOut; j++) {
    468                 outImage->data.F32[i][j] = inImage->data.F32[j][i];
     468                out->data.F32[i][j] = in->data.F32[j][i];
    469469            }
    470470        }
    … …  
    472472        for(i=0; i<numRowsOut; i++) {
    473473            for(j=0; j<numColsOut; j++) {
    474                 outImage->data.F64[i][j] = inImage->data.F64[j][i];
     474                out->data.F64[i][j] = in->data.F64[j][i];
    475475            }
    476476        }
    477477    }
    478478
    479     return outImage;
     479    return out;
    480480}
    481481
Note: See TracChangeset for help on using the changeset viewer.