diff options
author | Alberto Ruiz <aruiz@um.es> | 2007-06-11 10:46:39 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2007-06-11 10:46:39 +0000 |
commit | 473df6136476dfa07331dd25a6020260c4f02a9b (patch) | |
tree | 0639081371f7f0d3d03aba2a975921690c19f149 /lib/LAPACK | |
parent | f2cf177e93d4578b404909c68b24625a76466ee5 (diff) |
all eig
Diffstat (limited to 'lib/LAPACK')
-rw-r--r-- | lib/LAPACK/Internal.hs | 72 |
1 files changed, 62 insertions, 10 deletions
diff --git a/lib/LAPACK/Internal.hs b/lib/LAPACK/Internal.hs index e3c9927..ba50e6b 100644 --- a/lib/LAPACK/Internal.hs +++ b/lib/LAPACK/Internal.hs | |||
@@ -79,17 +79,48 @@ eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Dou | |||
79 | eigC (m@M {rows = r}) | 79 | eigC (m@M {rows = r}) |
80 | | r == 1 = (fromList [cdat m `at` 0], singleton 1) | 80 | | r == 1 = (fromList [cdat m `at` 0], singleton 1) |
81 | | otherwise = unsafePerformIO $ do | 81 | | otherwise = unsafePerformIO $ do |
82 | l <- createVector r | 82 | l <- createVector r |
83 | v <- createMatrix ColumnMajor r r | 83 | v <- createMatrix ColumnMajor r r |
84 | dummy <- createMatrix ColumnMajor 1 1 | 84 | dummy <- createMatrix ColumnMajor 1 1 |
85 | zgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigC" [fdat m] | 85 | zgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigC" [fdat m] |
86 | return (l,v) | 86 | return (l,v) |
87 | 87 | ||
88 | ----------------------------------------------------------------------------- | 88 | ----------------------------------------------------------------------------- |
89 | -- dgeev | 89 | -- dgeev |
90 | foreign import ccall "lapack-aux.h eig_l_R" | 90 | foreign import ccall "lapack-aux.h eig_l_R" |
91 | dgeev :: Double ::> Double ::> ((Complex Double) :> Double ::> IO Int) | 91 | dgeev :: Double ::> Double ::> ((Complex Double) :> Double ::> IO Int) |
92 | 92 | ||
93 | -- | Wrapper for LAPACK's /dgeev/, which computes the eigenvalues and right eigenvectors of a general real matrix: | ||
94 | -- | ||
95 | -- if @(l,v)=eigR m@ then @m \<\> v = v \<\> diag l@. | ||
96 | -- | ||
97 | -- The eigenvectors are the columns of v. | ||
98 | -- The eigenvalues are not sorted. | ||
99 | eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) | ||
100 | eigR (m@M {rows = r}) = (s', v'') | ||
101 | where (s,v) = eigRaux m | ||
102 | s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s)) | ||
103 | v' = toRows $ trans v | ||
104 | v'' = fromColumns $ fixeig (toList s') v' | ||
105 | |||
106 | eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) | ||
107 | eigRaux (m@M {rows = r}) | ||
108 | | r == 1 = (fromList [(cdat m `at` 0):+0], singleton 1) | ||
109 | | otherwise = unsafePerformIO $ do | ||
110 | l <- createVector r | ||
111 | v <- createMatrix ColumnMajor r r | ||
112 | dummy <- createMatrix ColumnMajor 1 1 | ||
113 | dgeev // mat fdat m // mat dat dummy // vec l // mat dat v // check "eigR" [fdat m] | ||
114 | return (l,v) | ||
115 | |||
116 | fixeig [] _ = [] | ||
117 | fixeig [r] [v] = [comp v] | ||
118 | fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | ||
119 | | r1 == r2 && i1 == (-i2) = toComplex (v1,v2) : toComplex (v1,scale (-1) v2) : fixeig r vs | ||
120 | | otherwise = comp v1 : fixeig ((r2:+i2):r) (v2:vs) | ||
121 | |||
122 | scale r v = fromList [r] `outer` v | ||
123 | |||
93 | ----------------------------------------------------------------------------- | 124 | ----------------------------------------------------------------------------- |
94 | -- dsyev | 125 | -- dsyev |
95 | foreign import ccall "lapack-aux.h eig_l_S" | 126 | foreign import ccall "lapack-aux.h eig_l_S" |
@@ -106,17 +137,38 @@ eigS m = (s', fliprl v) | |||
106 | where (s,v) = eigS' m | 137 | where (s,v) = eigS' m |
107 | s' = fromList . reverse . toList $ s | 138 | s' = fromList . reverse . toList $ s |
108 | 139 | ||
109 | eigS' (m@M {rows = r}) = unsafePerformIO $ do | 140 | eigS' (m@M {rows = r}) |
110 | l <- createVector r | 141 | | r == 1 = (fromList [cdat m `at` 0], singleton 1) |
111 | v <- createMatrix ColumnMajor r r | 142 | | otherwise = unsafePerformIO $ do |
112 | dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m] | 143 | l <- createVector r |
113 | return (l,v) | 144 | v <- createMatrix ColumnMajor r r |
145 | dsyev // mat fdat m // vec l // mat dat v // check "eigS" [fdat m] | ||
146 | return (l,v) | ||
114 | 147 | ||
115 | ----------------------------------------------------------------------------- | 148 | ----------------------------------------------------------------------------- |
116 | -- zheev | 149 | -- zheev |
117 | foreign import ccall "lapack-aux.h eig_l_H" | 150 | foreign import ccall "lapack-aux.h eig_l_H" |
118 | zheev :: (Complex Double) ::> (Double :> (Complex Double) ::> IO Int) | 151 | zheev :: (Complex Double) ::> (Double :> (Complex Double) ::> IO Int) |
119 | 152 | ||
153 | -- | Wrapper for LAPACK's /zheev/, which computes the eigenvalues and right eigenvectors of a hermitian complex matrix: | ||
154 | -- | ||
155 | -- if @(l,v)=eigH m@ then @m \<\> s v = v \<\> diag l@. | ||
156 | -- | ||
157 | -- The eigenvectors are the columns of v. | ||
158 | -- The eigenvalues are sorted in descending order. | ||
159 | eigH :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | ||
160 | eigH m = (s', fliprl v) | ||
161 | where (s,v) = eigH' m | ||
162 | s' = fromList . reverse . toList $ s | ||
163 | |||
164 | eigH' (m@M {rows = r}) | ||
165 | | r == 1 = (fromList [realPart (cdat m `at` 0)], singleton 1) | ||
166 | | otherwise = unsafePerformIO $ do | ||
167 | l <- createVector r | ||
168 | v <- createMatrix ColumnMajor r r | ||
169 | zheev // mat fdat m // vec l // mat dat v // check "eigH" [fdat m] | ||
170 | return (l,v) | ||
171 | |||
120 | ----------------------------------------------------------------------------- | 172 | ----------------------------------------------------------------------------- |
121 | -- dgesv | 173 | -- dgesv |
122 | foreign import ccall "lapack-aux.h linearSolveR_l" | 174 | foreign import ccall "lapack-aux.h linearSolveR_l" |