From cc737710b4878a387ea8090f9c2330b1e8d73f38 Mon Sep 17 00:00:00 2001
From: Maxim Koltsov <kolmax94@gmail.com>
Date: Tue, 13 Nov 2018 15:39:55 +0300
Subject: 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.
---
 packages/base/src/Internal/LAPACK.hs                | 10 +++++++++-
 packages/tests/src/Numeric/LinearAlgebra/Tests.hs   |  3 +++
 .../src/Numeric/LinearAlgebra/Tests/Instances.hs    | 21 +++++++++++++++++++++
 .../src/Numeric/LinearAlgebra/Tests/Properties.hs   |  7 +++++++
 4 files changed, 40 insertions(+), 1 deletion(-)

diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index b834b20..c9f5d0f 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -302,6 +302,14 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
     | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs)
 fixeig _ _ = error "fixeig with impossible inputs"
 
+-- For dggev alpha(i) / beta(i), alpha(i+1) / beta(i+1) form a complex conjugate pair when Im alpha(i) != 0.
+-- However, this does not lead to Re alpha(i) == Re alpha(i+1), since beta(i) and beta(i+1)
+-- can be different. Therefore old 'fixeig' would fail for 'eigG'.
+fixeigG  []  _  = []
+fixeigG [_] [v] = [comp' v]
+fixeigG ((ar1:+ai1) : an : as) (v1:v2:vs)
+    | abs ai1 > 1e-13 = toComplex' (v1, v2) : toComplex' (v1, mapVector negate v2) : fixeigG as vs
+    | otherwise = comp' v1 : fixeigG (an:as) (v2:vs)
 
 -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
 -- The eigenvalues are not sorted.
@@ -316,7 +324,7 @@ eigG a b = (alpha', beta, v'')
     (alpha, beta, v) = eigGaux dggev a b "eigG"
     alpha' = fixeig1 alpha
     v' = toRows $ trans v
-    v'' = fromColumns $ fixeig (toList alpha') v'
+    v'' = fromColumns $ fixeigG (toList alpha') v'
 
 eigGaux f ma mb st = unsafePerformIO $ do
     a <- copy ColumnMajor ma
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
-- 
cgit v1.2.3