Her, LU, LDL, Linear, Additive

3 files changed, 41 insertions, 37 deletions

diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 2ff1580..30480d7 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -127,8 +127,8 @@ expmTest2 = expm nd2 :~15~: (2><2)
 mbCholTest = utest "mbCholTest" (ok1 && ok2) where
     m1 = (2><2) [2,5,5,8 :: Double]
     m2 = (2><2) [3,5,5,9 :: Complex Double]
-    ok1 = mbCholSH m1 == Nothing
-    ok2 = mbCholSH m2 == Just (chol m2)
+    ok1 = mbChol (trustSym m1) == Nothing
+    ok2 = mbChol (trustSym m2) == Just (chol $ trustSym m2)
 
 ---------------------------------------------------------------------
 
@@ -403,8 +403,8 @@ indexProp g f x = a1 == g a2 && a2 == a3 && b1 == g b2 && b2 == b3
 --------------------------------------------------------------------------------
 
 sliceTest = utest "slice test" $ and
-    [ testSlice chol  (gen 5 :: Matrix R)
-    , testSlice chol  (gen 5 :: Matrix C)
+    [ testSlice (chol . trustSym)  (gen 5 :: Matrix R)
+    , testSlice (chol . trustSym)  (gen 5 :: Matrix C)
     , testSlice qr    (rec :: Matrix R)
     , testSlice qr    (rec :: Matrix C)
     , testSlice hess  (agen 5 :: Matrix R)
@@ -420,12 +420,12 @@ sliceTest = utest "slice test" $ and
 
     , testSlice eig   (agen 5 :: Matrix R)
     , testSlice eig   (agen 5 :: Matrix C)
-    , testSlice eigSH (gen 5 :: Matrix R)
-    , testSlice eigSH (gen 5 :: Matrix C)
+    , testSlice (eigSH . trustSym) (gen 5 :: Matrix R)
+    , testSlice (eigSH . trustSym) (gen 5 :: Matrix C)
     , testSlice eigenvalues   (agen 5 :: Matrix R)
     , testSlice eigenvalues   (agen 5 :: Matrix C)
-    , testSlice eigenvaluesSH (gen 5 :: Matrix R)
-    , testSlice eigenvaluesSH (gen 5 :: Matrix C)
+    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix R)
+    , testSlice (eigenvaluesSH . trustSym) (gen 5 :: Matrix C)
 
     , testSlice svd           (rec :: Matrix R)
     , testSlice thinSVD       (rec :: Matrix R)
@@ -489,10 +489,10 @@ sliceTest = utest "slice test" $ and
     , testSlice ((<>) (ogen 5:: Matrix (Z ./. 7))) (gen 5)
     , testSlice (flip (<>) (gen 5:: Matrix (Z ./. 7))) (ogen 5)
 
-    , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ gen 5)
-    , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ gen 5)
-    , testSlice (cholSolve (chol $ gen 5:: Matrix R)) (agen 5)
-    , testSlice (cholSolve (chol $ gen 5:: Matrix C)) (agen 5)
+    , testSlice (flip cholSolve (agen 5:: Matrix R)) (chol $ trustSym $ gen 5)
+    , testSlice (flip cholSolve (agen 5:: Matrix C)) (chol $ trustSym $ gen 5)
+    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix R)) (agen 5)
+    , testSlice (cholSolve (chol $ trustSym $ gen 5:: Matrix C)) (agen 5)
 
     , ok_qrgr        (rec :: Matrix R)
     , ok_qrgr        (rec :: Matrix C)
@@ -515,8 +515,8 @@ sliceTest = utest "slice test" $ and
 
     test_lus m = testSlice f lup
       where
-        f x = luSolve (x,p) m
-        (lup,p) = luPacked m
+        f x = luSolve (LU x p) m
+        (LU lup p) = luPacked m
 
     gen :: Numeric t => Int -> Matrix t
     gen n = diagRect 1 (konst 5 n) n n
@@ -588,11 +588,11 @@ runTests n = do
     test (linearSolveProp (luSolve.luPacked) . rSqWC)
     test (linearSolveProp (luSolve.luPacked) . cSqWC)
     putStrLn "------ ldlSolve"
-    test (linearSolveProp (ldlSolve.ldlPacked) . rSymWC)
-    test (linearSolveProp (ldlSolve.ldlPacked) . cSymWC)
+    test (linearSolvePropH (ldlSolve.ldlPacked) . rSymWC)
+    test (linearSolvePropH (ldlSolve.ldlPacked) . cSymWC)
     putStrLn "------ cholSolve"
-    test (linearSolveProp (cholSolve.chol) . rPosDef)
-    test (linearSolveProp (cholSolve.chol) . cPosDef)
+    test (linearSolveProp (cholSolve.chol.trustSym) . rPosDef)
+    test (linearSolveProp (cholSolve.chol.trustSym) . cPosDef)
     putStrLn "------ luSolveLS"
     test (linearSolveProp linearSolveLS . rSqWC)
     test (linearSolveProp linearSolveLS . cSqWC)
@@ -865,8 +865,8 @@ eigBench = do
     let m = reshape 1000 (randomVector 777 Uniform (1000*1000))
         s = m + tr m
     m `seq` s `seq` putStrLn ""
