summaryrefslogtreecommitdiff
path: root/packages/tests
diff options
context:
space:
mode:
Diffstat (limited to 'packages/tests')
-rw-r--r--packages/tests/hmatrix-tests.cabal8
-rw-r--r--packages/tests/src/Numeric/GSL/Tests.hs62
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests.hs420
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs109
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs170
5 files changed, 471 insertions, 298 deletions
diff --git a/packages/tests/hmatrix-tests.cabal b/packages/tests/hmatrix-tests.cabal
index 0514843..49e0640 100644
--- a/packages/tests/hmatrix-tests.cabal
+++ b/packages/tests/hmatrix-tests.cabal
@@ -1,5 +1,5 @@
1Name: hmatrix-tests 1Name: hmatrix-tests
2Version: 0.4.1.0 2Version: 0.5.0.0
3License: BSD3 3License: BSD3
4License-file: LICENSE 4License-file: LICENSE
5Author: Alberto Ruiz 5Author: Alberto Ruiz
@@ -26,11 +26,11 @@ flag gsl
26 26
27library 27library
28 28
29 Build-Depends: base >= 4 && < 5, 29 Build-Depends: base >= 4 && < 5, deepseq,
30 QuickCheck >= 2, HUnit, random, 30 QuickCheck >= 2, HUnit, random,
31 hmatrix >= 0.16 31 hmatrix >= 0.17
32 if flag(gsl) 32 if flag(gsl)
33 Build-Depends: hmatrix-gsl >= 0.16 33 Build-Depends: hmatrix-gsl >= 0.17
34 34
35 hs-source-dirs: src 35 hs-source-dirs: src
36 36
diff --git a/packages/tests/src/Numeric/GSL/Tests.hs b/packages/tests/src/Numeric/GSL/Tests.hs
index 9dff6f5..025427b 100644
--- a/packages/tests/src/Numeric/GSL/Tests.hs
+++ b/packages/tests/src/Numeric/GSL/Tests.hs
@@ -19,10 +19,11 @@ import System.Exit (exitFailure)
19 19
20import Test.HUnit (runTestTT, failures, Test(..), errors) 20import Test.HUnit (runTestTT, failures, Test(..), errors)
21 21
22import Numeric.LinearAlgebra 22import Numeric.LinearAlgebra.HMatrix
23import Numeric.GSL 23import Numeric.GSL
24import Numeric.GSL.SimulatedAnnealing
24import Numeric.LinearAlgebra.Tests (qCheck, utest) 25import Numeric.LinearAlgebra.Tests (qCheck, utest)
25import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~)) 26import Numeric.LinearAlgebra.Tests.Properties ((|~|), (~~), (~=))
26 27
27--------------------------------------------------------------------- 28---------------------------------------------------------------------
28 29
@@ -42,7 +43,7 @@ fittingTest = utest "levmar" (ok1 && ok2)
42 sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0] 43 sol = fst $ fitModel 1E-4 1E-4 20 (expModel, expModelDer) dat [1,0,0]
43 44
44 ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d 45 ok1 = and (zipWith f sols [5,0.1,1]) where f (x,d) r = abs (x-r)<2*d
45 ok2 = norm2 (fromList (map fst sols) - fromList sol) < 1E-5 46 ok2 = norm_2 (fromList (map fst sols) - fromList sol) < 1E-5
46 47
47--------------------------------------------------------------------- 48---------------------------------------------------------------------
48 49
@@ -66,6 +67,59 @@ rootFindingTest = TestList [ utest "root Hybrids" (fst sol1 ~~ [1,1])
66 jacobian a b [x,_y] = [ [-a , 0] 67 jacobian a b [x,_y] = [ [-a , 0]
67 , [-2*b*x, b] ] 68 , [-2*b*x, b] ]
68 69
70--------------------------------------------------------------------
71
72interpolationTest = TestList [
73 utest "interpolation evaluateV" (esol ~= ev)
74 , utest "interpolation evaluate" (esol ~= eval)
75 , utest "interpolation evaluateDerivativeV" (desol ~= dev)
76 , utest "interpolation evaluateDerivative" (desol ~= de)
77 , utest "interpolation evaluateDerivative2V" (d2esol ~= d2ev)
78 , utest "interpolation evaluateDerivative2" (d2esol ~= d2e)
79 , utest "interpolation evaluateIntegralV" (intesol ~= intev)
80 , utest "interpolation evaluateIntegral" (intesol ~= inte)
81 ]
82 where
83 xtest = 2.2
84 applyVec f = f Akima xs ys xtest
85 applyList f = f Akima (zip xs' ys') xtest
86
87 esol = xtest**2
88 ev = applyVec evaluateV
89 eval = applyList evaluate
90
91 desol = 2*xtest
92 dev = applyVec evaluateDerivativeV
93 de = applyList evaluateDerivative
94
95 d2esol = 2
96 d2ev = applyVec evaluateDerivative2V
97 d2e = applyList evaluateDerivative2
98
99 intesol = 1/3 * xtest**3
100 intev = evaluateIntegralV Akima xs ys 0 xtest
101 inte = evaluateIntegral Akima (zip xs' ys') (0, xtest)
102
103 xs' = [-1..10]
104 ys' = map (**2) xs'
105 xs = vector xs'
106 ys = vector ys'
107
108---------------------------------------------------------------------
109
110simanTest = TestList [
111 -- We use a slightly more relaxed tolerance here because the
112 -- simulated annealer is randomized
113 utest "simulated annealing manual example" $ abs (result - 1.3631300) < 1e-6
114 ]
115 where
116 -- This is the example from the GSL manual.
117 result = simanSolve 0 1 exampleParams 15.5 exampleE exampleM exampleS Nothing
118 exampleParams = SimulatedAnnealingParams 200 10000 1.0 1.0 0.008 1.003 2.0e-6
119 exampleE x = exp (-(x - 1)**2) * sin (8 * x)
120 exampleM x y = abs $ x - y
121 exampleS rands stepSize current = (rands ! 0) * 2 * stepSize - stepSize + current
122
69--------------------------------------------------------------------- 123---------------------------------------------------------------------
70 124
71minimizationTest = TestList 125minimizationTest = TestList
@@ -123,6 +177,8 @@ runTests n = do
123 , odeTest 177 , odeTest
124 , rootFindingTest 178 , rootFindingTest
125 , minimizationTest 179 , minimizationTest
180 , interpolationTest
181 , simanTest
126 , utest "deriv" derivTest 182 , utest "deriv" derivTest
127 , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5**3) < 1E-8) 183 , utest "integrate" (abs (volSphere 2.5 - 4/3*pi*2.5**3) < 1E-8)
128 , utest "polySolve" (polySolveProp [1,2,3,4]) 184 , utest "polySolve" (polySolveProp [1,2,3,4])
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 713af79..4b631cf 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -4,6 +4,8 @@
4{-# LANGUAGE TypeFamilies #-} 4{-# LANGUAGE TypeFamilies #-}
5{-# LANGUAGE FlexibleContexts #-} 5{-# LANGUAGE FlexibleContexts #-}
6{-# LANGUAGE RankNTypes #-} 6{-# LANGUAGE RankNTypes #-}
7{-# LANGUAGE TypeOperators #-}
8{-# LANGUAGE ViewPatterns #-}
7 9
8----------------------------------------------------------------------------- 10-----------------------------------------------------------------------------
9{- | 11{- |
@@ -28,12 +30,9 @@ module Numeric.LinearAlgebra.Tests(
28--, runBigTests 30--, runBigTests
29) where 31) where
30 32
31import Numeric.LinearAlgebra 33import Numeric.LinearAlgebra hiding (unitary)
32import Numeric.LinearAlgebra.HMatrix hiding ((<>),linearSolve) 34import Numeric.LinearAlgebra.Devel
33import Numeric.LinearAlgebra.Static(L) 35import Numeric.LinearAlgebra.Static(L)
34import Numeric.LinearAlgebra.Util(col,row)
35import Data.Packed
36import Numeric.LinearAlgebra.LAPACK
37import Numeric.LinearAlgebra.Tests.Instances 36import Numeric.LinearAlgebra.Tests.Instances
38import Numeric.LinearAlgebra.Tests.Properties 37import Numeric.LinearAlgebra.Tests.Properties
39import Test.HUnit hiding ((~:),test,Testable,State) 38import Test.HUnit hiding ((~:),test,Testable,State)
@@ -44,15 +43,13 @@ import qualified Prelude
44import System.CPUTime 43import System.CPUTime
45import System.Exit 44import System.Exit
46import Text.Printf 45import Text.Printf
47import Data.Packed.Development(unsafeFromForeignPtr,unsafeToForeignPtr) 46import Numeric.LinearAlgebra.Devel(unsafeFromForeignPtr,unsafeToForeignPtr)
48import Control.Arrow((***)) 47import Control.Arrow((***))
49import Debug.Trace 48import Debug.Trace
50import Control.Monad(when) 49import Control.Monad(when)
51import Numeric.LinearAlgebra.Util hiding (ones,row,col)
52import Control.Applicative 50import Control.Applicative
53import Control.Monad(ap) 51import Control.Monad(ap)
54 52import Control.DeepSeq ( NFData(..) )
55import Data.Packed.ST
56 53
57import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 54import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector
58 ,sized,classify,Testable,Property 55 ,sized,classify,Testable,Property
@@ -81,7 +78,7 @@ detTest1 = det m == 26
81 && det mc == 38 :+ (-3) 78 && det mc == 38 :+ (-3)
82 && det (feye 2) == -1 79 && det (feye 2) == -1
83 where 80 where
84 m = (3><3) 81 m = (3><3)
85 [ 1, 2, 3 82 [ 1, 2, 3
86 , 4, 5, 7 83 , 4, 5, 7
87 , 2, 8, 4 :: Double 84 , 2, 8, 4 :: Double
@@ -89,7 +86,7 @@ detTest1 = det m == 26
89 mc = (3><3) 86 mc = (3><3)
90 [ 1, 2, 3 87 [ 1, 2, 3
91 , 4, 5, 7 88 , 4, 5, 7
92 , 2, 8, i 89 , 2, 8, iC
93 ] 90 ]
94 91
95detTest2 = inv1 |~| inv2 && [det1] ~~ [det2] 92detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]
@@ -130,8 +127,8 @@ expmTest2 = expm nd2 :~15~: (2><2)
130mbCholTest = utest "mbCholTest" (ok1 && ok2) where 127mbCholTest = utest "mbCholTest" (ok1 && ok2) where
131 m1 = (2><2) [2,5,5,8 :: Double] 128 m1 = (2><2) [2,5,5,8 :: Double]
132 m2 = (2><2) [3,5,5,9 :: Complex Double] 129 m2 = (2><2) [3,5,5,9 :: Complex Double]
133 ok1 = mbCholSH m1 == Nothing 130 ok1 = mbChol (trustSym m1) == Nothing
134 ok2 = mbCholSH m2 == Just (chol m2) 131 ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2)
135 132
136--------------------------------------------------------------------- 133---------------------------------------------------------------------
137 134
@@ -140,7 +137,7 @@ randomTestGaussian = c :~1~: snd (meanCov dat) where
140 2,4,0, 137 2,4,0,
141 -2,2,1] 138 -2,2,1]
142 m = 3 |> [1,2,3] 139 m = 3 |> [1,2,3]
143 c = a <> trans a 140 c = a <> tr a
144 dat = gaussianSample 7 (10^6) m c 141 dat = gaussianSample 7 (10^6) m c
145 142
146randomTestUniform = c :~1~: snd (meanCov dat) where 143randomTestUniform = c :~1~: snd (meanCov dat) where
@@ -174,54 +171,54 @@ offsetTest = y == y' where
174 171
175normsVTest = TestList [ 172normsVTest = TestList [
176 utest "normv2CD" $ norm2PropC v 173 utest "normv2CD" $ norm2PropC v
177 , utest "normv2CF" $ norm2PropC (single v) 174-- , utest "normv2CF" $ norm2PropC (single v)
178#ifndef NONORMVTEST 175#ifndef NONORMVTEST
179 , utest "normv2D" $ norm2PropR x 176 , utest "normv2D" $ norm2PropR x
180 , utest "normv2F" $ norm2PropR (single x) 177-- , utest "normv2F" $ norm2PropR (single x)
181#endif 178#endif
182 , utest "normv1CD" $ norm1 v == 8 179 , utest "normv1CD" $ norm_1 v == 8
183 , utest "normv1CF" $ norm1 (single v) == 8 180-- , utest "normv1CF" $ norm_1 (single v) == 8
184 , utest "normv1D" $ norm1 x == 6 181 , utest "normv1D" $ norm_1 x == 6
185 , utest "normv1F" $ norm1 (single x) == 6 182-- , utest "normv1F" $ norm_1 (single x) == 6
186 183
187 , utest "normvInfCD" $ normInf v == 5 184 , utest "normvInfCD" $ norm_Inf v == 5
188 , utest "normvInfCF" $ normInf (single v) == 5 185-- , utest "normvInfCF" $ norm_Inf (single v) == 5
189 , utest "normvInfD" $ normInf x == 3 186 , utest "normvInfD" $ norm_Inf x == 3
190 , utest "normvInfF" $ normInf (single x) == 3 187-- , utest "normvInfF" $ norm_Inf (single x) == 3
191 188
192 ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double) 189 ] where v = fromList [1,-2,3:+4] :: Vector (Complex Double)
193 x = fromList [1,2,-3] :: Vector Double 190 x = fromList [1,2,-3] :: Vector Double
194#ifndef NONORMVTEST 191#ifndef NONORMVTEST
195 norm2PropR a = norm2 a =~= sqrt (udot a a) 192 norm2PropR a = norm_2 a =~= sqrt (udot a a)
196#endif 193#endif
197 norm2PropC a = norm2 a =~= realPart (sqrt (a <.> a)) 194 norm2PropC a = norm_2 a =~= realPart (sqrt (a `dot` a))
198 a =~= b = fromList [a] |~| fromList [b] 195 a =~= b = fromList [a] |~| fromList [b]
199 196
200normsMTest = TestList [ 197normsMTest = TestList [
201 utest "norm2mCD" $ pnorm PNorm2 v =~= 8.86164970498005 198 utest "norm2mCD" $ norm_2 v =~= 8.86164970498005
202 , utest "norm2mCF" $ pnorm PNorm2 (single v) =~= 8.86164970498005 199-- , utest "norm2mCF" $ norm_2 (single v) =~= 8.86164970498005
203 , utest "norm2mD" $ pnorm PNorm2 x =~= 5.96667765076216 200 , utest "norm2mD" $ norm_2 x =~= 5.96667765076216
204 , utest "norm2mF" $ pnorm PNorm2 (single x) =~= 5.96667765076216 201-- , utest "norm2mF" $ norm_2 (single x) =~= 5.96667765076216
205 202
206 , utest "norm1mCD" $ pnorm PNorm1 v == 9 203 , utest "norm1mCD" $ norm_1 v == 9
207 , utest "norm1mCF" $ pnorm PNorm1 (single v) == 9 204-- , utest "norm1mCF" $ norm_1 (single v) == 9
208 , utest "norm1mD" $ pnorm PNorm1 x == 7 205 , utest "norm1mD" $ norm_1 x == 7
209 , utest "norm1mF" $ pnorm PNorm1 (single x) == 7 206-- , utest "norm1mF" $ norm_1 (single x) == 7
210 207
211 , utest "normmInfCD" $ pnorm Infinity v == 12 208 , utest "normmInfCD" $ norm_Inf v == 12
212 , utest "normmInfCF" $ pnorm Infinity (single v) == 12 209-- , utest "normmInfCF" $ norm_Inf (single v) == 12
213 , utest "normmInfD" $ pnorm Infinity x == 8 210 , utest "normmInfD" $ norm_Inf x == 8
214 , utest "normmInfF" $ pnorm Infinity (single x) == 8 211-- , utest "normmInfF" $ norm_Inf (single x) == 8
215 212
216 , utest "normmFroCD" $ pnorm Frobenius v =~= 8.88819441731559 213 , utest "normmFroCD" $ norm_Frob v =~= 8.88819441731559
217 , utest "normmFroCF" $ pnorm Frobenius (single v) =~~= 8.88819441731559 214-- , utest "normmFroCF" $ norm_Frob (single v) =~~= 8.88819441731559
218 , utest "normmFroD" $ pnorm Frobenius x =~= 6.24499799839840 215 , utest "normmFroD" $ norm_Frob x =~= 6.24499799839840
219 , utest "normmFroF" $ pnorm Frobenius (single x) =~~= 6.24499799839840 216-- , utest "normmFroF" $ norm_Frob (single x) =~~= 6.24499799839840
220 217
221 ] where v = (2><2) [1,-2*i,3:+4,7] :: Matrix (Complex Double) 218 ] where v = (2><2) [1,-2*iC,3:+4,7] :: Matrix (Complex Double)
222 x = (2><2) [1,2,-3,5] :: Matrix Double 219 x = (2><2) [1,2,-3,5] :: Matrix Double
223 a =~= b = fromList [a] :~10~: fromList [b] 220 a =~= b = fromList [a] :~10~: fromList [b]
224 a =~~= b = fromList [a] :~5~: fromList [b] 221-- a =~~= b = fromList [a] :~5~: fromList [b]
225 222
226--------------------------------------------------------------------- 223---------------------------------------------------------------------
227 224
@@ -236,7 +233,7 @@ sumprodTest = TestList [
236 , utest "prodD" $ prodProp v 233 , utest "prodD" $ prodProp v
237 , utest "prodF" $ prodProp (single v) 234 , utest "prodF" $ prodProp (single v)
238 ] where v = fromList [1,2,3] :: Vector Double 235 ] where v = fromList [1,2,3] :: Vector Double
239 z = fromList [1,2-i,3+i] 236 z = fromList [1,2-iC,3+iC]
240 prodProp x = prodElements x == product (toList x) 237 prodProp x = prodElements x == product (toList x)
241 238
242--------------------------------------------------------------------- 239---------------------------------------------------------------------
@@ -250,7 +247,7 @@ chainTest = utest "chain" $ foldl1' (<>) ms |~| optimiseMult ms where
250 247
251--------------------------------------------------------------------- 248---------------------------------------------------------------------
252 249
253conjuTest m = mapVector conjugate (flatten (trans m)) == flatten (ctrans m) 250conjuTest m = cmap conjugate (flatten (conj (tr m))) == flatten (tr m)
254 251
255--------------------------------------------------------------------- 252---------------------------------------------------------------------
256 253
@@ -306,7 +303,7 @@ lift_maybe m = MaybeT $ do
306 303
307-- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs 304-- apply a test to successive elements of a vector, evaluates to true iff test passes for all pairs
308--successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool 305--successive_ :: Storable a => (a -> a -> Bool) -> Vector a -> Bool
309successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (dim v - 1) v))) (v @> 0) 306successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ stp (subVector 1 (size v - 1) v))) (v ! 0)
310 where stp e = do 307 where stp e = do
311 ep <- lift_maybe $ state_get 308 ep <- lift_maybe $ state_get
312 if t e ep 309 if t e ep
@@ -315,7 +312,7 @@ successive_ t v = maybe False (\_ -> True) $ evalState (runMaybeT (mapVectorM_ s
315 312
316-- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input 313-- operate on successive elements of a vector and return the resulting vector, whose length 1 less than that of the input
317--successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b 314--successive :: (Storable a, Storable b) => (a -> a -> b) -> Vector a -> Vector b
318successive f v = evalState (mapVectorM stp (subVector 1 (dim v - 1) v)) (v @> 0) 315successive f v = evalState (mapVectorM stp (subVector 1 (size v - 1) v)) (v ! 0)
319 where stp e = do 316 where stp e = do
320 ep <- state_get 317 ep <- state_get
321 state_put e 318 state_put e
@@ -362,7 +359,7 @@ accumTest = utest "accum" ok
362 ,0,1,7 359 ,0,1,7
363 ,0,0,4] 360 ,0,0,4]
364 && 361 &&
365 toList (flatten x) == [1,0,0,0,1,0,0,0,1] 362 toList (flatten x) == [1,0,0,0,1,0,0,0,1]
366 363
367-------------------------------------------------------------------------------- 364--------------------------------------------------------------------------------
368 365
@@ -377,28 +374,19 @@ convolutionTest = utest "convolution" ok
377 374
378-------------------------------------------------------------------------------- 375--------------------------------------------------------------------------------
379 376
380kroneckerTest = utest "kronecker" ok 377sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix))
381 where
382 a,x,b :: Matrix Double
383 a = (3><4) [1..]
384 x = (4><2) [3,5..]
385 b = (2><5) [0,5..]
386 v1 = vec (a <> x <> b)
387 v2 = (trans b `kronecker` a) <> vec x
388 s = trans b <> b
389 v3 = vec s
390 v4 = (dup 5 :: Matrix Double) <> vech s
391 ok = v1 == v2 && v3 == v4
392 && vtrans 1 a == trans a
393 && vtrans (rows a) a == asColumn (vec a)
394 378
395-------------------------------------------------------------------------------- 379--------------------------------------------------------------------------------
396 380
397sparseTest = utest "sparse" (fst $ checkT (undefined :: GMatrix)) 381staticTest = utest "static" (fst $ checkT (undefined :: L 3 5))
398 382
399-------------------------------------------------------------------------------- 383--------------------------------------------------------------------------------
400 384
401staticTest = utest "static" (fst $ checkT (undefined :: L 3 5)) 385intTest = utest "int ops" (fst $ checkT (undefined :: Matrix I))
386
387--------------------------------------------------------------------------------
388
389modularTest = utest "modular ops" (fst $ checkT (undefined :: Matrix (Mod 13 I)))
402 390
403-------------------------------------------------------------------------------- 391--------------------------------------------------------------------------------
404 392
@@ -414,6 +402,150 @@ indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3
414 402
415-------------------------------------------------------------------------------- 403--------------------------------------------------------------------------------
416 404
405sliceTest = utest "slice test" $ and
406 [ testSlice (chol . trustSym) (gen 5 :: Matrix R)
407 , testSlice (chol . trustSym) (gen 5 :: Matrix C)
408 , testSlice qr (rec :: Matrix R)
409 , testSlice qr (rec :: Matrix C)
410 , testSlice hess (agen 5 :: Matrix R)
411 , testSlice hess (agen 5 :: Matrix C)
412 , testSlice schur (agen 5 :: Matrix R)
413 , testSlice schur (agen 5 :: Matrix C)
414 , testSlice lu (agen 5 :: Matrix R)
415 , testSlice lu (agen 5 :: Matrix C)
416 , testSlice (luSolve (luPacked (agen 5 :: Matrix R))) (agen 5)
417 , testSlice (luSolve (luPacked (agen 5 :: Matrix C))) (agen 5)
418 , test_lus (agen 5 :: Matrix R)
419 , test_lus (agen 5 :: Matrix C)
420
421 , testSlice eig (agen 5 :: Matrix R)
422 , testSlice eig (agen 5 :: Matrix C)
423 , testSlice (eigSH . trustSym) (gen 5 :: Matrix R)
424 , testSlice (eigSH . trustSym) (gen 5 :: Matrix C)
425 , testSlice eigenvalues (agen 5 :: Matrix R)
426 , testSlice eigenvalues (agen 5 :: Matrix C)
427 , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix R)
428 , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix C)
429
430 , testSlice svd (rec :: Matrix R)
431 , testSlice thinSVD (rec :: Matrix R)
432 , testSlice compactSVD (rec :: Matrix R)
433 , testSlice leftSV (rec :: Matrix R)
434 , testSlice rightSV (rec :: Matrix R)
435 , testSlice singularValues (rec :: Matrix R)
436
437 , testSlice svd (rec :: Matrix C)
438 , testSlice thinSVD (rec :: Matrix C)
439 , testSlice compactSVD (rec :: Matrix C)
440 , testSlice leftSV (rec :: Matrix C)
441 , testSlice rightSV (rec :: Matrix C)
442 , testSlice singularValues (rec :: Matrix C)
443
444 , testSlice (linearSolve (agen 5:: Matrix R)) (agen 5)
445 , testSlice (flip linearSolve (agen 5:: Matrix R)) (agen 5)
446
447 , testSlice (linearSolve (agen 5:: Matrix C)) (agen 5)
448 , testSlice (flip linearSolve (agen 5:: Matrix C)) (agen 5)
449
450 , testSlice (linearSolveLS (ogen 5:: Matrix R)) (ogen 5)
451 , testSlice (flip linearSolveLS (ogen 5:: Matrix R)) (ogen 5)
452
453 , testSlice (linearSolveLS (ogen 5:: Matrix C)) (ogen 5)
454 , testSlice (flip linearSolveLS (ogen 5:: Matrix C)) (ogen 5)
455
456 , testSlice (linearSolveSVD (ogen 5:: Matrix R)) (ogen 5)
457 , testSlice (flip linearSolveSVD (ogen 5:: Matrix R)) (ogen 5)
458
459 , testSlice (linearSolveSVD (ogen 5:: Matrix C)) (ogen 5)
460 , testSlice (flip linearSolveSVD (ogen 5:: Matrix C)) (ogen 5)
461
462 , testSlice (linearSolveLS (ugen 5:: Matrix R)) (ugen 5)
463 , testSlice (flip linearSolveLS (ugen 5:: Matrix R)) (ugen 5)
464
465 , testSlice (linearSolveLS (ugen 5:: Matrix C)) (ugen 5)
466 , testSlice (flip linearSolveLS (ugen 5:: Matrix C)) (ugen 5)
467
468 , testSlice (linearSolveSVD (ugen 5:: Matrix R)) (ugen 5)
469 , testSlice (flip linearSolveSVD (ugen 5:: Matrix R)) (ugen 5)
470
471 , testSlice (linearSolveSVD (ugen 5:: Matrix C)) (ugen 5)
472 , testSlice (flip linearSolveSVD (ugen 5:: Matrix C)) (ugen 5)
473
474 , testSlice ((<>) (ogen 5:: Matrix R)) (gen 5)
475 , testSlice (flip (<>) (gen 5:: Matrix R)) (ogen 5)
476 , testSlice ((<>) (ogen 5:: Matrix C)) (gen 5)
477 , testSlice (flip (<>) (gen 5:: Matrix C)) (ogen 5)
478 , testSlice ((<>) (ogen 5:: Matrix Float)) (gen 5)
479 , testSlice (flip (<>) (gen 5:: Matrix Float)) (ogen 5)
480 , testSlice ((<>) (ogen 5:: Matrix (Complex Float))) (gen 5)
481 , testSlice (flip (<>) (gen 5:: Matrix (Complex Float))) (ogen 5)
482 , testSlice ((<>) (ogen 5:: Matrix I)) (gen 5)
483 , testSlice (flip (<>) (gen 5:: Matrix I)) (ogen 5)
484 , testSlice ((<>) (ogen 5:: Matrix Z)) (gen 5)
485 , testSlice (flip (<>) (gen 5:: Matrix Z)) (ogen 5)
486
487 , testSlice ((<>) (ogen 5:: Matrix (I ./. 7))) (gen 5)
488 , testSlice (flip (<>) (gen 5:: Matrix (I ./. 7))) (ogen 5)
489 , testSlice ((<>) (ogen 5:: Matrix (Z ./. 7))) (gen 5)
490 , testSlice (flip (<>) (gen 5:: Matrix (Z ./. 7))) (ogen 5)
491
492 , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ trustSym $ gen 5)
493 , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ trustSym $ gen 5)
494 , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix R)) (agen 5)
495 , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix C)) (agen 5)
496
497 , ok_qrgr (rec :: Matrix R)
498 , ok_qrgr (rec :: Matrix C)
499 , testSlice (test_qrgr 4 tau1) qrr1
500 , testSlice (test_qrgr 4 tau2) qrr2
501 ]
502 where
503 QR qrr1 tau1 = qrRaw (rec :: Matrix R)
504 QR qrr2 tau2 = qrRaw (rec :: Matrix C)
505
506 test_qrgr n t x = qrgr n (QR x t)
507
508 ok_qrgr x = simeq 1E-15 q q'
509 where
510 (q,_) = qr x
511 atau = qrRaw x
512 q' = qrgr (rows q) atau
513
514 simeq eps a b = not $ magnit eps (norm_1 $ flatten (a-b))
515
516 test_lus m = testSlice f lup
517 where
518 f x = luSolve (LU x p) m
519 (LU lup p) = luPacked m
520
521 gen :: Numeric t => Int -> Matrix t
522 gen n = diagRect 1 (konst 5 n) n n
523
524 agen :: (Numeric t, Num (Vector t))=> Int -> Matrix t
525 agen n = gen n + fromInt ((n><n)[0..])
526
527 ogen :: (Numeric t, Num (Vector t))=> Int -> Matrix t
528 ogen n = gen n === gen n
529
530 ugen :: (Numeric t, Num (Vector t))=> Int -> Matrix t
531 ugen n = takeRows 3 (gen n)
532
533
534 rec :: Numeric t => Matrix t
535 rec = subMatrix (0,0) (4,5) (gen 5)
536
537 testSlice f x@(size->sz@(r,c)) = all (==f x) (map f (g y1 ++ g y2))
538 where
539 subm = subMatrix
540 g y = [ subm (a*r,b*c) sz y | a <-[0..2], b <- [0..2]]
541 h z = fromBlocks (replicate 3 (replicate 3 z))
542 y1 = h x
543 y2 = (tr . h . tr) x
544
545
546
547--------------------------------------------------------------------------------
548
417-- | All tests must pass with a maximum dimension of about 20 549-- | All tests must pass with a maximum dimension of about 20
418-- (some tests may fail with bigger sizes due to precision loss). 550-- (some tests may fail with bigger sizes due to precision loss).
419runTests :: Int -- ^ maximum dimension 551runTests :: Int -- ^ maximum dimension
@@ -435,11 +567,11 @@ runTests n = do
435 test (multProp1 10 . cConsist) 567 test (multProp1 10 . cConsist)
436 test (multProp2 10 . rConsist) 568 test (multProp2 10 . rConsist)
437 test (multProp2 10 . cConsist) 569 test (multProp2 10 . cConsist)
438 putStrLn "------ mult Float" 570-- putStrLn "------ mult Float"
439 test (multProp1 6 . (single *** single) . rConsist) 571-- test (multProp1 6 . (single *** single) . rConsist)
440 test (multProp1 6 . (single *** single) . cConsist) 572-- test (multProp1 6 . (single *** single) . cConsist)
441 test (multProp2 6 . (single *** single) . rConsist) 573-- test (multProp2 6 . (single *** single) . rConsist)
442 test (multProp2 6 . (single *** single) . cConsist) 574-- test (multProp2 6 . (single *** single) . cConsist)
443 putStrLn "------ sub-trans" 575 putStrLn "------ sub-trans"
444 test (subProp . rM) 576 test (subProp . rM)
445 test (subProp . cM) 577 test (subProp . cM)
@@ -455,9 +587,12 @@ runTests n = do
455 putStrLn "------ luSolve" 587 putStrLn "------ luSolve"
456 test (linearSolveProp (luSolve.luPacked) . rSqWC) 588 test (linearSolveProp (luSolve.luPacked) . rSqWC)
457 test (linearSolveProp (luSolve.luPacked) . cSqWC) 589 test (linearSolveProp (luSolve.luPacked) . cSqWC)
590 putStrLn "------ ldlSolve"
591 test (linearSolvePropH (ldlSolve.ldlPacked) . rSymWC)
592 test (linearSolvePropH (ldlSolve.ldlPacked) . cSymWC)
458 putStrLn "------ cholSolve" 593 putStrLn "------ cholSolve"
459 test (linearSolveProp (cholSolve.chol) . rPosDef) 594 test (linearSolveProp (cholSolve.chol.trustSym) . rPosDef)
460 test (linearSolveProp (cholSolve.chol) . cPosDef) 595 test (linearSolveProp (cholSolve.chol.trustSym) . cPosDef)
461 putStrLn "------ luSolveLS" 596 putStrLn "------ luSolveLS"
462 test (linearSolveProp linearSolveLS . rSqWC) 597 test (linearSolveProp linearSolveLS . rSqWC)
463 test (linearSolveProp linearSolveLS . cSqWC) 598 test (linearSolveProp linearSolveLS . cSqWC)
@@ -472,16 +607,16 @@ runTests n = do
472 putStrLn "------ svd" 607 putStrLn "------ svd"
473 test (svdProp1 . rM) 608 test (svdProp1 . rM)
474 test (svdProp1 . cM) 609 test (svdProp1 . cM)
475 test (svdProp1a svdR) 610 test (svdProp1a svd . rM)
476 test (svdProp1a svdC) 611 test (svdProp1a svd . cM)
477 test (svdProp1a svdRd) 612-- test (svdProp1a svdRd)
478 test (svdProp1b svdR) 613 test (svdProp1b svd . rM)
479 test (svdProp1b svdC) 614 test (svdProp1b svd . cM)
480 test (svdProp1b svdRd) 615-- test (svdProp1b svdRd)
481 test (svdProp2 thinSVDR) 616 test (svdProp2 thinSVD . rM)
482 test (svdProp2 thinSVDC) 617 test (svdProp2 thinSVD . cM)
483 test (svdProp2 thinSVDRd) 618-- test (svdProp2 thinSVDRd)
484 test (svdProp2 thinSVDCd) 619-- test (svdProp2 thinSVDCd)
485 test (svdProp3 . rM) 620 test (svdProp3 . rM)
486 test (svdProp3 . cM) 621 test (svdProp3 . cM)
487 test (svdProp4 . rM) 622 test (svdProp4 . rM)
@@ -492,12 +627,12 @@ runTests n = do
492 test (svdProp6b) 627 test (svdProp6b)
493 test (svdProp7 . rM) 628 test (svdProp7 . rM)
494 test (svdProp7 . cM) 629 test (svdProp7 . cM)
495 putStrLn "------ svdCd" 630-- putStrLn "------ svdCd"
496#ifdef NOZGESDD 631#ifdef NOZGESDD
497 putStrLn "Omitted" 632-- putStrLn "Omitted"
498#else 633#else
499 test (svdProp1a svdCd) 634-- test (svdProp1a svdCd)
500 test (svdProp1b svdCd) 635-- test (svdProp1b svdCd)
501#endif 636#endif
502 putStrLn "------ eig" 637 putStrLn "------ eig"
503 test (eigSHProp . rHer) 638 test (eigSHProp . rHer)
@@ -515,10 +650,10 @@ runTests n = do
515 test (qrProp . rM) 650 test (qrProp . rM)
516 test (qrProp . cM) 651 test (qrProp . cM)
517 test (rqProp . rM) 652 test (rqProp . rM)
518 test (rqProp . cM) 653-- test (rqProp . cM)
519 test (rqProp1 . cM) 654 test (rqProp1 . cM)
520 test (rqProp2 . cM) 655 test (rqProp2 . cM)
521 test (rqProp3 . cM) 656-- test (rqProp3 . cM)
522 putStrLn "------ hess" 657 putStrLn "------ hess"
523 test (hessProp . rSq) 658 test (hessProp . rSq)
524 test (hessProp . cSq) 659 test (hessProp . cSq)
@@ -528,8 +663,8 @@ runTests n = do
528 putStrLn "------ chol" 663 putStrLn "------ chol"
529 test (cholProp . rPosDef) 664 test (cholProp . rPosDef)
530 test (cholProp . cPosDef) 665 test (cholProp . cPosDef)
531 test (exactProp . rPosDef) 666-- test (exactProp . rPosDef)
532 test (exactProp . cPosDef) 667-- test (exactProp . cPosDef)
533 putStrLn "------ expm" 668 putStrLn "------ expm"
534 test (expmDiagProp . complex. rSqWC) 669 test (expmDiagProp . complex. rSqWC)
535 test (expmDiagProp . cSqWC) 670 test (expmDiagProp . cSqWC)
@@ -539,12 +674,12 @@ runTests n = do
539 test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM)) 674 test (\u -> sin u ** 2 + cos u ** 2 |~| (1::RM))
540 test (\u -> cos u * tan u |~| sin (u::RM)) 675 test (\u -> cos u * tan u |~| sin (u::RM))
541 test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary 676 test $ (\u -> cos u * tan u |~| sin (u::CM)) . liftMatrix makeUnitary
542 putStrLn "------ vector operations - Float" 677-- putStrLn "------ vector operations - Float"
543 test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM)) 678-- test (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::FM))
544 test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary 679-- test $ (\u -> sin u ^ 2 + cos u ^ 2 |~~| (1::ZM)) . liftMatrix makeUnitary
545 test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM)) 680-- test (\u -> sin u ** 2 + cos u ** 2 |~~| (1::FM))
546 test (\u -> cos u * tan u |~~| sin (u::FM)) 681-- test (\u -> cos u * tan u |~~| sin (u::FM))
547 test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary 682-- test $ (\u -> cos u * tan u |~~| sin (u::ZM)) . liftMatrix makeUnitary
548 putStrLn "------ read . show" 683 putStrLn "------ read . show"
549 test (\m -> (m::RM) == read (show m)) 684 test (\m -> (m::RM) == read (show m))
550 test (\m -> (m::CM) == read (show m)) 685 test (\m -> (m::CM) == read (show m))
@@ -562,8 +697,8 @@ runTests n = do
562 , utest "expm1" (expmTest1) 697 , utest "expm1" (expmTest1)
563 , utest "expm2" (expmTest2) 698 , utest "expm2" (expmTest2)
564 , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) 699 , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)
565 , utest "arith2" $ ((scalar (1+i) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*i-51) :: CM) 700 , utest "arith2" $ ((scalar (1+iC) * ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| ( scalar (140*iC-51) :: CM)
566 , utest "arith3" $ exp (scalar i * ones(10,10)*pi) + 1 |~| 0 701 , utest "arith3" $ exp (scalar iC * ones(10,10)*pi) + 1 |~| 0
567 , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3] 702 , utest "<\\>" $ (3><2) [2,0,0,3,1,1::Double] <\> 3|>[4,9,5] |~| 2|>[2,3]
568-- , utest "gamma" (gamma 5 == 24.0) 703-- , utest "gamma" (gamma 5 == 24.0)
569-- , besselTest 704-- , besselTest
@@ -571,10 +706,10 @@ runTests n = do
571 , utest "randomGaussian" randomTestGaussian 706 , utest "randomGaussian" randomTestGaussian
572 , utest "randomUniform" randomTestUniform 707 , utest "randomUniform" randomTestUniform
573 , utest "buildVector/Matrix" $ 708 , utest "buildVector/Matrix" $
574 complex (10 |> [0::Double ..]) == buildVector 10 fromIntegral 709 complex (10 |> [0::Double ..]) == build 10 id
575 && ident 5 == buildMatrix 5 5 (\(r,c) -> if r==c then 1::Double else 0) 710 && ident 5 == build (5,5) (\r c -> if r==c then 1::Double else 0)
576 , utest "rank" $ rank ((2><3)[1,0,0,1,5*eps,0]) == 1 711 , utest "rank" $ rank ((2><3)[1,0,0,1,5*peps,0::Double]) == 1
577 && rank ((2><3)[1,0,0,1,7*eps,0]) == 2 712 && rank ((2><3)[1,0,0,1,7*peps,0::Double]) == 2
578 , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM) 713 , utest "block" $ fromBlocks [[ident 3,0],[0,ident 4]] == (ident 7 :: CM)
579 , mbCholTest 714 , mbCholTest
580 , utest "offset" offsetTest 715 , utest "offset" offsetTest
@@ -588,21 +723,23 @@ runTests n = do
588 , conformTest 723 , conformTest
589 , accumTest 724 , accumTest
590 , convolutionTest 725 , convolutionTest
591 , kroneckerTest
592 , sparseTest 726 , sparseTest
593 , staticTest 727 , staticTest
728 , intTest
729 , modularTest
730 , sliceTest
594 ] 731 ]
595 when (errors c + failures c > 0) exitFailure 732 when (errors c + failures c > 0) exitFailure
596 return () 733 return ()
597 734
598 735
599-- single precision approximate equality 736-- single precision approximate equality
600infixl 4 |~~| 737-- infixl 4 |~~|
601a |~~| b = a :~6~: b 738-- a |~~| b = a :~6~: b
602 739
603makeUnitary v | realPart n > 1 = v / scalar n 740makeUnitary v | realPart n > 1 = v / scalar n
604 | otherwise = v 741 | otherwise = v
605 where n = sqrt (v <.> v) 742 where n = sqrt (v `dot` v)
606 743
607-- -- | Some additional tests on big matrices. They take a few minutes. 744-- -- | Some additional tests on big matrices. They take a few minutes.
608-- runBigTests :: IO () 745-- runBigTests :: IO ()
@@ -625,6 +762,8 @@ runBenchmarks = do
625 mkVecBench 762 mkVecBench
626 multBench 763 multBench
627 cholBench 764 cholBench
765 luBench
766 luBench_2
628 svdBench 767 svdBench
629 eigBench 768 eigBench
630 putStrLn "" 769 putStrLn ""
@@ -668,9 +807,9 @@ manyvec5 xs = sumElements $ fromRows $ map (\x -> vec3 x (x**2) (x**3)) xs
668 807
669 808
670manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs 809manyvec2 xs = sum $ map (\x -> sqrt(x^2 + (x**2)^2 +(x**3)^2)) xs
671manyvec3 xs = sum $ map (pnorm PNorm2 . (\x -> fromList [x,x**2,x**3])) xs 810manyvec3 xs = sum $ map (norm_2 . (\x -> fromList [x,x**2,x**3])) xs
672 811
673manyvec4 xs = sum $ map (pnorm PNorm2 . (\x -> vec3 x (x**2) (x**3))) xs 812manyvec4 xs = sum $ map (norm_2 . (\x -> vec3 x (x**2) (x**3))) xs
674 813
675vec3 :: Double -> Double -> Double -> Vector Double 814vec3 :: Double -> Double -> Double -> Vector Double
676vec3 a b c = runSTVector $ do 815vec3 a b c = runSTVector $ do
@@ -695,11 +834,11 @@ mkVecBench = do
695 834
696subBench = do 835subBench = do
697 putStrLn "" 836 putStrLn ""
698 let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (dim v -1) v)) 837 let g = foldl1' (.) (replicate (10^5) (\v -> subVector 1 (size v -1) v))
699 time "0.1M subVector " (g (konst 1 (1+10^5) :: Vector Double) @> 0) 838 time "0.1M subVector " (g (konst 1 (1+10^5) :: Vector Double) ! 0)
700 let f = foldl1' (.) (replicate (10^5) (fromRows.toRows)) 839 let f = foldl1' (.) (replicate (10^5) (fromRows.toRows))
701 time "subVector-join 3" (f (ident 3 :: Matrix Double) @@>(0,0)) 840 time "subVector-join 3" (f (ident 3 :: Matrix Double) `atIndex` (0,0))
702 time "subVector-join 10" (f (ident 10 :: Matrix Double) @@>(0,0)) 841 time "subVector-join 10" (f (ident 10 :: Matrix Double) `atIndex` (0,0))
703 842
704-------------------------------- 843--------------------------------
705 844
@@ -724,10 +863,10 @@ multBench = do
724 863
725eigBench = do 864eigBench = do
726 let m = reshape 1000 (randomVector 777 Uniform (1000*1000)) 865 let m = reshape 1000 (randomVector 777 Uniform (1000*1000))
727 s = m + trans m 866 s = m + tr m
728 m `seq` s `seq` putStrLn "" 867 m `seq` s `seq` putStrLn ""
729 time "eigenvalues symmetric 1000x1000" (eigenvaluesSH' m) 868 time "eigenvalues symmetric 1000x1000" (eigenvaluesSH (trustSym m))
730 time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m) 869 time "eigenvectors symmetric 1000x1000" (snd $ eigSH (trustSym m))
731 time "eigenvalues general 1000x1000" (eigenvalues m) 870 time "eigenvalues general 1000x1000" (eigenvalues m)
732 time "eigenvectors general 1000x1000" (snd $ eig m) 871 time "eigenvectors general 1000x1000" (snd $ eig m)
733 872
@@ -736,7 +875,7 @@ eigBench = do
736svdBench = do 875svdBench = do
737 let a = reshape 500 (randomVector 777 Uniform (3000*500)) 876 let a = reshape 500 (randomVector 777 Uniform (3000*500))
738 b = reshape 1000 (randomVector 777 Uniform (1000*1000)) 877 b = reshape 1000 (randomVector 777 Uniform (1000*1000))
739 fv (_,_,v) = v@@>(0,0) 878 fv (_,_,v) = v `atIndex` (0,0)
740 a `seq` b `seq` putStrLn "" 879 a `seq` b `seq` putStrLn ""
741 time "singular values 3000x500" (singularValues a) 880 time "singular values 3000x500" (singularValues a)
742 time "thin svd 3000x500" (fv $ thinSVD a) 881 time "thin svd 3000x500" (fv $ thinSVD a)
@@ -748,26 +887,28 @@ svdBench = do
748 887
749solveBenchN n = do 888solveBenchN n = do
750 let x = uniformSample 777 (2*n) (replicate n (-1,1)) 889 let x = uniformSample 777 (2*n) (replicate n (-1,1))
751 a = trans x <> x 890 a = tr x <> x
752 b = asColumn $ randomVector 666 Uniform n 891 b = asColumn $ randomVector 666 Uniform n
753 a `seq` b `seq` putStrLn "" 892 a `seq` b `seq` putStrLn ""
754 time ("svd solve " ++ show n) (linearSolveSVD a b) 893 time ("svd solve " ++ show n) (linearSolveSVD a b)
755 time (" ls solve " ++ show n) (linearSolveLS a b) 894 time (" ls solve " ++ show n) (linearSolveLS a b)
756 time (" solve " ++ show n) (linearSolve a b) 895 time (" solve " ++ show n) (linearSolve a b)
757 time ("cholSolve " ++ show n) (cholSolve (chol a) b) 896-- time (" LU solve " ++ show n) (luSolve (luPacked a) b)
897 time ("LDL solve " ++ show n) (ldlSolve (ldlPacked (trustSym a)) b)
898 time ("cholSolve " ++ show n) (cholSolve (chol $ trustSym a) b)
758 899
759solveBench = do 900solveBench = do
760 solveBenchN 500 901 solveBenchN 500
761 solveBenchN 1000 902 solveBenchN 1000
762 -- solveBenchN 1500 903 solveBenchN 1500
763 904
764-------------------------------- 905--------------------------------
765 906
766cholBenchN n = do 907cholBenchN n = do
767 let x = uniformSample 777 (2*n) (replicate n (-1,1)) 908 let x = uniformSample 777 (2*n) (replicate n (-1,1))
768 a = trans x <> x 909 a = tr x <> x
769 a `seq` putStr "" 910 a `seq` putStr ""
770 time ("chol " ++ show n) (chol a) 911 time ("chol " ++ show n) (chol $ trustSym a)
771 912
772cholBench = do 913cholBench = do
773 putStrLn "" 914 putStrLn ""
@@ -776,3 +917,32 @@ cholBench = do
776 cholBenchN 300 917 cholBenchN 300
777-- cholBenchN 150 918-- cholBenchN 150
778-- cholBenchN 50 919-- cholBenchN 50
920
921--------------------------------------------------------------------------------
922
923luBenchN f n x msg = do
924 let m = diagRect 1 (fromList (replicate n x)) n n
925 m `seq` putStr ""
926 time (msg ++ " "++ show n) (rnf $ f m)
927
928luBench = do
929 putStrLn ""
930 luBenchN luPacked 1000 (5::R) "luPacked Double "
931 luBenchN luPacked' 1000 (5::R) "luPacked' Double "
932 luBenchN luPacked' 1000 (5::Mod 9973 I) "luPacked' I mod 9973"
933 luBenchN luPacked' 1000 (5::Mod 9973 Z) "luPacked' Z mod 9973"
934
935luBenchN_2 f g n x msg = do
936 let m = diagRect 1 (fromList (replicate n x)) n n
937 b = flipud m
938 m `seq` b `seq` putStr ""
939 time (msg ++ " "++ show n) (f (g m) b)
940
941luBench_2 = do
942 putStrLn ""
943 luBenchN_2 luSolve luPacked 500 (5::R) "luSolve .luPacked Double "
944 luBenchN_2 luSolve' luPacked' 500 (5::R) "luSolve'.luPacked' Double "
945 luBenchN_2 luSolve' luPacked' 500 (5::Mod 9973 I) "luSolve'.luPacked' I mod 9973"
946 luBenchN_2 luSolve' luPacked' 500 (5::Mod 9973 Z) "luSolve'.luPacked' Z mod 9973"
947
948
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 53fc4d2..3d5441d 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -1,5 +1,4 @@
1{-# LANGUAGE FlexibleContexts, UndecidableInstances, CPP, FlexibleInstances #-} 1{-# LANGUAGE FlexibleContexts, UndecidableInstances, FlexibleInstances #-}
2{-# OPTIONS_GHC -fno-warn-unused-imports #-}
3----------------------------------------------------------------------------- 2-----------------------------------------------------------------------------
4{- | 3{- |
5Module : Numeric.LinearAlgebra.Tests.Instances 4Module : Numeric.LinearAlgebra.Tests.Instances
@@ -15,9 +14,9 @@ Arbitrary instances for vectors, matrices.
15module Numeric.LinearAlgebra.Tests.Instances( 14module Numeric.LinearAlgebra.Tests.Instances(
16 Sq(..), rSq,cSq, 15 Sq(..), rSq,cSq,
17 Rot(..), rRot,cRot, 16 Rot(..), rRot,cRot,
18 Her(..), rHer,cHer, 17 rHer,cHer,
19 WC(..), rWC,cWC, 18 WC(..), rWC,cWC,
20 SqWC(..), rSqWC, cSqWC, 19 SqWC(..), rSqWC, cSqWC, rSymWC, cSymWC,
21 PosDef(..), rPosDef, cPosDef, 20 PosDef(..), rPosDef, cPosDef,
22 Consistent(..), rConsist, cConsist, 21 Consistent(..), rConsist, cConsist,
23 RM,CM, rM,cM, 22 RM,CM, rM,cM,
@@ -26,15 +25,11 @@ module Numeric.LinearAlgebra.Tests.Instances(
26 25
27import System.Random 26import System.Random
28 27
29import Numeric.LinearAlgebra 28import Numeric.LinearAlgebra.HMatrix hiding (vector)
30import Numeric.LinearAlgebra.Devel
31import Numeric.Container
32import Control.Monad(replicateM) 29import Control.Monad(replicateM)
33import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 30import Test.QuickCheck(Arbitrary,arbitrary,choose,vector,sized,shrink)
34 ,sized,classify,Testable,Property 31
35 ,quickCheckWith,maxSize,stdArgs,shrink)
36 32
37#if MIN_VERSION_QuickCheck(2,0,0)
38shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]] 33shrinkListElementwise :: (Arbitrary a) => [a] -> [[a]]
39shrinkListElementwise [] = [] 34shrinkListElementwise [] = []
40shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ] 35shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ]
@@ -42,25 +37,6 @@ shrinkListElementwise (x:xs) = [ y:xs | y <- shrink x ]
42 37
43shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)] 38shrinkPair :: (Arbitrary a, Arbitrary b) => (a,b) -> [(a,b)]
44shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ] 39shrinkPair (a,b) = [ (a,x) | x <- shrink b ] ++ [ (x,b) | x <- shrink a ]
45#endif
46
47#if MIN_VERSION_QuickCheck(2,1,1)
48#else
49instance (Arbitrary a, RealFloat a) => Arbitrary (Complex a) where
50 arbitrary = do
51 re <- arbitrary
52 im <- arbitrary
53 return (re :+ im)
54
55#if MIN_VERSION_QuickCheck(2,0,0)
56 shrink (re :+ im) =
57 [ u :+ v | (u,v) <- shrinkPair (re,im) ]
58#else
59 -- this has been moved to the 'Coarbitrary' class in QuickCheck 2
60 coarbitrary = undefined
61#endif
62
63#endif
64 40
65chooseDim = sized $ \m -> choose (1,max 1 m) 41chooseDim = sized $ \m -> choose (1,max 1 m)
66 42
@@ -68,15 +44,9 @@ instance (Field a, Arbitrary a) => Arbitrary (Vector a) where
68 arbitrary = do m <- chooseDim 44 arbitrary = do m <- chooseDim
69 l <- vector m 45 l <- vector m
70 return $ fromList l 46 return $ fromList l
71
72#if MIN_VERSION_QuickCheck(2,0,0)
73 -- shrink any one of the components 47 -- shrink any one of the components
74 shrink = map fromList . shrinkListElementwise . toList 48 shrink = map fromList . shrinkListElementwise . toList
75 49
76#else
77 coarbitrary = undefined
78#endif
79
80instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where 50instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
81 arbitrary = do 51 arbitrary = do
82 m <- chooseDim 52 m <- chooseDim
@@ -84,16 +54,11 @@ instance (Element a, Arbitrary a) => Arbitrary (Matrix a) where
84 l <- vector (m*n) 54 l <- vector (m*n)
85 return $ (m><n) l 55 return $ (m><n) l
86 56
87#if MIN_VERSION_QuickCheck(2,0,0)
88 -- shrink any one of the components 57 -- shrink any one of the components
89 shrink a = map (rows a >< cols a) 58 shrink a = map (rows a >< cols a)
90 . shrinkListElementwise 59 . shrinkListElementwise
91 . concat . toLists 60 . concat . toLists
92 $ a 61 $ a
93#else
94 coarbitrary = undefined
95#endif
96
97 62
98-- a square matrix 63-- a square matrix
99newtype (Sq a) = Sq (Matrix a) deriving Show 64newtype (Sq a) = Sq (Matrix a) deriving Show
@@ -103,11 +68,7 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
103 l <- vector (n*n) 68 l <- vector (n*n)
104 return $ Sq $ (n><n) l 69 return $ Sq $ (n><n) l
105 70
106#if MIN_VERSION_QuickCheck(2,0,0)
107 shrink (Sq a) = [ Sq b | b <- shrink a ] 71 shrink (Sq a) = [ Sq b | b <- shrink a ]
108#else
109 coarbitrary = undefined
110#endif
111 72
112 73
113-- a unitary matrix 74-- a unitary matrix
@@ -118,24 +79,14 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
118 let (q,_) = qr m 79 let (q,_) = qr m
119 return (Rot q) 80 return (Rot q)
120 81
121#if MIN_VERSION_QuickCheck(2,0,0)
122#else
123 coarbitrary = undefined
124#endif
125
126 82
127-- a complex hermitian or real symmetric matrix 83-- a complex hermitian or real symmetric matrix
128newtype (Her a) = Her (Matrix a) deriving Show 84instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Herm a) where
129instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
130 arbitrary = do 85 arbitrary = do
131 Sq m <- arbitrary 86 Sq m <- arbitrary
132 let m' = m/2 87 let m' = m/2
133 return $ Her (m' + ctrans m') 88 return $ sym m'
134 89
135#if MIN_VERSION_QuickCheck(2,0,0)
136#else
137 coarbitrary = undefined
138#endif
139 90
140class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a 91class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
141instance ArbitraryField Double 92instance ArbitraryField Double
@@ -144,7 +95,7 @@ instance ArbitraryField (Complex Double)
144 95
145-- a well-conditioned general matrix (the singular values are between 1 and 100) 96-- a well-conditioned general matrix (the singular values are between 1 and 100)
146newtype (WC a) = WC (Matrix a) deriving Show 97newtype (WC a) = WC (Matrix a) deriving Show
147instance (ArbitraryField a) => Arbitrary (WC a) where 98instance (Numeric a, ArbitraryField a) => Arbitrary (WC a) where
148 arbitrary = do 99 arbitrary = do
149 m <- arbitrary 100 m <- arbitrary
150 let (u,_,v) = svd m 101 let (u,_,v) = svd m
@@ -153,48 +104,33 @@ instance (ArbitraryField a) => Arbitrary (WC a) where
153 n = min r c 104 n = min r c
154 sv' <- replicateM n (choose (1,100)) 105 sv' <- replicateM n (choose (1,100))
155 let s = diagRect 0 (fromList sv') r c 106 let s = diagRect 0 (fromList sv') r c
156 return $ WC (u `mXm` real s `mXm` trans v) 107 return $ WC (u <> real s <> tr v)
157
158#if MIN_VERSION_QuickCheck(2,0,0)
159#else
160 coarbitrary = undefined
161#endif
162 108
163 109
164-- a well-conditioned square matrix (the singular values are between 1 and 100) 110-- a well-conditioned square matrix (the singular values are between 1 and 100)
165newtype (SqWC a) = SqWC (Matrix a) deriving Show 111newtype (SqWC a) = SqWC (Matrix a) deriving Show
166instance (ArbitraryField a) => Arbitrary (SqWC a) where 112instance (ArbitraryField a, Numeric a) => Arbitrary (SqWC a) where
167 arbitrary = do 113 arbitrary = do
168 Sq m <- arbitrary 114 Sq m <- arbitrary
169 let (u,_,v) = svd m 115 let (u,_,v) = svd m
170 n = rows m 116 n = rows m
171 sv' <- replicateM n (choose (1,100)) 117 sv' <- replicateM n (choose (1,100))
172 let s = diag (fromList sv') 118 let s = diag (fromList sv')
173 return $ SqWC (u `mXm` real s `mXm` trans v) 119 return $ SqWC (u <> real s <> tr v)
174
175#if MIN_VERSION_QuickCheck(2,0,0)
176#else
177 coarbitrary = undefined
178#endif
179 120
180 121
181-- a positive definite square matrix (the eigenvalues are between 0 and 100) 122-- a positive definite square matrix (the eigenvalues are between 0 and 100)
182newtype (PosDef a) = PosDef (Matrix a) deriving Show 123newtype (PosDef a) = PosDef (Matrix a) deriving Show
183instance (ArbitraryField a, Num (Vector a)) 124instance (Numeric a, ArbitraryField a, Num (Vector a))
184 => Arbitrary (PosDef a) where 125 => Arbitrary (PosDef a) where
185 arbitrary = do 126 arbitrary = do
186 Her m <- arbitrary 127 m <- arbitrary
187 let (_,v) = eigSH m 128 let (_,v) = eigSH m
188 n = rows m 129 n = rows (unSym m)
189 l <- replicateM n (choose (0,100)) 130 l <- replicateM n (choose (0,100))
190 let s = diag (fromList l) 131 let s = diag (fromList l)
191 p = v `mXm` real s `mXm` ctrans v 132 p = v <> real s <> tr v
192 return $ PosDef (0.5 * p + 0.5 * ctrans p) 133 return $ PosDef (0.5 * p + 0.5 * tr p)
193
194#if MIN_VERSION_QuickCheck(2,0,0)
195#else
196 coarbitrary = undefined
197#endif
198 134
199 135
200-- a pair of matrices that can be multiplied 136-- a pair of matrices that can be multiplied
@@ -208,11 +144,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Consistent a) where
208 lb <- vector (k*m) 144 lb <- vector (k*m)
209 return $ Consistent ((n><k) la, (k><m) lb) 145 return $ Consistent ((n><k) la, (k><m) lb)
210 146
211#if MIN_VERSION_QuickCheck(2,0,0)
212 shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ] 147 shrink (Consistent (x,y)) = [ Consistent (u,v) | (u,v) <- shrinkPair (x,y) ]
213#else
214 coarbitrary = undefined
215#endif
216 148
217 149
218 150
@@ -228,8 +160,8 @@ fM m = m :: FM
228zM m = m :: ZM 160zM m = m :: ZM
229 161
230 162
231rHer (Her m) = m :: RM 163rHer m = unSym m :: RM
232cHer (Her m) = m :: CM 164cHer m = unSym m :: CM
233 165
234rRot (Rot m) = m :: RM 166rRot (Rot m) = m :: RM
235cRot (Rot m) = m :: CM 167cRot (Rot m) = m :: CM
@@ -243,6 +175,9 @@ cWC (WC m) = m :: CM
243rSqWC (SqWC m) = m :: RM 175rSqWC (SqWC m) = m :: RM
244cSqWC (SqWC m) = m :: CM 176cSqWC (SqWC m) = m :: CM
245 177
178rSymWC (SqWC m) = sym m :: Herm R
179cSymWC (SqWC m) = sym m :: Herm C
180
246rPosDef (PosDef m) = m :: RM 181rPosDef (PosDef m) = m :: RM
247cPosDef (PosDef m) = m :: CM 182cPosDef (PosDef m) = m :: CM
248 183
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index a5c37f4..046644f 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -1,5 +1,4 @@
1{-# LANGUAGE CPP, FlexibleContexts #-} 1{-# LANGUAGE FlexibleContexts #-}
2{-# OPTIONS_GHC -fno-warn-unused-imports #-}
3{-# LANGUAGE TypeFamilies #-} 2{-# LANGUAGE TypeFamilies #-}
4 3
5----------------------------------------------------------------------------- 4-----------------------------------------------------------------------------
@@ -15,7 +14,7 @@ Testing properties.
15-} 14-}
16 15
17module Numeric.LinearAlgebra.Tests.Properties ( 16module Numeric.LinearAlgebra.Tests.Properties (
18 dist, (|~|), (~~), (~:), Aprox((:~)), 17 dist, (|~|), (~~), (~:), Aprox((:~)), (~=),
19 zeros, ones, 18 zeros, ones,
20 square, 19 square,
21 unitary, 20 unitary,
@@ -29,7 +28,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
29 pinvProp, 28 pinvProp,
30 detProp, 29 detProp,
31 nullspaceProp, 30 nullspaceProp,
32 bugProp, 31-- bugProp,
33 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, 32 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
34 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, 33 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
35 eigProp, eigSHProp, eigProp2, eigSHProp2, 34 eigProp, eigSHProp, eigProp2, eigSHProp2,
@@ -40,23 +39,21 @@ module Numeric.LinearAlgebra.Tests.Properties (
40 expmDiagProp, 39 expmDiagProp,
41 multProp1, multProp2, 40 multProp1, multProp2,
42 subProp, 41 subProp,
43 linearSolveProp, linearSolveProp2 42 linearSolveProp, linearSolvePropH, linearSolveProp2
44) where 43) where
45 44
46import Numeric.Container 45import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary)
47import Numeric.LinearAlgebra --hiding (real,complex) 46import Test.QuickCheck
48import Numeric.LinearAlgebra.LAPACK 47
49import Debug.Trace 48(~=) :: Double -> Double -> Bool
50import Test.QuickCheck(Arbitrary,arbitrary,coarbitrary,choose,vector 49a ~= b = abs (a - b) < 1e-10
51 ,sized,classify,Testable,Property
52 ,quickCheckWith,maxSize,stdArgs,shrink)
53 50
54trivial :: Testable a => Bool -> a -> Property 51trivial :: Testable a => Bool -> a -> Property
55trivial = (`classify` "trivial") 52trivial = (`classify` "trivial")
56 53
57-- relative error 54-- relative error
58dist :: (Normed c t, Num (c t)) => c t -> c t -> Double 55dist :: (Num a, Normed a) => a -> a -> Double
59dist = relativeError Infinity 56dist = relativeError norm_Inf
60 57
61infixl 4 |~| 58infixl 4 |~|
62a |~| b = a :~10~: b 59a |~| b = a :~10~: b
@@ -73,11 +70,11 @@ a :~n~: b = dist a b < 10^^(-n)
73square m = rows m == cols m 70square m = rows m == cols m
74 71
75-- orthonormal columns 72-- orthonormal columns
76orthonormal m = ctrans m <> m |~| ident (cols m) 73orthonormal m = tr m <> m |~| ident (cols m)
77 74
78unitary m = square m && orthonormal m 75unitary m = square m && orthonormal m
79 76
80hermitian m = square m && m |~| ctrans m 77hermitian m = square m && m |~| tr m
81 78
82wellCond m = rcond m > 1/100 79wellCond m = rcond m > 1/100
83 80
@@ -85,12 +82,12 @@ positiveDefinite m = minimum (toList e) > 0
85 where (e,_v) = eigSH m 82 where (e,_v) = eigSH m
86 83
87upperTriang m = rows m == 1 || down == z 84upperTriang m = rows m == 1 || down == z
88 where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m)) 85 where down = fromList $ concat $ zipWith drop [1..] (toLists (tr m))
89 z = konst 0 (dim down) 86 z = konst 0 (size down)
90 87
91upperHessenberg m = rows m < 3 || down == z 88upperHessenberg m = rows m < 3 || down == z
92 where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m)) 89 where down = fromList $ concat $ zipWith drop [2..] (toLists (tr m))
93 z = konst 0 (dim down) 90 z = konst 0 (size down)
94 91
95zeros (r,c) = reshape c (konst 0 (r*c)) 92zeros (r,c) = reshape c (konst 0 (r*c))
96 93
@@ -118,81 +115,94 @@ detProp m = s d1 |~| s d2
118 s x = fromList [x] 115 s x = fromList [x]
119 116
120nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c) 117nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c)
121 && orthonormal (fromColumns nl)) 118 && orthonormal n)
122 where nl = nullspacePrec 1 m 119 where n = nullspaceSVD (Left (1*peps)) m (rightSV m)
123 n = fromColumns nl 120 nl = toColumns n
124 r = rows m 121 r = rows m
125 c = cols m - rank m 122 c = cols m - rank m
126 123
127------------------------------------------------------------------ 124------------------------------------------------------------------
128 125{-
129-- testcase for nonempty fpu stack 126-- testcase for nonempty fpu stack
130-- uncommenting unitary' signature eliminates the problem 127-- uncommenting unitary' signature eliminates the problem
131bugProp m = m |~| u <> real d <> trans v && unitary' u && unitary' v 128bugProp m = m |~| u <> real d <> tr v && unitary' u && unitary' v
132 where (u,d,v) = fullSVD m 129 where (u,d,v) = svd m
133 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool 130 -- unitary' :: (Num (Vector t), Field t) => Matrix t -> Bool
134 unitary' a = unitary a 131 unitary' a = unitary a
135 132-}
136------------------------------------------------------------------ 133------------------------------------------------------------------
137 134
138-- fullSVD 135-- fullSVD
139svdProp1 m = m |~| u <> real d <> trans v && unitary u && unitary v 136svdProp1 m = m |~| u <> real d <> tr v && unitary u && unitary v
140 where (u,d,v) = fullSVD m 137 where
138 (u,s,v) = svd m
139 d = diagRect 0 s (rows m) (cols m)
141 140
142svdProp1a svdfun m = m |~| u <> real d <> trans v && unitary u && unitary v where 141svdProp1a svdfun m = m |~| u <> real d <> tr v && unitary u && unitary v
142 where
143 (u,s,v) = svdfun m 143 (u,s,v) = svdfun m
144 d = diagRect 0 s (rows m) (cols m) 144 d = diagRect 0 s (rows m) (cols m)
145 145
146svdProp1b svdfun m = unitary u && unitary v where 146svdProp1b svdfun m = unitary u && unitary v
147 where
147 (u,_,v) = svdfun m 148 (u,_,v) = svdfun m
148 149
149-- thinSVD 150-- thinSVD
150svdProp2 thinSVDfun m = m |~| u <> diag (real s) <> trans v && orthonormal u && orthonormal v && dim s == min (rows m) (cols m) 151svdProp2 thinSVDfun m
151 where (u,s,v) = thinSVDfun m 152 = m |~| u <> diag (real s) <> tr v
153 && orthonormal u && orthonormal v
154 && size s == min (rows m) (cols m)
155 where
156 (u,s,v) = thinSVDfun m
152 157
153-- compactSVD 158-- compactSVD
154svdProp3 m = (m |~| u <> real (diag s) <> trans v 159svdProp3 m = (m |~| u <> real (diag s) <> tr v
155 && orthonormal u && orthonormal v) 160 && orthonormal u && orthonormal v)
156 where (u,s,v) = compactSVD m 161 where
162 (u,s,v) = compactSVD m
157 163
158svdProp4 m' = m |~| u <> real (diag s) <> trans v 164svdProp4 m' = m |~| u <> real (diag s) <> tr v
159 && orthonormal u && orthonormal v 165 && orthonormal u && orthonormal v
160 && (dim s == r || r == 0 && dim s == 1) 166 && (size s == r || r == 0 && size s == 1)
161 where (u,s,v) = compactSVD m 167 where
162 m = fromBlocks [[m'],[m']] 168 (u,s,v) = compactSVD m
163 r = rank m' 169 m = fromBlocks [[m'],[m']]
164 170 r = rank m'
165svdProp5a m = all (s1|~|) [s2,s3,s4,s5,s6] where 171
166 s1 = svR m 172svdProp5a m = all (s1|~|) [s3,s5] where
167 s2 = svRd m 173 s1 = singularValues (m :: Matrix Double)
168 (_,s3,_) = svdR m 174-- s2 = svRd m
169 (_,s4,_) = svdRd m 175 (_,s3,_) = svd m
170 (_,s5,_) = thinSVDR m 176-- (_,s4,_) = svdRd m
171 (_,s6,_) = thinSVDRd m 177 (_,s5,_) = thinSVD m
172 178-- (_,s6,_) = thinSVDRd m
173svdProp5b m = all (s1|~|) [s2,s3,s4,s5,s6] where 179
174 s1 = svC m 180svdProp5b m = all (s1|~|) [s3,s5] where
175 s2 = svCd m 181 s1 = singularValues (m :: Matrix (Complex Double))
176 (_,s3,_) = svdC m 182-- s2 = svCd m
177 (_,s4,_) = svdCd m 183 (_,s3,_) = svd m
178 (_,s5,_) = thinSVDC m 184-- (_,s4,_) = svdCd m
179 (_,s6,_) = thinSVDCd m 185 (_,s5,_) = thinSVD m
186-- (_,s6,_) = thinSVDCd m
180 187
181svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 188svdProp6a m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
182 where (u,s,v) = svdR m 189 where
183 (s',v') = rightSVR m 190 (u,s,v) = svd (m :: Matrix Double)
184 (u',s'') = leftSVR m 191 (s',v') = rightSV m
192 (u',s'') = leftSV m
185 193
186svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u' 194svdProp6b m = s |~| s' && v |~| v' && s |~| s'' && u |~| u'
187 where (u,s,v) = svdC m 195 where
188 (s',v') = rightSVC m 196 (u,s,v) = svd (m :: Matrix (Complex Double))
189 (u',s'') = leftSVC m 197 (s',v') = rightSV m
198 (u',s'') = leftSV m
190 199
191svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s''' 200svdProp7 m = s |~| s' && u |~| u' && v |~| v' && s |~| s'''
192 where (u,s,v) = svd m 201 where
193 (s',v') = rightSV m 202 (u,s,v) = svd m
194 (u',_s'') = leftSV m 203 (s',v') = rightSV m
195 s''' = singularValues m 204 (u',_s'') = leftSV m
205 s''' = singularValues m
196 206
197------------------------------------------------------------------ 207------------------------------------------------------------------
198 208
@@ -201,12 +211,12 @@ eigProp m = complex m <> v |~| v <> diag s
201 211
202eigSHProp m = m <> v |~| v <> real (diag s) 212eigSHProp m = m <> v |~| v <> real (diag s)
203 && unitary v 213 && unitary v
204 && m |~| v <> real (diag s) <> ctrans v 214 && m |~| v <> real (diag s) <> tr v
205 where (s, v) = eigSH m 215 where (s, v) = eigSH' m
206 216
207eigProp2 m = fst (eig m) |~| eigenvalues m 217eigProp2 m = fst (eig m) |~| eigenvalues m
208 218
209eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m 219eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m
210 220
211------------------------------------------------------------------ 221------------------------------------------------------------------
212 222
@@ -226,22 +236,22 @@ rqProp3 m = upperTriang' r
226 where (r,_q) = rq m 236 where (r,_q) = rq m
227 237
228upperTriang' r = upptr (rows r) (cols r) * r |~| r 238upperTriang' r = upptr (rows r) (cols r) * r |~| r
229 where upptr f c = buildMatrix f c $ \(r',c') -> if r'-t > c' then 0 else 1 239 where upptr f c = build (f,c) $ \r' c' -> if r'-t > c' then 0 else 1
230 where t = f-c 240 where t = fromIntegral (f-c)
231 241
232hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h 242hessProp m = m |~| p <> h <> tr p && unitary p && upperHessenberg h
233 where (p,h) = hess m 243 where (p,h) = hess m
234 244
235schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s 245schurProp1 m = m |~| u <> s <> tr u && unitary u && upperTriang s
236 where (u,s) = schur m 246 where (u,s) = schur m
237 247
238schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme 248schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
239 where (u,s) = schur m 249 where (u,s) = schur m
240 250
241cholProp m = m |~| ctrans c <> c && upperTriang c 251cholProp m = m |~| tr c <> c && upperTriang c
242 where c = chol m 252 where c = chol (trustSym m)
243 253
244exactProp m = chol m == chol (m+0) 254exactProp m = chol (trustSym m) == chol (trustSym (m+0))
245 255
246expmDiagProp m = expm (logm m) :~ 7 ~: complex m 256expmDiagProp m = expm (logm m) :~ 7 ~: complex m
247 where logm = matFunc log 257 where logm = matFunc log
@@ -252,14 +262,16 @@ mulH a b = fromLists [[ doth ai bj | bj <- toColumns b] | ai <- toRows a ]
252 262
253multProp1 p (a,b) = (a <> b) :~p~: (mulH a b) 263multProp1 p (a,b) = (a <> b) :~p~: (mulH a b)
254 264
255multProp2 p (a,b) = (ctrans (a <> b)) :~p~: (ctrans b <> ctrans a) 265multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)
256 266
257linearSolveProp f m = f m m |~| ident (rows m) 267linearSolveProp f m = f m m |~| ident (rows m)
258 268
269linearSolvePropH f m = f m (unSym m) |~| ident (rows (unSym m))
270
259linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b) 271linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
260 where q = min (rows a) (cols a) 272 where q = min (rows a) (cols a)
261 b = a <> x 273 b = a <> x
262 wc = rank a == q 274 wc = rank a == q
263 275
264subProp m = m == (trans . fromColumns . toRows) m 276subProp m = m == (conj . tr . fromColumns . toRows) m
265 277
generated by cgit v1.2.3 (git 2.39.1) at 2024-11-10 04:16:42 +0000