2 files changed, 127 insertions, 6 deletions

diff --git a/lib/Numeric/LinearAlgebra/Tests/Instances.hs b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
index 583143a..af486c8 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -14,10 +14,13 @@ Arbitrary instances for vectors, matrices.
 -}
 
 module Numeric.LinearAlgebra.Tests.Instances(
-    Sq(..),
-    Rot(..),
-    Her(..),
-    WC(..)
+    Sq(..),     rSq,cSq,
+    Rot(..),    rRot,cRot,
+    Her(..),    rHer,cHer,
+    WC(..),     rWC,cWC,
+    SqWC(..),   rSqWC, cSqWC,
+    PosDef(..), rPosDef, cPosDef,
+    RM,CM, rM,cM
 ) where
 
 import Numeric.LinearAlgebra
@@ -74,7 +77,7 @@ instance (Field a, Arbitrary a) => Arbitrary (Her a) where
         return $ Her (m' + ctrans m')
     coarbitrary = undefined
 
--- a well-conditioned matrix (the singular values are between 1 and 100)
+-- a well-conditioned general matrix (the singular values are between 1 and 100)
 newtype (WC a) = WC (Matrix a) deriving Show
 instance (Field a, Arbitrary a) => Arbitrary (WC a) where
     arbitrary = do
@@ -87,3 +90,53 @@ instance (Field a, Arbitrary a) => Arbitrary (WC a) where
         let s = diagRect (fromList sv) r c
         return $ WC (u <> real s <> trans v)
     coarbitrary = undefined
+
+-- a well-conditioned square matrix (the singular values are between 1 and 100)
+newtype (SqWC a) = SqWC (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (SqWC a) where
+    arbitrary = do
+        Sq m <- arbitrary
+        let (u,_,v) = svd m
+            n = rows m
+        sv <- replicateM n (choose (1,100))
+        let s = diag (fromList sv)
+        return $ SqWC (u <> real s <> trans v)
+    coarbitrary = undefined
+
+-- a positive definite square matrix (the eigenvalues are between 0 and 100)
+newtype (PosDef a) = PosDef (Matrix a) deriving Show
+instance (Field a, Arbitrary a) => Arbitrary (PosDef a) where
+    arbitrary = do
+        Her m <- arbitrary
+        let (_,v) = eigSH m
+            n = rows m
+        l <- replicateM n (choose (0,100))
+        let s = diag (fromList l)
+            p = v <> real s <> ctrans v
+        return $ PosDef (0.5 .* p + 0.5 .* ctrans p)
+    coarbitrary = undefined
+
+type RM = Matrix Double
+type CM = Matrix (Complex Double)
+
+rM m = m :: RM
+cM m = m :: CM
+
+rHer (Her m) = m :: RM
+cHer (Her m) = m :: CM
+
+rRot (Rot m) = m :: RM
+cRot (Rot m) = m :: CM
+
+rSq  (Sq m)  = m :: RM
+cSq  (Sq m)  = m :: CM
+
+rWC (WC m) = m :: RM
+cWC (WC m) = m :: CM
+
+rSqWC (SqWC m) = m :: RM
+cSqWC (SqWC m) = m :: CM
+
+rPosDef (PosDef m) = m :: RM
+cPosDef (PosDef m) = m :: CM
+
diff --git a/lib/Numeric/LinearAlgebra/Tests/Properties.hs b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
index 351615b..0317469 100644
--- a/lib/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/lib/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -19,6 +19,7 @@ where
 
 import Numeric.LinearAlgebra
 import Numeric.LinearAlgebra.Tests.Instances(Sq(..),Her(..),Rot(..))
+import Test.QuickCheck
 
 -- relative error
 dist :: (Normed t, Num t) => t -> t -> Double
@@ -53,7 +54,74 @@ degenerate m = rank m < min (rows m) (cols m)
 
 wellCond m = rcond m > 1/100
 
+positiveDefinite m = minimum (toList e) > 0
+    where (e,v) = eigSH m
+
+upperTriang m = rows m == 1 || down == z
+    where down = fromList $ concat $ zipWith drop [1..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
+upperHessenberg m = rows m < 3 || down == z
+    where down = fromList $ concat $ zipWith drop [2..] (toLists (ctrans m))
+          z = constant 0 (dim down)
+
+zeros (r,c) = reshape c (constant 0 (r*c))
+
+ones (r,c) = zeros (r,c) + 1
+
 -----------------------------------------------------
 
-luTest m = m |~| p <> l <> u && det p == s
+luProp m = m |~| p <> l <> u && det p == s
     where (l,u,p,s) = lu m
+
+invProp m = m <> inv m |~| ident (rows m)
+
+pinvProp m =  m <> p <> m |~| m
+           && p <> m <> p |~| p
+           && hermitian (m<>p)
+           && hermitian (p<>m)
+    where p = pinv m
+
+detProp m = s d1 |~| s d2
+    where d1 = det m
+          d2 = det' m * det q
+          det' m = product $ toList $ takeDiag r
+          (q,r) = qr m
+          s x = fromList [x]
+
+nullspaceProp m = null nl `trivial` (null nl || m <> n |~| zeros (r,c))
+    where nl = nullspacePrec 1 m
+          n = fromColumns nl
+          r = rows m
+          c = cols m - rank m
+
+svdProp1 m = u <> real d <> trans v |~| m
+          && unitary u && unitary v
+    where (u,d,v) = full svd m
+
+svdProp2 m = (m |~| 0) `trivial` ((m |~| 0) || u <> real (diag s) <> trans v |~| m)
+    where (u,s,v) = economy svd m
+
+eigProp m = complex m <> v |~| v <> diag s
+    where (s, v) = eig m
+
+eigSHProp m = m <> v |~| v <> real (diag s)
+              && unitary v
+              && m |~| v <> real (diag s) <> ctrans v
+    where (s, v) = eigSH m
+
+qrProp m = q <> r |~| m && unitary q && upperTriang r
+    where (q,r) = qr m
+
+hessProp m = m |~| p <> h <> ctrans p && unitary p && upperHessenberg h
+    where (p,h) = hess m
+
+schurProp1 m = m |~| u <> s <> ctrans u && unitary u && upperTriang s
+    where (u,s) = schur m
+
+schurProp2 m = m |~| u <> s <> ctrans u && unitary u && upperHessenberg s -- fixme
+    where (u,s) = schur m
+
+cholProp m = m |~| ctrans c <> c && upperTriang c
+    where c = chol m
+          pos = positiveDefinite m