Index: trunk/psLib/test/math/Makefile.am =================================================================== --- trunk/psLib/test/math/Makefile.am (revision 24086) +++ trunk/psLib/test/math/Makefile.am (revision 24087) @@ -48,8 +48,9 @@ tap_psMatrixVectorArithmetic01 \ tap_psMatrixVectorArithmetic04 \ - tap_psRandom \ tap_psMinimizePowell \ tap_psSpline1D \ tap_psPolynomialMD + +# tap_psRandom if BUILD_TESTS Index: trunk/psLib/test/math/data/Agj.fits =================================================================== --- trunk/psLib/test/math/data/Agj.fits (revision 24087) +++ trunk/psLib/test/math/data/Agj.fits (revision 24087) @@ -0,0 +1,3 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = -32 / number of bits per data pixel NAXIS = 2 / number of data axes NAXIS1 = 9 / length of data axis 1 NAXIS2 = 9 / length of data axis 2 EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H BZERO = 0.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK BSCALE = 1.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK END =þÓŒE¢u>~Cºœ¥NOŒE¢v;€}&œ¥NP<Öñ?€>~C·œ¥NN@žÈÿ¿$ +kÀÙcÈ>FŸèœ¥NN<Öñ¿$ +k?õhŸ>Š_‡>FŸèÀN›>¿Œµ7>Š_‡?$*¿Œµ7ÀzþÀÙcÇ>FŸæBMI«>GK>FŸîÀN›A¿Œµ8>GOA奈@ÄŸ­¿Œµ8Àzÿ@ÄŸ®A€Ö Index: trunk/psLib/test/math/data/Bgj.fits =================================================================== --- trunk/psLib/test/math/data/Bgj.fits (revision 24087) +++ trunk/psLib/test/math/data/Bgj.fits (revision 24087) @@ -0,0 +1,2 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = -32 / number of bits per data pixel NAXIS = 2 / number of data axes NAXIS1 = 1 / length of data axis 1 NAXIS2 = 9 / length of data axis 2 EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H BZERO = 0.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK BSCALE = 1.000000000000E+00 / Scaling: TRUE = BZERO + BSCALE * DISK END ?ÃÒ+ŸþžšBÄTÖÀ +B+ÁŒWECtššClBÒ? Index: trunk/psLib/test/math/data/input =================================================================== --- trunk/psLib/test/math/data/input (revision 24087) +++ trunk/psLib/test/math/data/input (revision 24087) @@ -0,0 +1,9 @@ + +macro go + rd A Agj.fits + rd Bv Bgj.fits + mget Bv B -y 0 + set Ai = A + set Bi = B + gaussj A B +end Index: trunk/psLib/test/math/tap_psMatrix08.c =================================================================== --- trunk/psLib/test/math/tap_psMatrix08.c (revision 24086) +++ 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"); }