summaryrefslogtreecommitdiff
path: root/packages
diff options
context:
space:
mode:
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.hs51
-rw-r--r--packages/base/src/Numeric/LinearAlgebra.hs4
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,
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 == diag (alphas / betas) \<> b \<> v@
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..b834b20 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
@@ -304,6 +308,53 @@ fixeig _ _ = error "fixeig with impossible inputs"
304eigOnlyR :: Matrix Double -> Vector (Complex Double) 308eigOnlyR :: Matrix Double -> Vector (Complex Double)
305eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" 309eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR"
306 310
311-- | Generalized eigenvalues and right eigenvectors of a pair of real matrices, using LAPACK's /dggev/.
312-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
313eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double))
314eigG a b = (alpha', beta, v'')
315 where
316 (alpha, beta, v) = eigGaux dggev a b "eigG"
317 alpha' = fixeig1 alpha
318 v' = toRows $ trans v
319 v'' = fromColumns $ fixeig (toList alpha') v'
320
321eigGaux f ma mb st = unsafePerformIO $ do
322 a <- copy ColumnMajor ma
323 b <- copy ColumnMajor mb
324 alpha <- createVector r
325 beta <- createVector r
326 vr <- createMatrix ColumnMajor r r
327
328 (a # b # alpha # beta #! vr) g #| st
329
330 return (alpha, beta, vr)
331 where
332 r = rows ma
333 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
334
335eigGOnlyAux f ma mb st = unsafePerformIO $ do
336 a <- copy ColumnMajor ma
337 b <- copy ColumnMajor mb
338 alpha <- createVector r
339 beta <- createVector r
340
341 (a # b # alpha #! beta) g #| st
342
343 return (alpha, beta)
344 where
345 r = rows ma
346 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
347
348-- | Generalized eigenvalues and right eigenvectors of a pair of complex matrices, using LAPACK's /zggev/.
349-- The eigenvectors are the columns of v. The eigenvalues are represented as alphas / betas and not sorted.
350eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double))
351eigGC a b = eigGaux zggev a b "eigGC"
352
353eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double)
354eigOnlyG a b = first fixeig1 $ eigGOnlyAux dggev a b "eigOnlyG"
355
356eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double))
357eigOnlyGC a b = eigGOnlyAux zggev a b "eigOnlyGC"
307 358
308----------------------------------------------------------------------------- 359-----------------------------------------------------------------------------
309 360
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