Merge pull request #285 from maksbotan/add-geig

Add generalized eigenvalues via dggev and zggev

7 files changed, 201 insertions, 2 deletions

diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs
index 6027c46..f5bddc6 100644
--- a/packages/base/src/Internal/Algorithms.hs
+++ b/packages/base/src/Internal/Algorithms.hs
@@ -70,8 +70,10 @@ class (Numeric t,
     linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
     linearSolveLS'  :: Matrix t -> Matrix t -> Matrix t
     eig'         :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
+    geig'        :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double))
     eigSH''      :: Matrix t -> (Vector Double, Matrix t)
     eigOnly      :: Matrix t -> Vector (Complex Double)
+    geigOnly     :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)
     eigOnlySH    :: Matrix t -> Vector Double
     cholSH'      :: Matrix t -> Matrix t
     mbCholSH'    :: Matrix t -> Maybe (Matrix t)
@@ -96,7 +98,9 @@ instance Field Double where
     linearSolveSVD' = linearSolveSVDR Nothing
     eig' = eigR
     eigSH'' = eigS
+    geig' = eigG
     eigOnly = eigOnlyR
+    geigOnly = eigOnlyG
     eigOnlySH = eigOnlyS
     cholSH' = cholS
     mbCholSH' = mbCholS
@@ -126,7 +130,9 @@ instance Field (Complex Double) where
     linearSolveLS' = linearSolveLSC
     linearSolveSVD' = linearSolveSVDC Nothing
     eig' = eigC
+    geig' = eigGC
     eigOnly = eigOnlyC
+    geigOnly = eigOnlyGC
     eigSH'' = eigH
     eigOnlySH = eigOnlyH
     cholSH' = cholH
@@ -512,10 +518,25 @@ a = (3><3)
 eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
 eig = {-# SCC "eig" #-} eig'
 
+-- | Generalized eigenvalues (not ordered) and eigenvectors (as columns) of a pair of nonsymmetric matrices.
+-- Eigenvalues are represented as pairs of alpha, beta, where eigenvalue = alpha / beta. Alpha is always
+-- complex, but betas has the same type as the input matrix.
+--
+-- If @(alphas, betas, v) = geig a b@, then @a \<> v == b \<> v \<> diag (alphas / betas)@
+--
+-- Note that beta can be 0 and that has reasonable interpretation.
+geig :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double))
+geig = {-# SCC "geig" #-} geig'
+
 -- | Eigenvalues (not ordered) of a general square matrix.
 eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
 eigenvalues = {-# SCC "eigenvalues" #-} eigOnly
 
+-- | Generalized eigenvalues of a pair of matrices. Represented as pairs of alpha, beta,
+-- where eigenvalue is alpha / beta as in 'geig'.
+geigenvalues :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)
+geigenvalues = {-# SCC "geigenvalues" #-} geigOnly
+
 -- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
 eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
 eigSH' = {-# SCC "eigSH'" #-} eigSH''
diff --git a/packages/base/src/Internal/C/lapack-aux.c b/packages/base/src/Internal/C/lapack-aux.c
index 5018a45..dddbefb 100644
--- a/packages/base/src/Internal/C/lapack-aux.c
+++ b/packages/base/src/Internal/C/lapack-aux.c
@@ -415,6 +415,93 @@ int eig_l_R(ODMAT(a),ODMAT(u), CVEC(s),ODMAT(v)) {
     OK
 }
 
+//////////////////// generalized real eigensystem ////////////
+
+int dggev_(char *jobvl, char *jobvr, integer *n,
+    doublereal *a, integer *lda, doublereal *b, integer *ldb,
+    doublereal *alphar, doublereal *alphai, doublereal *beta,
+    doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr,
+    doublereal *work,
+        integer *lwork, integer *info);
+
+int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) {
+    integer n = ar;
+    REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);
+    REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);
+    char jobvl = vlp==NULL?'N':'V';
+    REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);
+    char jobvr = vrp==NULL?'N':'V';
+    DEBUGMSG("eig_l_G");
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    double ans;
+    dggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             (double*)alphap, (double*)alphap+n, betap,
+             vlp, &n, vrp, &n,
+             &ans, &lwork,
+             &res);
+    lwork = ceil(ans);
+    double * work = (double*)malloc(lwork*sizeof(double));
+    CHECK(!work,MEM);
+    dggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             (double*)alphap, (double*)alphap+n, betap,
+             vlp, &n, vrp, &n,
+             work, &lwork,
+             &res);
+    CHECK(res,res);
+    free(work);
+    OK
+}
+
+//////////////////// generalized complex eigensystem ////////////
+
+int zggev_(char *jobvl, char *jobvr, integer *n,
+    doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
+    doublecomplex *alphar, doublecomplex *beta,
+    doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr,
+    doublecomplex *work, integer *lwork,
+    doublereal *rwork, integer *info);
+
+int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) {
+    integer n = ar;
+    REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);
+    REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);
+    char jobvl = vlp==NULL?'N':'V';
+    REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);
+    char jobvr = vrp==NULL?'N':'V';
+    DEBUGMSG("eig_l_GC");
+    double *rwork = (double*) malloc(8*n*sizeof(double));
+    CHECK(!rwork,MEM);
+    integer lwork = -1;
+    integer res;
+    // ask for optimal lwork
+    doublecomplex ans;
+    zggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             alphap, betap,
+             vlp, &n, vrp, &n,
+             &ans, &lwork,
+             rwork, &res);
+    lwork = ceil(ans.r);
+    doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
+    CHECK(!work,MEM);
+    zggev_  (&jobvl,&jobvr,
+             &n,
+             ap,&n,bp,&n,
+             alphap, betap,
+             vlp, &n, vrp, &n,
+             work, &lwork,
+             rwork, &res);
+    CHECK(res,res);
+    free(work);
+    OK
+}
 
 //////////////////// symmetric real eigensystem ////////////
 
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs
index 64cf2f5..ff55688 100644
--- a/packages/base/src/Internal/LAPACK.hs
+++ b/packages/base/src/Internal/LAPACK.hs
@@ -18,6 +18,8 @@
 
 module Internal.LAPACK where
 