-    time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH' m)
-    time "eigenvectors symmetric 1000x1000" (snd $ eigSH' m)
+    time "eigenvalues  symmetric 1000x1000" (eigenvaluesSH (trustSym m))
+    time "eigenvectors symmetric 1000x1000" (snd $ eigSH (trustSym m))
     time "eigenvalues  general   1000x1000" (eigenvalues m)
     time "eigenvectors general   1000x1000" (snd $ eig m)
 
@@ -893,12 +893,14 @@ solveBenchN n = do
     time ("svd solve " ++ show n) (linearSolveSVD a b)
     time (" ls solve " ++ show n) (linearSolveLS a b)
     time ("    solve " ++ show n) (linearSolve a b)
-    time ("cholSolve " ++ show n) (cholSolve (chol a) b)
+--    time (" LU solve " ++ show n) (luSolve (luPacked a) b)
+    time ("LDL solve " ++ show n) (ldlSolve (ldlPacked (trustSym a)) b)
+    time ("cholSolve " ++ show n) (cholSolve (chol $ trustSym a) b)
 
 solveBench = do
     solveBenchN 500
     solveBenchN 1000
-    -- solveBenchN 1500
+    solveBenchN 1500
 
 --------------------------------
 
@@ -906,7 +908,7 @@ cholBenchN n = do
     let x = uniformSample 777 (2*n) (replicate n (-1,1))
         a = tr x <> x
     a `seq` putStr ""
-    time ("chol " ++ show n) (chol a)
+    time ("chol " ++ show n) (chol $ trustSym a)
 
 cholBench = do
     putStrLn ""
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 7c54535..4704989 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -14,7 +14,7 @@ Arbitrary instances for vectors, matrices.
 module Numeric.LinearAlgebra.Tests.Instances(
     Sq(..),     rSq,cSq,
     Rot(..),    rRot,cRot,
-    Her(..),    rHer,cHer,
+                rHer,cHer,
     WC(..),     rWC,cWC,
     SqWC(..),   rSqWC, cSqWC, rSymWC, cSymWC,
     PosDef(..), rPosDef, cPosDef,
@@ -81,12 +81,12 @@ instance (Field a, Arbitrary a) => Arbitrary (Rot a) where
 
 
 -- a complex hermitian or real symmetric matrix
-newtype (Her a) = Her (Matrix a) deriving Show
+--newtype (Her a) = Her (Matrix a) deriving Show
 instance (Field a, Arbitrary a, Num (Vector a)) => Arbitrary (Her a) where
     arbitrary = do
         Sq m <- arbitrary
         let m' = m/2
-        return $ Her (m' + tr m')
+        return $ sym m'
 
 
 class (Field a, Arbitrary a, Element (RealOf a), Random (RealOf a)) => ArbitraryField a
@@ -125,9 +125,9 @@ newtype (PosDef a) = PosDef (Matrix a) deriving Show
 instance (Numeric a, ArbitraryField a, Num (Vector a)) 
     => Arbitrary (PosDef a) where
     arbitrary = do
-        Her m <- arbitrary
+        m <- arbitrary
         let (_,v) = eigSH m
-            n = rows m
+            n = rows (her m)
         l <- replicateM n (choose (0,100))
         let s = diag (fromList l)
             p = v <> real s <> tr v
@@ -161,8 +161,8 @@ fM m = m :: FM
 zM m = m :: ZM
 
 
-rHer (Her m) = m :: RM
-cHer (Her m) = m :: CM
+rHer m = her m :: RM
+cHer m = her m :: CM
 
 rRot (Rot m) = m :: RM
 cRot (Rot m) = m :: CM
@@ -176,8 +176,8 @@ cWC (WC m) = m :: CM
 rSqWC (SqWC m) = m :: RM
 cSqWC (SqWC m) = m :: CM
 
-rSymWC (SqWC m) = m + tr m :: RM
-cSymWC (SqWC m) = m + tr m :: CM
+rSymWC (SqWC m) = sym m :: Her R
+cSymWC (SqWC m) = sym m :: Her C
 
 rPosDef (PosDef m) = m :: RM
 cPosDef (PosDef m) = m :: CM
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index 207a303..2ac3588 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -39,7 +39,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
     expmDiagProp,
     multProp1, multProp2,
     subProp,
-    linearSolveProp, linearSolveProp2
+    linearSolveProp, linearSolvePropH, linearSolveProp2
 ) where
 
 import Numeric.LinearAlgebra.HMatrix hiding (Testable,unitary)
@@ -209,11 +209,11 @@ eigProp m = complex m <> v |~| v <> diag s
 eigSHProp m = m <> v |~| v <> real (diag s)
               && unitary v
               && m |~| v <> real (diag s) <> tr v
-    where (s, v) = eigSH m
+    where (s, v) = eigSH' m
 
 eigProp2 m = fst (eig m) |~| eigenvalues m
 
-eigSHProp2 m = fst (eigSH m) |~| eigenvaluesSH m
+eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m
 
 ------------------------------------------------------------------
 
@@ -246,9 +246,9 @@ schurProp2 m = m |~| u <> s <> tr u && unitary u && upperHessenberg s -- fixme
     where (u,s) = schur m
 
 cholProp m = m |~| tr c <> c && upperTriang c
-    where c = chol m
+    where c = chol (trustSym m)
 
-exactProp m = chol m == chol (m+0)
+exactProp m = chol (trustSym m) == chol (trustSym (m+0))
 
 expmDiagProp m = expm (logm m) :~ 7 ~: complex m
     where logm = matFunc log
@@ -263,6 +263,8 @@ multProp2 p (a,b) = (tr (a <> b)) :~p~: (tr b <> tr a)
 
 linearSolveProp f m = f m m |~| ident (rows m)
 
+linearSolvePropH f m = f m (her m) |~| ident (rows (her m))
+
 linearSolveProp2 f (a,x) = not wc `trivial` (not wc || a <> f a b |~| b)
     where q = min (rows a) (cols a)
           b = a <> x