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 From cc737710b4878a387ea8090f9c2330b1e8d73f38 Mon Sep 17 00:00:00 2001 From: Maxim Koltsov 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> 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 From 9b37d60cf68adf1201163e80ab67f6c7706b5d76 Mon Sep 17 00:00:00 2001 From: Maxim Koltsov Date: Tue, 13 Nov 2018 15:44:42 +0300 Subject: Fix doc for geig, fix warning --- packages/base/src/Internal/Algorithms.hs | 2 +- packages/base/src/Internal/LAPACK.hs | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/packages/base/src/Internal/Algorithms.hs b/packages/base/src/Internal/Algorithms.hs index 1ca373f..f5bddc6 100644 --- a/packages/base/src/Internal/Algorithms.hs +++ b/packages/base/src/Internal/Algorithms.hs @@ -522,7 +522,7 @@ eig = {-# SCC "eig" #-} eig' -- 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@ +-- 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)) diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs index c9f5d0f..ff55688 100644 --- a/packages/base/src/Internal/LAPACK.hs +++ b/packages/base/src/Internal/LAPACK.hs @@ -307,9 +307,10 @@ fixeig _ _ = error "fixeig with impossible inputs" -- can be different. Therefore old 'fixeig' would fail for 'eigG'. fixeigG [] _ = [] fixeigG [_] [v] = [comp' v] -fixeigG ((ar1:+ai1) : an : as) (v1:v2:vs) +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. -- cgit v1.2.3