generalized symmetric eigensystem

1 files changed, 18 insertions, 0 deletions

diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs
index a6b3174..bea33ea 100644
--- a/lib/Numeric/LinearAlgebra/Algorithms.hs
+++ b/lib/Numeric/LinearAlgebra/Algorithms.hs
@@ -44,6 +44,7 @@ module Numeric.LinearAlgebra.Algorithms (
 -- ** Eigensystems
     eig, eigSH, eigSH',
     eigenvalues, eigenvaluesSH, eigenvaluesSH',
+    geigSH',
 -- ** QR
     qr, rq,
 -- ** Cholesky
@@ -712,3 +713,20 @@ relativeError x y = dig (norm (x `sub` y) / norm x)
     where norm = pnorm Infinity
           dig r = round $ -logBase 10 (realToFrac r :: Double)
 
+----------------------------------------------------------------------
+
+-- | Generalized symmetric positive definite eigensystem Av = lBv,
+-- for A and B symmetric, B positive definite (conditions not checked).
+geigSH' :: Field t
+        => Matrix t -- ^ A
+        -> Matrix t -- ^ B
+        -> (Vector Double, Matrix t)
+geigSH' a b = (l,v')
+  where
+    u = cholSH b
+    iu = inv u
+    c = ctrans iu <> a <> iu
+    (l,v) = eigSH' c
+    v' = iu <> v
+    (<>) = mXm
+