3 files changed, 71 insertions, 17 deletions

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