Merge pull request #285 from maksbotan/add-geig

Add generalized eigenvalues via dggev and zggev

4 files changed, 170 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