/** @file  tst_psMatrix_03.c
 *
 *  @brief Test driver for psMatrix LU functions
 *
 *  This test driver contains the following tests for psMatrix test point 3:
 *     A)  Create input and output images and vectors
 *     B)  Calculate LU matrix
 *     C)  Determine solution to matrix equation
 *     D)  Free input and output images and vectors
 *     E)  Attempt to use null image input argument
 *     F)  Attempt to use null input vector argument
 *     G)  ttempt to use null LU image argument
 *
 *  @author  Ross Harman, MHPCC
 *
 *  @version $Revision: 1.13 $  $Name: not supported by cvs2svn $
 *  @date  $Date: 2005-01-17 22:17:29 $
 *
 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
 *
 */

#include "pslib.h"
#include "psTest.h"

#define TOLERANCE 0.000001

#define CHECK_MATRIX(IMAGE)                                                                                  \
for(psU32 i=0; i<IMAGE->numRows; i++) {                                                                  \
    for(psU32 j=0; j<IMAGE->numCols; j++) {                                                              \
        if(IMAGE->type.type == PS_TYPE_F64) {                                                            \
            if(fabs(IMAGE->data.F64[i][j]-truthMatrix[i][j]) > TOLERANCE) {                              \
                printf("Matrix values at element %d, %d don't agree %lf vs %lf\n", i, j,                 \
                       IMAGE->data.F64[i][j], truthMatrix[i][j]);                                        \
            }                                                                                            \
        } else if(IMAGE->type.type == PS_TYPE_F32){                                                      \
            if(fabs(IMAGE->data.F32[i][j]-truthMatrix[i][j]) > TOLERANCE) {                              \
                printf("Matrix values at element %d, %d don't agree %f vs %lf\n", i, j,                  \
                       IMAGE->data.F32[i][j], truthMatrix[i][j]);                                        \
            }                                                                                            \
        }                                                                                                \
    }                                                                                                    \
}

#define CHECK_VECTOR(VECTOR)                                                                                 \
for(psU32 i=0; i<VECTOR->n; i++) {                                                                       \
    if(VECTOR->type.type == PS_TYPE_F64) {                                                               \
        if(fabs(VECTOR->data.F64[i]-truthVector[i]) > TOLERANCE) {                                       \
            printf("Vector values at element %d don't agree %lf vs %lf\n", i,                            \
                   VECTOR->data.F64[i], truthVector[i]);                                                 \
        }                                                                                                \
    } else if(VECTOR->type.type == PS_TYPE_F32){                                                         \
        if(fabs(VECTOR->data.F32[i]-truthVector[i]) > TOLERANCE) {                                       \
            printf("Vector values at element %d don't agree %f vs %lf\n", i,                             \
                   VECTOR->data.F32[i], truthVector[i]);                                                 \
        }                                                                                                \
    }                                                                                                    \
}