+import Data.Bifunctor (first)
+
 import Internal.Devel
 import Internal.Vector
 import Internal.Matrix hiding ((#), (#!))
@@ -234,7 +236,9 @@ leftSVAux f st x = unsafePerformIO $ do
 -----------------------------------------------------------------------------
 
 foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok
+foreign import ccall unsafe "eig_l_G" dggev :: R ::> R ::> C :> R :> R ::> R ::> Ok
 foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok
+foreign import ccall unsafe "eig_l_GC" zggev :: C ::> C ::> C :> C :> C ::> C ::> Ok
 foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok
 foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok
 
@@ -298,12 +302,68 @@ 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 ((_:+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)
+fixeigG _ _ = error "fixeigG with impossible inputs"
 
 -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
 -- The eigenvalues are not sorted.
 eigOnlyR :: Matrix Double -> Vector (Complex Double)
 eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR"
 
+-- | Generalized eigenvalues and right eigenvectors of a pair of real matrices, using LAPACK's /dggev/.
+-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
+eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double))
+eigG a b = (alpha', beta, v'')
+  where
+    (alpha, beta, v) = eigGaux dggev a b "eigG"
+    alpha' = fixeig1 alpha
+    v' = toRows $ trans v
+    v'' = fromColumns $ fixeigG (toList alpha') v'
+
+eigGaux f ma mb st = unsafePerformIO $ do
+    a <- copy ColumnMajor ma
+    b <- copy ColumnMajor mb
+    alpha <- createVector r
+    beta <- createVector r
+    vr <- createMatrix ColumnMajor r r
+
+    (a # b # alpha # beta #! vr) g #| st
+
+    return (alpha, beta, vr)
+  where
+    r = rows ma
+    g ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta = f ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta 0 0 0 0 nullPtr 
+
+eigGOnlyAux f ma mb st = unsafePerformIO $ do
+    a <- copy ColumnMajor ma
+    b <- copy ColumnMajor mb
+    alpha <- createVector r
+    beta <- createVector r
+
+    (a # b # alpha #! beta) g #| st
+
+    return (alpha, beta)
+  where
+    r = rows ma
+    g ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta = f ar ac xra xca pa br bc xrb xcb pb alphan palpha betan pbeta 0 0 0 0 nullPtr 0 0 0 0 nullPtr
+
+-- | Generalized eigenvalues and right eigenvectors of a pair of complex matrices, using LAPACK's /zggev/.
+-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
+eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double))
+eigGC a b = eigGaux zggev a b "eigGC"
+
+eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double)
+eigOnlyG a b = first fixeig1 $ eigGOnlyAux dggev a b "eigOnlyG"
+
+eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double))
+eigOnlyGC a b = eigGOnlyAux zggev a b "eigOnlyGC"
 
 -----------------------------------------------------------------------------
 
diff --git a/packages/base/src/Numeric/LinearAlgebra.hs b/packages/base/src/Numeric/LinearAlgebra.hs
index 91923e9..9670187 100644
--- a/packages/base/src/Numeric/LinearAlgebra.hs
+++ b/packages/base/src/Numeric/LinearAlgebra.hs
@@ -131,8 +131,8 @@ module Numeric.LinearAlgebra (
     leftSV, rightSV,
 
     -- * Eigendecomposition
-    eig, eigSH,
-    eigenvalues, eigenvaluesSH,
+    eig, geig, eigSH,
+    eigenvalues, geigenvalues, eigenvaluesSH,
     geigSH,
 
     -- * QR
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