summaryrefslogtreecommitdiff
path: root/packages/base/src/Internal/LAPACK.hs
diff options
context:
space:
mode:
authorMaxim Koltsov <kolmax94@gmail.com>2018-11-08 11:07:33 +0300
committerMaxim Koltsov <kolmax94@gmail.com>2018-11-08 11:07:33 +0300
commitacb57ec8dd3011534813bdbee84563f1a830f499 (patch)
treed2378dac72cd2d6d8c0b73375c12929a9c97daf2 /packages/base/src/Internal/LAPACK.hs
parent94c7f89929e09ca6e98976cb59dec574ee164e20 (diff)
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.
Diffstat (limited to 'packages/base/src/Internal/LAPACK.hs')
-rw-r--r--packages/base/src/Internal/LAPACK.hs51
1 files changed, 51 insertions, 0 deletions
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