summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
authoridontgetoutmuch <dominic@steinitz.org>2018-11-17 09:31:17 +0000
committerGitHub <noreply@github.com>2018-11-17 09:31:17 +0000
commit480c4db1585ed122dfc491bb15421565966dad66 (patch)
treef2270a62ef1d07a47e18a8c7c355714284f209c2 /packages
parent94c7f89929e09ca6e98976cb59dec574ee164e20 (diff)
parent9b37d60cf68adf1201163e80ab67f6c7706b5d76 (diff)
Merge pull request #285 from maksbotan/add-geig
Add generalized eigenvalues via dggev and zggev
Diffstat (limited to 'packages')
-rw-r--r--packages/base/src/Internal/Algorithms.hs21
-rw-r--r--packages/base/src/Internal/C/lapack-aux.c87
-rw-r--r--packages/base/src/Internal/LAPACK.hs60
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs4
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests.hs3
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Instances.hs21
-rw-r--r--packages/tests/src/Numeric/LinearAlgebra/Tests/Properties.hs7
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,
70 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t 70 linearSolveSVD' :: Matrix t -> Matrix t -> Matrix t
71 linearSolveLS' :: Matrix t -> Matrix t -> Matrix t 71 linearSolveLS' :: Matrix t -> Matrix t -> Matrix t
72 eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 72 eig' :: Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
73 geig' :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double))
73 eigSH'' :: Matrix t -> (Vector Double, Matrix t) 74 eigSH'' :: Matrix t -> (Vector Double, Matrix t)
74 eigOnly :: Matrix t -> Vector (Complex Double) 75 eigOnly :: Matrix t -> Vector (Complex Double)
76 geigOnly :: Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)
75 eigOnlySH :: Matrix t -> Vector Double 77 eigOnlySH :: Matrix t -> Vector Double
76 cholSH' :: Matrix t -> Matrix t 78 cholSH' :: Matrix t -> Matrix t
77 mbCholSH' :: Matrix t -> Maybe (Matrix t) 79 mbCholSH' :: Matrix t -> Maybe (Matrix t)
@@ -96,7 +98,9 @@ instance Field Double where
96 linearSolveSVD' = linearSolveSVDR Nothing 98 linearSolveSVD' = linearSolveSVDR Nothing
97 eig' = eigR 99 eig' = eigR
98 eigSH'' = eigS 100 eigSH'' = eigS
101 geig' = eigG
99 eigOnly = eigOnlyR 102 eigOnly = eigOnlyR
103 geigOnly = eigOnlyG
100 eigOnlySH = eigOnlyS 104 eigOnlySH = eigOnlyS
101 cholSH' = cholS 105 cholSH' = cholS
102 mbCholSH' = mbCholS 106 mbCholSH' = mbCholS
@@ -126,7 +130,9 @@ instance Field (Complex Double) where
126 linearSolveLS' = linearSolveLSC 130 linearSolveLS' = linearSolveLSC
127 linearSolveSVD' = linearSolveSVDC Nothing 131 linearSolveSVD' = linearSolveSVDC Nothing
128 eig' = eigC 132 eig' = eigC
133 geig' = eigGC
129 eigOnly = eigOnlyC 134 eigOnly = eigOnlyC
135 geigOnly = eigOnlyGC
130 eigSH'' = eigH 136 eigSH'' = eigH
131 eigOnlySH = eigOnlyH 137 eigOnlySH = eigOnlyH
132 cholSH' = cholH 138 cholSH' = cholH
@@ -512,10 +518,25 @@ a = (3><3)
512eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double)) 518eig :: Field t => Matrix t -> (Vector (Complex Double), Matrix (Complex Double))
513eig = {-# SCC "eig" #-} eig' 519eig = {-# SCC "eig" #-} eig'
514 520
521-- | Generalized eigenvalues (not ordered) and eigenvectors (as columns) of a pair of nonsymmetric matrices.
522-- Eigenvalues are represented as pairs of alpha, beta, where eigenvalue = alpha / beta. Alpha is always
523-- complex, but betas has the same type as the input matrix.
524--
525-- If @(alphas, betas, v) = geig a b@, then @a \<> v == b \<> v \<> diag (alphas / betas)@
526--
527-- Note that beta can be 0 and that has reasonable interpretation.
528geig :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t, Matrix (Complex Double))
529geig = {-# SCC "geig" #-} geig'
530
515-- | Eigenvalues (not ordered) of a general square matrix. 531-- | Eigenvalues (not ordered) of a general square matrix.
516eigenvalues :: Field t => Matrix t -> Vector (Complex Double) 532eigenvalues :: Field t => Matrix t -> Vector (Complex Double)
517eigenvalues = {-# SCC "eigenvalues" #-} eigOnly 533eigenvalues = {-# SCC "eigenvalues" #-} eigOnly
518 534
535-- | Generalized eigenvalues of a pair of matrices. Represented as pairs of alpha, beta,
536-- where eigenvalue is alpha / beta as in 'geig'.
537geigenvalues :: Field t => Matrix t -> Matrix t -> (Vector (Complex Double), Vector t)
538geigenvalues = {-# SCC "geigenvalues" #-} geigOnly
539
519-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part. 540-- | Similar to 'eigSH' without checking that the input matrix is hermitian or symmetric. It works with the upper triangular part.
520eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t) 541eigSH' :: Field t => Matrix t -> (Vector Double, Matrix t)
521eigSH' = {-# SCC "eigSH'" #-} eigSH'' 542eigSH' = {-# 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)) {
415 OK 415 OK
416} 416}
417 417
418//////////////////// generalized real eigensystem ////////////
419
420int dggev_(char *jobvl, char *jobvr, integer *n,
421 doublereal *a, integer *lda, doublereal *b, integer *ldb,
422 doublereal *alphar, doublereal *alphai, doublereal *beta,
423 doublereal *vl, integer *ldvl, doublereal *vr, integer *ldvr,
424 doublereal *work,
425 integer *lwork, integer *info);
426
427int eig_l_G(ODMAT(a), ODMAT(b), CVEC(alpha), DVEC(beta), ODMAT(vl), ODMAT(vr)) {
428 integer n = ar;
429 REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);
430 REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);
431 char jobvl = vlp==NULL?'N':'V';
432 REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);
433 char jobvr = vrp==NULL?'N':'V';
434 DEBUGMSG("eig_l_G");
435 integer lwork = -1;
436 integer res;
437 // ask for optimal lwork
438 double ans;
439 dggev_ (&jobvl,&jobvr,
440 &n,
441 ap,&n,bp,&n,
442 (double*)alphap, (double*)alphap+n, betap,
443 vlp, &n, vrp, &n,
444 &ans, &lwork,
445 &res);
446 lwork = ceil(ans);
447 double * work = (double*)malloc(lwork*sizeof(double));
448 CHECK(!work,MEM);
449 dggev_ (&jobvl,&jobvr,
450 &n,
451 ap,&n,bp,&n,
452 (double*)alphap, (double*)alphap+n, betap,
453 vlp, &n, vrp, &n,
454 work, &lwork,
455 &res);
456 CHECK(res,res);
457 free(work);
458 OK
459}
460
461//////////////////// generalized complex eigensystem ////////////
462
463int zggev_(char *jobvl, char *jobvr, integer *n,
464 doublecomplex *a, integer *lda, doublecomplex *b, integer *ldb,
465 doublecomplex *alphar, doublecomplex *beta,
466 doublecomplex *vl, integer *ldvl, doublecomplex *vr, integer *ldvr,
467 doublecomplex *work, integer *lwork,
468 doublereal *rwork, integer *info);
469
470int eig_l_GC(OCMAT(a), OCMAT(b), CVEC(alpha), CVEC(beta), OCMAT(vl), OCMAT(vr)) {
471 integer n = ar;
472 REQUIRES(ac == n && br == n && bc == n && alphan == n && betan == n, BAD_SIZE);
473 REQUIRES(vlp==NULL || (vlr==n && vlc==n), BAD_SIZE);
474 char jobvl = vlp==NULL?'N':'V';
475 REQUIRES(vrp==NULL || (vrr==n && vrc==n), BAD_SIZE);
476 char jobvr = vrp==NULL?'N':'V';
477 DEBUGMSG("eig_l_GC");
478 double *rwork = (double*) malloc(8*n*sizeof(double));
479 CHECK(!rwork,MEM);
480 integer lwork = -1;
481 integer res;
482 // ask for optimal lwork
483 doublecomplex ans;
484 zggev_ (&jobvl,&jobvr,
485 &n,
486 ap,&n,bp,&n,
487 alphap, betap,
488 vlp, &n, vrp, &n,
489 &ans, &lwork,
490 rwork, &res);
491 lwork = ceil(ans.r);
492 doublecomplex * work = (doublecomplex*)malloc(lwork*sizeof(doublecomplex));
493 CHECK(!work,MEM);
494 zggev_ (&jobvl,&jobvr,
495 &n,
496 ap,&n,bp,&n,
497 alphap, betap,
498 vlp, &n, vrp, &n,
499 work, &lwork,
500 rwork, &res);
501 CHECK(res,res);
502 free(work);
503 OK
504}
418 505
419//////////////////// symmetric real eigensystem //////////// 506//////////////////// symmetric real eigensystem ////////////
420 507
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 @@
18 18
19module Internal.LAPACK where 19module Internal.LAPACK where
20 20
21import Data.Bifunctor (first)
22
21import Internal.Devel 23import Internal.Devel
22import Internal.Vector 24import Internal.Vector
23import Internal.Matrix hiding ((#), (#!)) 25import Internal.Matrix hiding ((#), (#!))
@@ -234,7 +236,9 @@ leftSVAux f st x = unsafePerformIO $ do
234----------------------------------------------------------------------------- 236-----------------------------------------------------------------------------
235 237
236foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok 238foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok
239foreign import ccall unsafe "eig_l_G" dggev :: R ::> R ::> C :> R :> R ::> R ::> Ok
237foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok 240foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok
241foreign import ccall unsafe "eig_l_GC" zggev :: C ::> C ::> C :> C :> C ::> C ::> Ok
238foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok 242foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok
239foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok 243foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok
240 244
@@ -298,12 +302,68 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs)
298 | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) 302 | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs)
299fixeig _ _ = error "fixeig with impossible inputs" 303fixeig _ _ = error "fixeig with impossible inputs"
300 304
305-- For dggev alpha(i) / beta(i), alpha(i+1) / beta(i+1) form a complex conjugate pair when Im alpha(i) != 0.
306-- However, this does not lead to Re alpha(i) == Re alpha(i+1), since beta(i) and beta(i+1)
307-- can be different. Therefore old 'fixeig' would fail for 'eigG'.
308fixeigG [] _ = []
309fixeigG [_] [v] = [comp' v]
310fixeigG ((_:+ai1) : an : as) (v1:v2:vs)
311 | abs ai1 > 1e-13 = toComplex' (v1, v2) : toComplex' (v1, mapVector negate v2) : fixeigG as vs
312 | otherwise = comp' v1 : fixeigG (an:as) (v2:vs)
313fixeigG _ _ = error "fixeigG with impossible inputs"
301 314
302-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. 315-- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'.
303-- The eigenvalues are not sorted. 316-- The eigenvalues are not sorted.
304eigOnlyR :: Matrix Double -> Vector (Complex Double) 317eigOnlyR :: Matrix Double -> Vector (Complex Double)
305eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" 318eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR"
306 319
320-- | Generalized eigenvalues and right eigenvectors of a pair of real matrices, using LAPACK's /dggev/.
321-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
322eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double))
323eigG a b = (alpha', beta, v'')
324 where
325 (alpha, beta, v) = eigGaux dggev a b "eigG"
326 alpha' = fixeig1 alpha
327 v' = toRows $ trans v
328 v'' = fromColumns $ fixeigG (toList alpha') v'
329
330eigGaux f ma mb st = unsafePerformIO $ do
331 a <- copy ColumnMajor ma
332 b <- copy ColumnMajor mb
333 alpha <- createVector r
334 beta <- createVector r
335 vr <- createMatrix ColumnMajor r r
336
337 (a # b # alpha # beta #! vr) g #| st
338
339 return (alpha, beta, vr)
340 where
341 r = rows ma
342 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
343
344eigGOnlyAux f ma mb st = unsafePerformIO $ do
345 a <- copy ColumnMajor ma
346 b <- copy ColumnMajor mb
347 alpha <- createVector r
348 beta <- createVector r
349
350 (a # b # alpha #! beta) g #| st
351
352 return (alpha, beta)
353 where
354 r = rows ma
355 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
356
357-- | Generalized eigenvalues and right eigenvectors of a pair of complex matrices, using LAPACK's /zggev/.
358-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
359eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double))
360eigGC a b = eigGaux zggev a b "eigGC"
361
362eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double)
363eigOnlyG a b = first fixeig1 $ eigGOnlyAux dggev a b "eigOnlyG"
364
365eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double))
366eigOnlyGC a b = eigGOnlyAux zggev a b "eigOnlyGC"
307 367
308----------------------------------------------------------------------------- 368-----------------------------------------------------------------------------
309 369
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 (
131 leftSV, rightSV, 131 leftSV, rightSV,
132 132
133 -- * Eigendecomposition 133 -- * Eigendecomposition
134 eig, eigSH, 134 eig, geig, eigSH,
135 eigenvalues, eigenvaluesSH, 135 eigenvalues, geigenvalues, eigenvaluesSH,
136 geigSH, 136 geigSH,
137 137
138 -- * QR 138 -- * 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
756 test (eigSHProp2 . cHer) 756 test (eigSHProp2 . cHer)
757 test (eigProp2 . rSq) 757 test (eigProp2 . rSq)
758 test (eigProp2 . cSq) 758 test (eigProp2 . cSq)
759 putStrLn "------ geig"
760 test (uncurry geigProp . rSq2WC)
761 test (uncurry geigProp . cSq2WC)
759 putStrLn "------ nullSpace" 762 putStrLn "------ nullSpace"
760 test (nullspaceProp . rM) 763 test (nullspaceProp . rM)
761 test (nullspaceProp . cM) 764 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.
17 17
18module Numeric.LinearAlgebra.Tests.Instances( 18module Numeric.LinearAlgebra.Tests.Instances(
19 Sq(..), rSq,cSq, 19 Sq(..), rSq,cSq,
20 Sq2WC(..), rSq2WC,cSq2WC,
20 Rot(..), rRot,cRot, 21 Rot(..), rRot,cRot,
21 rHer,cHer, 22 rHer,cHer,
22 WC(..), rWC,cWC, 23 WC(..), rWC,cWC,
@@ -105,6 +106,23 @@ instance (Element a, Arbitrary a) => Arbitrary (Sq a) where
105 106
106 shrink (Sq a) = [ Sq b | b <- shrink a ] 107 shrink (Sq a) = [ Sq b | b <- shrink a ]
107 108
109-- a pair of square matrices
110newtype (Sq2WC a) = Sq2WC (Matrix a, Matrix a) deriving Show
111instance (ArbitraryField a, Numeric a) => Arbitrary (Sq2WC a) where
112 arbitrary = do
113 n <- chooseDim
114 l <- vector (n*n)
115 r <- vector (n*n)
116 l' <- makeWC $ (n><n) l
117 r' <- makeWC $ (n><n) r
118 return $ Sq2WC (l', r')
119 where
120 makeWC m = do
121 let (u,_,v) = svd m
122 n = rows m
123 sv' <- replicateM n (choose (1,100))
124 let s = diag (fromList sv')
125 return $ u <> real s <> tr v
108 126
109-- a unitary matrix 127-- a unitary matrix
110newtype (Rot a) = Rot (Matrix a) deriving Show 128newtype (Rot a) = Rot (Matrix a) deriving Show
@@ -204,6 +222,9 @@ cRot (Rot m) = m :: CM
204rSq (Sq m) = m :: RM 222rSq (Sq m) = m :: RM
205cSq (Sq m) = m :: CM 223cSq (Sq m) = m :: CM
206 224
225rSq2WC (Sq2WC (a, b)) = (a, b) :: (RM, RM)
226cSq2WC (Sq2WC (a, b)) = (a, b) :: (CM, CM)
227
207rWC (WC m) = m :: RM 228rWC (WC m) = m :: RM
208cWC (WC m) = m :: CM 229cWC (WC m) = m :: CM
209 230
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 (
36 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4, 36 svdProp1, svdProp1a, svdProp1b, svdProp2, svdProp3, svdProp4,
37 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7, 37 svdProp5a, svdProp5b, svdProp6a, svdProp6b, svdProp7,
38 eigProp, eigSHProp, eigProp2, eigSHProp2, 38 eigProp, eigSHProp, eigProp2, eigSHProp2,
39 geigProp,
39 qrProp, rqProp, rqProp1, rqProp2, rqProp3, 40 qrProp, rqProp, rqProp1, rqProp2, rqProp3,
40 hessProp, 41 hessProp,
41 schurProp1, schurProp2, 42 schurProp1, schurProp2,
@@ -237,6 +238,12 @@ eigProp2 m = fst (eig m) |~| eigenvalues m
237 238
238eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m 239eigSHProp2 m = fst (eigSH' m) |~| eigenvaluesSH' m
239 240
241geigProp a b = a' <> v <> diag betas' |~| b' <> v <> diag alphas
242 where (alphas, betas, v) = geig a b
243 betas' = complex betas
244 a' = complex a
245 b' = complex b
246
240------------------------------------------------------------------ 247------------------------------------------------------------------
241 248
242qrProp m = q <> r |~| m && unitary q && upperTriang r 249qrProp m = q <> r |~| m && unitary q && upperTriang r