diff options
author | Alberto Ruiz <aruiz@um.es> | 2011-03-23 11:25:21 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2011-03-23 11:25:21 +0000 |
commit | a5a7d08b73223a9bd6022aebfcd32dfd42f6633a (patch) | |
tree | 7f1912f2324a91287a7c74bacb2793420b9bf322 /lib | |
parent | a325832d82be14b4ca7be2baa2a943f79629b194 (diff) |
findNaN
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK/lapack-aux.c | 62 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests.hs | 9 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests/Properties.hs | 17 |
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 @@ | |||
48 | void asm_finit() { | 48 | void asm_finit() { |
49 | #ifdef i386 | 49 | #ifdef i386 |
50 | 50 | ||
51 | #if FPUDEBUG | 51 | // asm("finit"); |
52 | uint val=0; | ||
53 | asm("fstsw" | ||
54 | : "=a" (val) | ||
55 | : "a" (val) | ||
56 | ); | ||
57 | 52 | ||
58 | val = (val%16384)/2048; // bits 13-11 | 53 | static unsigned char buf[108]; |
54 | asm("FSAVE %0":"=m" (buf)); | ||
59 | 55 | ||
60 | if (val != 0) { | 56 | #if FPUDEBUG |
61 | printf("Warning: FPU Stack: %d\n",val); | 57 | if(buf[8]!=255 || buf[9]!=255) { |
62 | asm("finit"); | 58 | printf("Warning: Nonempty FPU Stack. TAG = %x %x\n",buf[8],buf[9]); |
63 | } | 59 | } |
64 | #else | 60 | #endif |
65 | asm("finit"); | 61 | |
62 | #if NANDEBUG | ||
63 | asm("FRSTOR %0":"=m" (buf)); | ||
66 | #endif | 64 | #endif |
67 | 65 | ||
68 | #endif | 66 | #endif |
69 | } | 67 | } |
68 | |||
69 | //--------------------------------------- | ||
70 | |||
71 | #if NANDEBUG | ||
72 | |||
73 | #define CHECKNANR(M,msg) \ | ||
74 | { int k; \ | ||
75 | for(k=0; k<(M##r * M##c); k++) { \ | ||
76 | if(M##p[k] != M##p[k]) { \ | ||
77 | printf(msg); \ | ||
78 | TRACEMAT(M) \ | ||
79 | /*exit(1);*/ \ | ||
80 | } \ | ||
81 | } \ | ||
82 | } | ||
83 | |||
84 | #define CHECKNANC(M,msg) \ | ||
85 | { int k; \ | ||
86 | for(k=0; k<(M##r * M##c); k++) { \ | ||
87 | if( M##p[k].r != M##p[k].r \ | ||
88 | || M##p[k].i != M##p[k].i) { \ | ||
89 | printf(msg); \ | ||
90 | /*exit(1);*/ \ | ||
91 | } \ | ||
92 | } \ | ||
93 | } | ||
94 | |||
95 | #else | ||
96 | #define CHECKNANC(M,msg) | ||
97 | #define CHECKNANR(M,msg) | ||
98 | #endif | ||
99 | |||
70 | //--------------------------------------- | 100 | //--------------------------------------- |
71 | 101 | ||
72 | //////////////////// real svd //////////////////////////////////// | 102 | //////////////////// real svd //////////////////////////////////// |
@@ -1031,6 +1061,8 @@ void dgemm_(char *, char *, integer *, integer *, integer *, | |||
1031 | int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) { | 1061 | int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) { |
1032 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); | 1062 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); |
1033 | DEBUGMSG("dgemm_"); | 1063 | DEBUGMSG("dgemm_"); |
1064 | CHECKNANR(a,"NaN multR Input\n") | ||
1065 | CHECKNANR(b,"NaN multR Input\n") | ||
1034 | integer m = ta?ac:ar; | 1066 | integer m = ta?ac:ar; |
1035 | integer n = tb?br:bc; | 1067 | integer n = tb?br:bc; |
1036 | integer k = ta?ar:ac; | 1068 | integer k = ta?ar:ac; |
@@ -1040,6 +1072,7 @@ int multiplyR(int ta, int tb, KDMAT(a),KDMAT(b),DMAT(r)) { | |||
1040 | double alpha = 1; | 1072 | double alpha = 1; |
1041 | double beta = 0; | 1073 | double beta = 0; |
1042 | dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc); | 1074 | dgemm_(ta?"T":"N",tb?"T":"N",&m,&n,&k,&alpha,ap,&lda,bp,&ldb,&beta,rp,&ldc); |
1075 | CHECKNANR(r,"NaN multR Output\n") | ||
1043 | OK | 1076 | OK |
1044 | } | 1077 | } |
1045 | 1078 | ||
@@ -1050,6 +1083,8 @@ void zgemm_(char *, char *, integer *, integer *, integer *, | |||
1050 | int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) { | 1083 | int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) { |
1051 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); | 1084 | //REQUIRES(ac==br && ar==rr && bc==rc,BAD_SIZE); |
1052 | DEBUGMSG("zgemm_"); | 1085 | DEBUGMSG("zgemm_"); |
1086 | CHECKNANC(a,"NaN multC Input\n") | ||
1087 | CHECKNANC(b,"NaN multC Input\n") | ||
1053 | integer m = ta?ac:ar; | 1088 | integer m = ta?ac:ar; |
1054 | integer n = tb?br:bc; | 1089 | integer n = tb?br:bc; |
1055 | integer k = ta?ar:ac; | 1090 | integer k = ta?ar:ac; |
@@ -1062,6 +1097,7 @@ int multiplyC(int ta, int tb, KCMAT(a),KCMAT(b),CMAT(r)) { | |||
1062 | ap,&lda, | 1097 | ap,&lda, |
1063 | bp,&ldb,&beta, | 1098 | bp,&ldb,&beta, |
1064 | rp,&ldc); | 1099 | rp,&ldc); |
1100 | CHECKNANC(r,"NaN multC Output\n") | ||
1065 | OK | 1101 | OK |
1066 | } | 1102 | } |
1067 | 1103 | ||
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. | |||
17 | module Numeric.LinearAlgebra.Tests( | 17 | module Numeric.LinearAlgebra.Tests( |
18 | -- module Numeric.LinearAlgebra.Tests.Instances, | 18 | -- module Numeric.LinearAlgebra.Tests.Instances, |
19 | -- module Numeric.LinearAlgebra.Tests.Properties, | 19 | -- module Numeric.LinearAlgebra.Tests.Properties, |
20 | qCheck, runTests, runBenchmarks | 20 | qCheck, runTests, runBenchmarks, findNaN |
21 | --, runBigTests | 21 | --, runBigTests |
22 | ) where | 22 | ) where |
23 | 23 | ||
@@ -512,6 +512,8 @@ runTests n = do | |||
512 | putStrLn "------ chol" | 512 | putStrLn "------ chol" |
513 | test (cholProp . rPosDef) | 513 | test (cholProp . rPosDef) |
514 | test (cholProp . cPosDef) | 514 | test (cholProp . cPosDef) |
515 | test (exactProp . rPosDef) | ||
516 | test (exactProp . cPosDef) | ||
515 | putStrLn "------ expm" | 517 | putStrLn "------ expm" |
516 | test (expmDiagProp . complex. rSqWC) | 518 | test (expmDiagProp . complex. rSqWC) |
517 | test (expmDiagProp . cSqWC) | 519 | test (expmDiagProp . cSqWC) |
@@ -594,6 +596,11 @@ makeUnitary v | realPart n > 1 = v / scalar n | |||
594 | -- runBigTests :: IO () | 596 | -- runBigTests :: IO () |
595 | -- runBigTests = undefined | 597 | -- runBigTests = undefined |
596 | 598 | ||
599 | -- testcase for nonempty fpu stack | ||
600 | findNaN :: Int -> Bool | ||
601 | findNaN n = all (bugProp . eye) (take n $ cycle [1..20]) | ||
602 | where eye m = ident m :: Matrix ( Double) | ||
603 | |||
597 | -------------------------------------------------------------------------------- | 604 | -------------------------------------------------------------------------------- |
598 | 605 | ||
599 | -- | Performance measurements. | 606 | -- | 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 @@ | |||
1 | {-# LANGUAGE CPP #-} | 1 | {-# LANGUAGE CPP, FlexibleContexts #-} |
2 | {-# OPTIONS_GHC -fno-warn-unused-imports #-} | 2 | {-# OPTIONS_GHC -fno-warn-unused-imports #-} |
3 | ----------------------------------------------------------------------------- | 3 | ----------------------------------------------------------------------------- |
4 | {- | | 4 | {- | |
@@ -29,13 +29,14 @@ module Numeric.LinearAlgebra.Tests.Properties ( | |||
29 | pinvProp, | 29 | pinvProp, |
30 | detProp, | 30 | detProp, |
31 | nullspaceProp, | 31 | nullspaceProp, |
32 | bugProp, | ||
32 | svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, | 33 | svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, |
33 | svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, | 34 | svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, |
34 | eigProp, eigSHProp, eigProp2, eigSHProp2, | 35 | eigProp, eigSHProp, eigProp2, eigSHProp2, |
35 | qrProp, rqProp, rqProp1, rqProp2, rqProp3, | 36 | qrProp, rqProp, rqProp1, rqProp2, rqProp3, |
36 | hessProp, | 37 | hessProp, |
37 | schurProp1, schurProp2, | 38 | schurProp1, schurProp2, |
38 | cholProp, | 39 | cholProp, exactProp, |
39 | expmDiagProp, | 40 | expmDiagProp, |
40 | multProp1, multProp2, | 41 | multProp1, multProp2, |
41 | subProp, | 42 | subProp, |
@@ -132,6 +133,15 @@ nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) | |||
132 | 133 | ||
133 | ------------------------------------------------------------------ | 134 | ------------------------------------------------------------------ |
134 | 135 | ||
136 | -- testcase for nonempty fpu stack | ||
137 | -- uncommenting unitary' signature eliminates the problem | ||
138 | bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v | ||
139 | where (u,d,v) = fullSVD m | ||
140 | -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool | ||
141 | unitary' a = unitary a | ||
142 | |||
143 | ------------------------------------------------------------------ | ||
144 | |||
135 | -- fullSVD | 145 | -- fullSVD |
136 | svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v | 146 | svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v |
137 | where (u,d,v) = fullSVD m | 147 | where (u,d,v) = fullSVD m |
@@ -237,7 +247,8 @@ schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fix | |||
237 | 247 | ||
238 | cholProp m = m |~| ctrans c <> c && upperTriang c | 248 | cholProp m = m |~| ctrans c <> c && upperTriang c |
239 | where c = chol m | 249 | where c = chol m |
240 | -- pos = positiveDefinite m | 250 | |
251 | exactProp m = chol m == chol (m+0) | ||
241 | 252 | ||
242 | expmDiagProp m = expm (logm m) :~ 7 ~: complex m | 253 | expmDiagProp m = expm (logm m) :~ 7 ~: complex m |
243 | where logm = matFunc log | 254 | where logm = matFunc log |