From a5a7d08b73223a9bd6022aebfcd32dfd42f6633a Mon Sep 17 00:00:00 2001 From: Alberto Ruiz Date: Wed, 23 Mar 2011 11:25:21 +0000 Subject: findNaN --- lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 62 +++++++++++++++++++++------ lib/Numeric/LinearAlgebra/Tests.hs | 9 +++- lib/Numeric/LinearAlgebra/Tests/Properties.hs | 17 ++++++-- 3 files changed, 71 insertions(+), 17 deletions(-) (limited to 'lib') diff --git a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c index 77a8847..4bf29b7 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c +++ b/lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c @@ -48,25 +48,55 @@ void asm_finit() { #ifdef i386 - #if FPUDEBUG - uint val=0; - asm("fstsw" - : "=a" (val) - : "a" (val) - ); +// asm("finit"); - val = (val%16384)/2048; // bits 13-11 + static unsigned char buf[108]; + asm("FSAVE %0":"=m" (buf)); - if (val != 0) { - printf("Warning: FPU Stack: %d\n",val); - asm("finit"); - } - #else - asm("finit"); + #if FPUDEBUG + if(buf[8]!=255 || buf[9]!=255) { + printf("Warning: Nonempty FPU Stack. TAG = %x %x\n",buf[8],buf[9]); + } + #endif + + #if NANDEBUG + asm("FRSTOR %0":"=m" (buf)); #endif #endif } + +//--------------------------------------- + +#if NANDEBUG + +#define CHECKNANR(M,msg) \ +{ int k; \ +for(k=0; k<(M##r * M##c); k++) { \ + if(M##p[k] != M##p[k]) { \ + printf(msg); \ + TRACEMAT(M) \ + /*exit(1);*/ \ + } \ +} \ +} + +#define CHECKNANC(M,msg) \ +{ int k; \ +for(k=0; k<(M##r * M##c); k++) { \ + if( M##p[k].r != M##p[k].r \ + || M##p[k].i != M##p[k].i) { \ + printf(msg); \ + /*exit(1);*/ \ + } \ +} \ +} + +#else +#define CHECKNANC(M,msg) +#define CHECKNANR(M,msg) +#endif + //--------------------------------------- //////////////////// real svd //////////////////////////////////// @@ -1031,6 +1061,8 @@ void dgemm_(char *, char *, integer *, integer *, integer *, int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) { //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); DEBUGMSG("dgemm_"); + CHECKNANR(a,"NaN multR Input\n") + CHECKNANR(b,"NaN multR Input\n") integer m = ta?ac:ar; integer n = tb?br:bc; integer k = ta?ar:ac; @@ -1040,6 +1072,7 @@ int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) { double alpha = 1; double beta = 0; dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc); + CHECKNANR(r,"NaN multR Output\n") OK } @@ -1050,6 +1083,8 @@ void zgemm_(char *, char *, integer *, integer *, integer *, int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) { //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); DEBUGMSG("zgemm_"); + CHECKNANC(a,"NaN multC Input\n") + CHECKNANC(b,"NaN multC Input\n") integer m = ta?ac:ar; integer n = tb?br:bc; integer k = ta?ar:ac; @@ -1062,6 +1097,7 @@ int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) { ap,&lda, bp,&ldb,&beta, rp,&ldc); + CHECKNANC(r,"NaN multC Output\n") OK } diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index 5036e3f..e859450 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs @@ -17,7 +17,7 @@ Some tests. module Numeric.LinearAlgebra.Tests( -- module Numeric.LinearAlgebra.Tests.Instances, -- module Numeric.LinearAlgebra.Tests.Properties, - qCheck, runTests, runBenchmarks + qCheck, runTests, runBenchmarks, findNaN --, runBigTests ) where @@ -512,6 +512,8 @@ runTests n = do putStrLn "------ chol" test (cholProp . rPosDef) test (cholProp . cPosDef) + test (exactProp . rPosDef) + test (exactProp . cPosDef) putStrLn "------ expm" test (expmDiagProp . complex. rSqWC) test (expmDiagProp . cSqWC) @@ -594,6 +596,11 @@ makeUnitary v | realPart n > 1 = v / scalar n -- runBigTests :: IO () -- runBigTests = undefined +-- testcase for nonempty fpu stack +findNaN :: Int -> Bool +findNaN n = all (bugProp . eye) (take n $ cycle [1..20]) + where eye m = ident m :: Matrix ( Double) + -------------------------------------------------------------------------------- -- | Performance measurements. diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs index b5f0501..fe13544 100644 --- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs +++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE CPP #-} +{-# LANGUAGE CPP, FlexibleContexts #-} {-# OPTIONS_GHC -fno-warn-unused-imports #-} ----------------------------------------------------------------------------- {- | @@ -29,13 +29,14 @@ module Numeric.LinearAlgebra.Tests.Properties ( pinvProp, detProp, nullspaceProp, + bugProp, svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, eigProp, eigSHProp, eigProp2, eigSHProp2, qrProp, rqProp, rqProp1, rqProp2, rqProp3, hessProp, schurProp1, schurProp2, - cholProp, + cholProp, exactProp, expmDiagProp, multProp1, multProp2, subProp, @@ -132,6 +133,15 @@ nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) ------------------------------------------------------------------ +-- testcase for nonempty fpu stack +-- uncommenting unitary' signature eliminates the problem +bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v + where (u,d,v) = fullSVD m + -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool + unitary' a = unitary a + +------------------------------------------------------------------ + -- fullSVD svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v where (u,d,v) = fullSVD m @@ -237,7 +247,8 @@ schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fix cholProp m = m |~| ctrans c <> c && upperTriang c where c = chol m - -- pos = positiveDefinite m + +exactProp m = chol m == chol (m+0) expmDiagProp m = expm (logm m) :~ 7 ~: complex m where logm = matFunc log -- cgit v1.2.3