From acb57ec8dd3011534813bdbee84563f1a830f499 Mon Sep 17 00:00:00 2001 From: Maxim Koltsov Date: Thu, 8 Nov 2018 11:07:33 +0300 Subject: Add generalized eigenvalues via dggev and zggev These lapack functions generalize dgeev and zgeev. Interface for them was added similarly to eig* functions already present in hmatrix. --- packages/base/src/Internal/Algorithms.hs | 21 ++++++++ packages/base/src/Internal/C/lapack-aux.c | 87 ++++++++++++++++++++++++++++++ packages/base/src/Internal/LAPACK.hs | 51 ++++++++++++++++++ packages/base/src/Numeric/LinearAlgebra.hs | 4 +- 4 files changed, 161 insertions(+), 2 deletions(-) diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs index 6027c46..1ca373f 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 == diag (alphas / betas) \<> b \<> v@ +-- +-- 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..b834b20 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 @@ -304,6 +308,53 @@ fixeig _ _ = error "fixeig with impossible inputs" 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 $ fixeig (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 -- cgit v1.2.3