4 files changed, 157 insertions, 39 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index 79cc64d..1de018c 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -142,7 +142,7 @@ eigSH m | m `equal` ctrans m = eigSH' m
 
 -- | Cholesky factorization of a positive definite hermitian or symmetric matrix using lapack's dpotrf or zportrf.
 --
--- If @c = chol m@ then @m == c \<> ctrans c@.
+-- If @c = chol m@ then @m == ctrans c \<> c@.
 chol :: Field t => Matrix t ->  Matrix t
 chol m | m `equal` ctrans m = cholSH m
        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs
index 534eb04..7ecf812 100644
--- a/lib/Numeric/LinearAlgebra/Tests.hs
+++ b/lib/Numeric/LinearAlgebra/Tests.hs
@@ -15,8 +15,7 @@ Some tests.
 module Numeric.LinearAlgebra.Tests(
   module Numeric.LinearAlgebra.Tests.Instances,
   module Numeric.LinearAlgebra.Tests.Properties,
-  qCheck,RM,CM, rM,cM, rHer,cHer,rRot,cRot,rSq,cSq,
-    runTests
+  qCheck, runTests
 --, runBigTests
 ) where
 
@@ -24,43 +23,41 @@ import Numeric.LinearAlgebra
 import Numeric.LinearAlgebra.Tests.Instances
 import Numeric.LinearAlgebra.Tests.Properties
 import Test.QuickCheck
-import Debug.Trace
+import Numeric.GSL(setErrorHandlerOff)
 
 qCheck n = check defaultConfig {configSize = const n}
 
-debug x = trace (show x) x
-
-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
-
 -- | It runs all the tests.
 runTests :: Int  -- ^ maximum dimension
          -> IO ()
 runTests n = do
-    qCheck n (hermitian . rHer)
-    qCheck n (hermitian . cHer)
-    qCheck n (unitary . rRot)
-    qCheck n (unitary . cRot)
-    qCheck n (wellCond . rWC)
-    qCheck n (wellCond . cWC)
-    --------------------------------
-    qCheck n (luTest . rM)
-    qCheck n (luTest . cM)
+    setErrorHandlerOff
+    let test p = qCheck n p
+    test (luProp    . rM)
+    test (luProp    . cM)
+    test (invProp   . rSqWC)
+    test (invProp   . cSqWC)
+    test (pinvProp  . rM)
+    test (pinvProp  . cM)
+    test (detProp   . cSqWC)
+    test (svdProp1  . rM)
+    test (svdProp1  . cM)
+    test (svdProp2  . rM)
+    test (svdProp2  . cM)
+    test (eigSHProp . rHer)
+    test (eigSHProp . cHer)
+    test (eigProp   . rSq)
+    test (eigProp   . cSq)
+    test (nullspaceProp . rM)
+    test (nullspaceProp . cM)
+    test (qrProp     . rM)
+    test (qrProp     . cM)
+    test (hessProp   . rSq)
+    test (hessProp   . cSq)
+    test (schurProp2 . rSq)
+    test (schurProp1 . cSq)
+    test (cholProp   . rPosDef)
+    test (cholProp   . cPosDef)
 
 -- | Some additional tests on big matrices. They take a few minutes.
 runBigTests :: IO ()
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