Fix geig for complex eigenvalues, add tests

eigG was incorrectly returning eigenvectors corresponding to complex
eigenvalues. This was discovered with tests for geig, which are also
added to the repo.

3 files changed, 31 insertions, 0 deletions

diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
index 4ed1462..c0c151a 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests.hs
@@ -756,6 +756,9 @@ runTests n = do
     test (eigSHProp2 . cHer)
     test (eigProp2   . rSq)
     test (eigProp2   . cSq)
+    putStrLn "------ geig"
+    test (uncurry geigProp . rSq2WC)
+    test (uncurry geigProp . cSq2WC)
     putStrLn "------ nullSpace"
     test (nullspaceProp . rM)
     test (nullspaceProp . cM)
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
index 59230e0..97cfd01 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs
@@ -17,6 +17,7 @@ Arbitrary instances for vectors, matrices.
 
 module Numeric.LinearAlgebra.Tests.Instances(
     Sq(..),     rSq,cSq,
+    Sq2WC(..),  rSq2WC,cSq2WC,
     Rot(..),    rRot,cRot,
                 rHer,cHer,
     WC(..),     rWC,cWC,
@@ -105,6 +106,23 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
 
     shrink (Sq a) = [ Sq b | b <- shrink a ]
 
+-- a pair of square matrices
+newtype (Sq2WC a) = Sq2WC (Matrix a, Matrix a) deriving Show
+instance (ArbitraryField a, Numeric a) => Arbitrary (Sq2WC a) where
+    arbitrary = do
+        n <- chooseDim
+        l <- vector (n*n)
+        r <- vector (n*n)
+        l' <- makeWC $ (n><n) l
+        r' <- makeWC $ (n><n) r
+        return $ Sq2WC (l', r')
+        where
+            makeWC m = do
+              let (u,_,v) = svd m
+                  n = rows m
+              sv' <- replicateM n (choose (1,100))
+              let s = diag (fromList sv')
+              return $ u <> real s <> tr v
 
 -- a unitary matrix
 newtype (Rot a) = Rot (Matrix a) deriving Show
@@ -204,6 +222,9 @@ cRot (Rot m) = m :: CM
 rSq  (Sq m)  = m :: RM
 cSq  (Sq m)  = m :: CM
 
+rSq2WC (Sq2WC (a, b)) = (a, b) :: (RM, RM)
+cSq2WC (Sq2WC (a, b)) = (a, b) :: (CM, CM)
+
 rWC (WC m) = m :: RM
 cWC (WC m) = m :: CM
 
diff --git a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
index 6cd3a9e..38aa977 100644
--- a/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
+++ b/packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs
@@ -36,6 +36,7 @@ module Numeric.LinearAlgebra.Tests.Properties (
     svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
     svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
     eigProp, eigSHProp, eigProp2, eigSHProp2,
+    geigProp,
     qrProp, rqProp, rqProp1, rqProp2, rqProp3,
     hessProp,
     schurProp1, schurProp2,
@@ -237,6 +238,12 @@ eigProp2 m = fst (eig m) |~| eigenvalues m
 
 eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m
 
+geigProp a b = a' <> v <> diag betas' |~| b' <> v <> diag alphas
+    where (alphas, betas, v) = geig a b
+          betas' = complex betas
+          a' = complex a
+          b' = complex b
+
 ------------------------------------------------------------------
 
 qrProp m = q <> r |~| m && unitary q && upperTriang r