diff options
Diffstat (limited to 'packages/base/src/Internal/LAPACK.hs')
-rw-r--r-- | packages/base/src/Internal/LAPACK.hs | 10 |
1 files changed, 9 insertions, 1 deletions
diff --git a/packages/base/src/Internal/LAPACK.hs b/packages/base/src/Internal/LAPACK.hs index b834b20..c9f5d0f 100644 --- a/packages/base/src/Internal/LAPACK.hs +++ b/packages/base/src/Internal/LAPACK.hs | |||
@@ -302,6 +302,14 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | |||
302 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) | 302 | | otherwise = comp' v1 : fixeig ((r2:+i2):r) (v2:vs) |
303 | fixeig _ _ = error "fixeig with impossible inputs" | 303 | fixeig _ _ = error "fixeig with impossible inputs" |
304 | 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 ((ar1:+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) | ||
305 | 313 | ||
306 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. | 314 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. |
307 | -- The eigenvalues are not sorted. | 315 | -- The eigenvalues are not sorted. |
@@ -316,7 +324,7 @@ eigG a b = (alpha', beta, v'') | |||
316 | (alpha, beta, v) = eigGaux dggev a b "eigG" | 324 | (alpha, beta, v) = eigGaux dggev a b "eigG" |
317 | alpha' = fixeig1 alpha | 325 | alpha' = fixeig1 alpha |
318 | v' = toRows $ trans v | 326 | v' = toRows $ trans v |
319 | v'' = fromColumns $ fixeig (toList alpha') v' | 327 | v'' = fromColumns $ fixeigG (toList alpha') v' |
320 | 328 | ||
321 | eigGaux f ma mb st = unsafePerformIO $ do | 329 | eigGaux f ma mb st = unsafePerformIO $ do |
322 | a <- copy ColumnMajor ma | 330 | a <- copy ColumnMajor ma |