Index: trunk/psLib/test/math/tap_psMatrix08.c =================================================================== --- trunk/psLib/test/math/tap_psMatrix08.c (revision 24020) +++ trunk/psLib/test/math/tap_psMatrix08.c (revision 24087) @@ -1,16 +1,11 @@ /** @file tap_psMatrix_08.c * -* @brief Test driver for psMatrix transpose function -* -* This test driver contains the following tests: -* Transpose input image into output image -* Transpose input image into auto allocated NULL output image -* -* @author Ross Harman, MHPCC +* @brief psMatrixLUSolve, psMatrixGJSolve tests (for ill-conditioned matrix) +* @author Eugene Magnier, IfA * * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ * @date $Date: 2007-05-02 04:20:06 $ * -* Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii +* Copyright 2004-2005 Institute for Astronomy, University of Hawaii * */ @@ -21,33 +16,10 @@ #include "pstap.h" -# if (0) -bool check_matrix(psImage *img) -{ - bool errorFlag = false; - - for(psU32 i=0; inumRows; i++) { - for(psU32 j=0; jnumCols; j++) { - if(img->type.type == PS_TYPE_F64) { - if(fabs(img->data.F64[i][j]-truthMatrix[i][j]) > TOLERANCE) { - diag("Matrix values at element %d, %d don't agree %lf vs %lf\n", i, j, - img->data.F64[i][j], truthMatrix[i][j]); - errorFlag = true; - } - } else if(img->type.type == PS_TYPE_F32) { - if(fabs(img->data.F32[i][j]-truthMatrix[i][j]) > TOLERANCE) { - diag("Matrix values at element %d, %d don't agree %f vs %lf\n", i, j, - img->data.F32[i][j], truthMatrix[i][j]); - errorFlag = true; - } - } - } - } - return(errorFlag); -} -# endif +# define DEBUG 1 psS32 main( psS32 argc, char* argv[] ) { - plan_tests(14); + plan_tests(23); + // psTraceSetLevel("psLib.math.psMatrixGJSolve", 4); // Transpose input image into output image @@ -55,20 +27,185 @@ psMemId id = psMemGetId(); - // we have a specific image which gave us trouble elsewhere: - psImage *inImage = psImageAlloc(2, 2, PS_TYPE_F64); - inImage->data.F64[0][0] = 0.0; - inImage->data.F64[0][1] = -0.000100588316855; - inImage->data.F64[1][0] = +0.000100588316855; - inImage->data.F64[1][1] = 4.4276353417e-11; - - psVector *inVector = psVectorAlloc (2, PS_TYPE_F64); - inVector->data.F64[0] = 0.0; - inVector->data.F64[1] = 9.06388443347e-08; + psFits *fits = NULL; + + // we have a specific image and vector pair which gave us trouble elsewhere: + // XXX this is an ill-conditioned matrix. LU Decomposition does not inform us that it is ill-conditioned. + // the result solves the equation, but what are the errors on the values? + fits = psFitsOpen ("data/Agj.fits", "r"); + ok(fits, "opened test image Agj.fits"); + + psImage *Aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Aimage, "loaded test image Agj.fits"); + + psImage *aimage = psImageCopy (NULL, Aimage, Aimage->type.type); + ok(aimage, "copied test image Agj.fits"); + + psFitsClose (fits); + + fits = psFitsOpen ("data/Bgj.fits", "r"); + ok(fits, "opened test image Bgj.fits"); + + psImage *Bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Aimage, "loaded test image Bgj.fits"); + + psFitsClose (fits); + + psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type); + ok(Bvector, "allocated B vector"); + + for (int i = 0; i < Bimage->numRows; i++) { + double value = psImageGet (Bimage, 0, i); + psVectorSet (Bvector, i, value); + } bool status; - status = psMatrixGJSolveF32(inImage, inVector); - ok(status, "psMatrixGJSolve returns true exit status"); - psFree(inImage); - psFree(inVector); + status = psMatrixLUSolve(Aimage, Bvector); + ok(!status, "psMatrixLUSolve correctly returns false for ill-conditioned matrix"); + +# if (DEBUG) + fprintf (stderr, "LU Solution:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = psVectorGet (Bvector, i); + double valerr = psImageGet (Aimage, i, i); + fprintf (stderr, "%f +/- %f\n", value, valerr); + } + + // calculate Ax and compare with B: + fprintf (stderr, "result:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = 0; + for (int j = 0; j < Bvector->n; j++) { + double tmpV = psVectorGet (Bvector, j); + double tmpI = psImageGet (aimage, j, i); + value += tmpV*tmpI; + } + double actual = psImageGet (Bimage, 0, i); + fprintf (stderr, "%f vs %f (delta: %f)\n", value, actual, actual - value); + } +# endif + + psFree(Aimage); + psFree(Bimage); + psFree(Bvector); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // Transpose input image into output image + { + psMemId id = psMemGetId(); + + psFits *fits = NULL; + + // we have a specific ill-conditioned matrix in Agj.fits. psMatrixGJSolve detects this and reports a failure. + fits = psFitsOpen ("data/Agj.fits", "r"); + ok(fits, "opened test image Agj.fits"); + + psImage *Aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Aimage, "loaded test image Agj.fits"); + + psFitsClose (fits); + + fits = psFitsOpen ("data/Bgj.fits", "r"); + ok(fits, "opened test image Bgj.fits"); + + psImage *Bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(Bimage, "loaded test image Bgj.fits"); + + psFitsClose (fits); + + psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type); + ok(Bvector, "allocated B vector"); + + for (int i = 0; i < Bimage->numRows; i++) { + double value = psImageGet (Bimage, 0, i); + psVectorSet (Bvector, i, value); + } + + bool status; + status = psMatrixGJSolve(Aimage, Bvector); + ok(!status, "psMatrixGJSolve correctly returns false for ill-conditioned matrix"); + +# if (DEBUG) + fprintf (stderr, "GJ Solution:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = psVectorGet (Bvector, i); + fprintf (stderr, "%f\n", value); + } +# endif + + psFree(Aimage); + psFree(Bimage); + psFree(Bvector); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } + + // Transpose input image into output image + { + psMemId id = psMemGetId(); + + psFits *fits = NULL; + + // we have a specific ill-conditioned matrix in Agj.fits. psMatrixGJSolve detects this and reports a failure. + fits = psFitsOpen ("data/Agj.fits", "r"); + ok(fits, "opened test image Agj.fits"); + + psImage *aimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(aimage, "loaded test image Agj.fits"); + + psImage *Aimage = psImageCopy (NULL, aimage, PS_TYPE_F64); + ok(Aimage, "converted test image to F64"); + + psFitsClose (fits); + + fits = psFitsOpen ("data/Bgj.fits", "r"); + ok(fits, "opened test image Bgj.fits"); + + psImage *bimage = psFitsReadImage (fits, psRegionSet(0,0,0,0), 0); + ok(bimage, "loaded test image Bgj.fits"); + + psImage *Bimage = psImageCopy (NULL, bimage, PS_TYPE_F64); + ok(Bimage, "converted test image to F64"); + + psFitsClose (fits); + + psVector *Bvector = psVectorAlloc (Bimage->numRows, Bimage->type.type); + ok(Bvector, "allocated B vector"); + + for (int i = 0; i < Bimage->numRows; i++) { + double value = psImageGet (Bimage, 0, i); + psVectorSet (Bvector, i, value); + } + + bool status; + status = psMatrixGJSolve(Aimage, Bvector); + ok(!status, "psMatrixGJSolve correctly returns false for ill-conditioned matrix"); + +# if (DEBUG) + fprintf (stderr, "GJ Solution:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = psVectorGet (Bvector, i); + double valerr = psImageGet (Aimage, i, i); + fprintf (stderr, "%f +/- %f\n", value, valerr); + } + + // calculate Ax and compare with B: + fprintf (stderr, "result:\n"); + for (int i = 0; i < Bimage->numRows; i++) { + double value = 0; + for (int j = 0; j < Bvector->n; j++) { + double tmpV = psVectorGet (Bvector, j); + double tmpI = psImageGet (aimage, j, i); + value += tmpV*tmpI; + } + double actual = psImageGet (Bimage, 0, i); + fprintf (stderr, "%f vs %f (delta: %f)\n", value, actual, actual - value); + } +# endif + + psFree(Aimage); + psFree(Bimage); + psFree(aimage); + psFree(bimage); + psFree(Bvector); ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); }