diff options
Diffstat (limited to 'packages/base/src/Internal/LAPACK.hs')
-rw-r--r-- | packages/base/src/Internal/LAPACK.hs | 60 |
1 files changed, 60 insertions, 0 deletions
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 | ||
19 | module Internal.LAPACK where | 19 | module Internal.LAPACK where |
20 | 20 | ||
21 | import Data.Bifunctor (first) | ||
22 | |||
21 | import Internal.Devel | 23 | import Internal.Devel |
22 | import Internal.Vector | 24 | import Internal.Vector |
23 | import Internal.Matrix hiding ((#), (#!)) | 25 | import Internal.Matrix hiding ((#), (#!)) |
@@ -234,7 +236,9 @@ leftSVAux f st x = unsafePerformIO $ do | |||
234 | ----------------------------------------------------------------------------- | 236 | ----------------------------------------------------------------------------- |
235 | 237 | ||
236 | foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok | 238 | foreign import ccall unsafe "eig_l_R" dgeev :: R ::> R ::> C :> R ::> Ok |
239 | foreign import ccall unsafe "eig_l_G" dggev :: R ::> R ::> C :> R :> R ::> R ::> Ok | ||
237 | foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok | 240 | foreign import ccall unsafe "eig_l_C" zgeev :: C ::> C ::> C :> C ::> Ok |
241 | foreign import ccall unsafe "eig_l_GC" zggev :: C ::> C ::> C :> C :> C ::> C ::> Ok | ||
238 | foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok | 242 | foreign import ccall unsafe "eig_l_S" dsyev :: CInt -> R :> R ::> Ok |
239 | foreign import ccall unsafe "eig_l_H" zheev :: CInt -> R :> C ::> Ok | 243 | foreign 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) |
299 | fixeig _ _ = error "fixeig with impossible inputs" | 303 | fixeig _ _ = 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'. | ||
308 | fixeigG [] _ = [] | ||
309 | fixeigG [_] [v] = [comp' v] | ||
310 | fixeigG ((_:+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) | ||
313 | fixeigG _ _ = 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. |
304 | eigOnlyR :: Matrix Double -> Vector (Complex Double) | 317 | eigOnlyR :: Matrix Double -> Vector (Complex Double) |
305 | eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" | 318 | eigOnlyR = 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. | ||
322 | eigG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double, Matrix (Complex Double)) | ||
323 | eigG 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 | |||
330 | eigGaux 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 | |||
344 | eigGOnlyAux 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. | ||
359 | eigGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double), Matrix (Complex Double)) | ||
360 | eigGC a b = eigGaux zggev a b "eigGC" | ||
361 | |||
362 | eigOnlyG :: Matrix Double -> Matrix Double -> (Vector (Complex Double), Vector Double) | ||
363 | eigOnlyG a b = first fixeig1 $ eigGOnlyAux dggev a b "eigOnlyG" | ||
364 | |||
365 | eigOnlyGC :: Matrix (Complex Double) -> Matrix (Complex Double) -> (Vector (Complex Double), Vector (Complex Double)) | ||
366 | eigOnlyGC a b = eigGOnlyAux zggev a b "eigOnlyGC" | ||
307 | 367 | ||
308 | ----------------------------------------------------------------------------- | 368 | ----------------------------------------------------------------------------- |
309 | 369 | ||