psS32 main(psS32 argc, char* argv[])
{
    psImage *luImage = NULL;
    psImage *inImage = NULL;
    psImage *tempImage = NULL;
    psVector *tempVector = NULL;
    psVector *perm = NULL;
    psVector *outVector = NULL;
    psVector *inVector = NULL;
    psImage *luImage32 = NULL;
    psImage *inImage32 = NULL;
    psImage *tempImage32 = NULL;
    psVector *tempVector32 = NULL;
    psVector *perm32 = NULL;
    psVector *outVector32 = NULL;
    psVector *inVector32 = NULL;

    double truthVector[3] = {
                                4.000000,
                                -2.000000,
                                3.000000
                            };

    double truthMatrix[3][3] = {{4.000000,  5.000000,  6.000000},
                                {0.750000, -2.750000, -6.500000},
                                {0.500000, -0.545455, -0.545455}};


    // Test A - Create input and output images and vectors
    printPositiveTestHeader(stdout, "psMatrix", "Create input and output images and vectors");
    luImage = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64);
    perm = (psVector*)psVectorAlloc(3, PS_TYPE_F64);
    outVector = (psVector*)psVectorAlloc(3, PS_TYPE_F64);
    inVector = (psVector*)psVectorAlloc(3, PS_TYPE_F64);
    inImage = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64);
    inImage->data.F64[0][0] =  2;
    inImage->data.F64[0][1] =  4;
    inImage->data.F64[0][2] =  6;
    inImage->data.F64[1][0] =  4;
    inImage->data.F64[1][1] =  5;
    inImage->data.F64[1][2] =  6;
    inImage->data.F64[2][0] =  3;
    inImage->data.F64[2][1] =  1;
    inImage->data.F64[2][2] = -2;
    inVector->data.F64[0] = 18.0;
    inVector->data.F64[1] = 24.0;
    inVector->data.F64[2] =  4.0;
    inVector->n = 3;
    luImage32 = (psImage*)psImageAlloc(3, 3, PS_TYPE_F32);
    outVector32 = (psVector*)psVectorAlloc(3, PS_TYPE_F32);
    inVector32 = (psVector*)psVectorAlloc(3, PS_TYPE_F32);
    inImage32 = (psImage*)psImageAlloc(3, 3, PS_TYPE_F32);
    inImage32->data.F32[0][0] =  2;
    inImage32->data.F32[0][1] =  4;
    inImage32->data.F32[0][2] =  6;
    inImage32->data.F32[1][0] =  4;
    inImage32->data.F32[1][1] =  5;
    inImage32->data.F32[1][2] =  6;
    inImage32->data.F32[2][0] =  3;
    inImage32->data.F32[2][1] =  1;
    inImage32->data.F32[2][2] = -2;
    inVector32->data.F32[0] = 18.0;
    inVector32->data.F32[1] = 24.0;
    inVector32->data.F32[2] =  4.0;
    inVector32->n = 3;
    printFooter(stdout, "psMatrix", "Create input and output images and vectors", true);


    // Test B - Calculate LU matrix
    printPositiveTestHeader(stdout, "psMatrix", "Calculate LU matrix");
    tempImage = luImage;
    luImage = psMatrixLUD(luImage, &perm, inImage);
    CHECK_MATRIX(luImage);
    if(luImage->type.dimen != PS_DIMEN_IMAGE) {
        printf("Error: Resulting image is not PS_DIMEN_IMAGE\n");
    } else if(luImage != tempImage) {
        printf("Error: Return pointer not equal to output argument pointer\n");
    }

    tempImage32 = luImage32;
    luImage32 = psMatrixLUD(luImage32, &perm32, inImage32);
    CHECK_MATRIX(luImage32);
    if(luImage32->type.dimen != PS_DIMEN_IMAGE) {
        printf("Error: Resulting image is not PS_DIMEN_IMAGE\n");
    } else if(luImage32 != tempImage32) {
        printf("Error: Return pointer not equal to output argument pointer\n");
    }
    printFooter(stdout, "psMatrix", "Calculate LU matrix", true);


    // Test C - Determine solution to matrix equation
    printPositiveTestHeader(stdout, "psMatrix", "Determine solution to matrix equation");
    tempVector = outVector;
    outVector = psMatrixLUSolve(outVector, luImage, inVector, perm);
    CHECK_VECTOR(outVector);
    if(outVector->type.dimen != PS_DIMEN_VECTOR) {
        printf("Error: Resulting image is not PS_DIMEN_VECTOR\n");
    } else if(outVector != tempVector) {
        printf("Error: Return pointer not equal to output argument pointer\n");
    }

    tempVector32 = outVector32;
    outVector32 = psMatrixLUSolve(outVector32, luImage32, inVector32, perm32);
    CHECK_VECTOR(outVector32);
    if(outVector32->type.dimen != PS_DIMEN_VECTOR) {
        printf("Error: Resulting image is not PS_DIMEN_VECTOR\n");
    } else if(outVector32 != tempVector32) {
        printf("Error: Return pointer not equal to output argument pointer\n");
    }
    printFooter(stdout, "psMatrix", "Determine solution to matrix equation", true);


    // Test D - Free input and output images and vectors
    printPositiveTestHeader(stdout, "psMatrix", "Free input and output images and vectors");
    psFree(inImage);
    psFree(luImage);
    psFree(perm);
    psFree(outVector);
    psFree(inVector);
    psFree(inImage32);
    psFree(luImage32);
    psFree(perm32);
    psFree(outVector32);
    psFree(inVector32);
    if( psMemCheckLeaks(0, NULL, stdout, false) ) {
        psError(PS_ERR_UNKNOWN,true,"Memory leaks detected");
        return 10;
    }
    psS32 nBad = psMemCheckCorruption(0);
    if(nBad) {
        printf("ERROR: Found %d bad memory blocks\n", nBad);
    }
    printFooter(stdout, "psMatrix" ,"Free input and output images and vectors", true);


    // Test E - Attempt to use null image input argument
    printNegativeTestHeader(stdout,"psMatrix", "Attempt to use null image input argument",
                            "Invalid operation: inImage or its data is NULL.", 0);
    psImage *imageTest = (psImage*)psImageAlloc(3, 3, PS_TYPE_F64);
    psMatrixLUD(imageTest, NULL, NULL);
    printFooter(stdout, "psMatrix", "Attempt to use null image input argument", true);


    // Test F - Attempt to use null input vector argument
    printNegativeTestHeader(stdout,"psMatrix", "Attempt to use null input vector argument",
                            "Invalid operation: inVector or its data is NULL.", 0);
    psVector *vectorBad = NULL;
    psVector *vectorBadOut = (psVector*)psVectorAlloc(3, PS_TYPE_F64);
    psVector *permBad = (psVector*)psVectorAlloc(3, PS_TYPE_F64);
    psMatrixLUSolve(vectorBadOut, imageTest, vectorBad, permBad);
    printFooter(stdout, "psMatrix", "Attempt to use null input vector argument", true);


    // Test G - Attempt to use null LU image argument
    printNegativeTestHeader(stdout,"psMatrix", "Attempt to use null LU image argument",
                            "Invalid operation: inImage or its data is NULL.", 0);
    psMatrixLUSolve(vectorBadOut, NULL, vectorBad, permBad);
    printFooter(stdout, "psMatrix", "Attempt to use null LU image argument", true);

    return 0;
}
