diff options
Diffstat (limited to 'lib/Numeric/LinearAlgebra/LAPACK.hs')
-rw-r--r-- | lib/Numeric/LinearAlgebra/LAPACK.hs | 85 |
1 files changed, 54 insertions, 31 deletions
diff --git a/lib/Numeric/LinearAlgebra/LAPACK.hs b/lib/Numeric/LinearAlgebra/LAPACK.hs index 2eae9b7..a1ac1cf 100644 --- a/lib/Numeric/LinearAlgebra/LAPACK.hs +++ b/lib/Numeric/LinearAlgebra/LAPACK.hs | |||
@@ -27,6 +27,7 @@ module Numeric.LinearAlgebra.LAPACK ( | |||
27 | rightSVR, rightSVC, leftSVR, leftSVC, | 27 | rightSVR, rightSVC, leftSVR, leftSVC, |
28 | -- * Eigensystems | 28 | -- * Eigensystems |
29 | eigR, eigC, eigS, eigS', eigH, eigH', | 29 | eigR, eigC, eigS, eigS', eigH, eigH', |
30 | eigOnlyR, eigOnlyC, eigOnlyS, eigOnlyH, | ||
30 | -- * LU | 31 | -- * LU |
31 | luR, luC, | 32 | luR, luC, |
32 | -- * Cholesky | 33 | -- * Cholesky |
@@ -198,18 +199,16 @@ leftSVAux f st x = unsafePerformIO $ do | |||
198 | 199 | ||
199 | foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM | 200 | foreign import ccall "LAPACK/lapack-aux.h eig_l_R" dgeev :: TMMCVM |
200 | foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM | 201 | foreign import ccall "LAPACK/lapack-aux.h eig_l_C" zgeev :: TCMCMCVCM |
201 | foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: TMVM | 202 | foreign import ccall "LAPACK/lapack-aux.h eig_l_S" dsyev :: CInt -> TMVM |
202 | foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: TCMVCM | 203 | foreign import ccall "LAPACK/lapack-aux.h eig_l_H" zheev :: CInt -> TCMVCM |
203 | 204 | ||
204 | eigAux f st m | 205 | eigAux f st m = unsafePerformIO $ do |
205 | | r == 1 = (fromList [flatten m `at` 0], singleton 1) | ||
206 | | otherwise = unsafePerformIO $ do | ||
207 | l <- createVector r | 206 | l <- createVector r |
208 | v <- createMatrix ColumnMajor r r | 207 | v <- createMatrix ColumnMajor r r |
209 | dummy <- createMatrix ColumnMajor 1 1 | 208 | app3 g mat m vec l mat v st |
210 | app4 f mat m mat dummy vec l mat v st | ||
211 | return (l,v) | 209 | return (l,v) |
212 | where r = rows m | 210 | where r = rows m |
211 | g ra ca pa = f ra ca pa 0 0 nullPtr | ||
213 | 212 | ||
214 | 213 | ||
215 | -- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. | 214 | -- | Eigenvalues and right eigenvectors of a general complex matrix, using LAPACK's /zgeev/. |
@@ -217,26 +216,39 @@ eigAux f st m | |||
217 | eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) | 216 | eigC :: Matrix (Complex Double) -> (Vector (Complex Double), Matrix (Complex Double)) |
218 | eigC = eigAux zgeev "eigC" . fmat | 217 | eigC = eigAux zgeev "eigC" . fmat |
219 | 218 | ||
219 | eigOnlyAux f st m = unsafePerformIO $ do | ||
220 | l <- createVector r | ||
221 | app2 g mat m vec l st | ||
222 | return l | ||
223 | where r = rows m | ||
224 | g ra ca pa nl pl = f ra ca pa 0 0 nullPtr nl pl 0 0 nullPtr | ||
225 | |||
226 | -- | Eigenvalues of a general complex matrix, using LAPACK's /zgeev/ with jobz == \'N\'. | ||
227 | -- The eigenvalues are not sorted. | ||
228 | eigOnlyC :: Matrix (Complex Double) -> Vector (Complex Double) | ||
229 | eigOnlyC = eigOnlyAux zgeev "eigOnlyC" . fmat | ||
230 | |||
220 | -- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. | 231 | -- | Eigenvalues and right eigenvectors of a general real matrix, using LAPACK's /dgeev/. |
221 | -- The eigenvectors are the columns of v. The eigenvalues are not sorted. | 232 | -- The eigenvectors are the columns of v. The eigenvalues are not sorted. |
222 | eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) | 233 | eigR :: Matrix Double -> (Vector (Complex Double), Matrix (Complex Double)) |
223 | eigR m = (s', v'') | 234 | eigR m = (s', v'') |
224 | where (s,v) = eigRaux (fmat m) | 235 | where (s,v) = eigRaux (fmat m) |
225 | s' = toComplex (subVector 0 r (asReal s), subVector r r (asReal s)) | 236 | s' = fixeig1 s |
226 | v' = toRows $ trans v | 237 | v' = toRows $ trans v |
227 | v'' = fromColumns $ fixeig (toList s') v' | 238 | v'' = fromColumns $ fixeig (toList s') v' |
228 | r = rows m | 239 | r = rows m |
229 | 240 | ||
230 | eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) | 241 | eigRaux :: Matrix Double -> (Vector (Complex Double), Matrix Double) |
231 | eigRaux m | 242 | eigRaux m = unsafePerformIO $ do |
232 | | r == 1 = (fromList [(flatten m `at` 0):+0], singleton 1) | ||
233 | | otherwise = unsafePerformIO $ do | ||
234 | l <- createVector r | 243 | l <- createVector r |
235 | v <- createMatrix ColumnMajor r r | 244 | v <- createMatrix ColumnMajor r r |
236 | dummy <- createMatrix ColumnMajor 1 1 | 245 | app3 g mat m vec l mat v "eigR" |
237 | app4 dgeev mat m mat dummy vec l mat v "eigR" | ||
238 | return (l,v) | 246 | return (l,v) |
239 | where r = rows m | 247 | where r = rows m |
248 | g ra ca pa = dgeev ra ca pa 0 0 nullPtr | ||
249 | |||
250 | fixeig1 s = toComplex (subVector 0 r (asReal s), subVector r r (asReal s)) | ||
251 | where r = dim s | ||
240 | 252 | ||
241 | fixeig [] _ = [] | 253 | fixeig [] _ = [] |
242 | fixeig [_] [v] = [comp v] | 254 | fixeig [_] [v] = [comp v] |
@@ -246,8 +258,22 @@ fixeig ((r1:+i1):(r2:+i2):r) (v1:v2:vs) | |||
246 | where scale = vectorMapValR Scale | 258 | where scale = vectorMapValR Scale |
247 | fixeig _ _ = error "fixeig with impossible inputs" | 259 | fixeig _ _ = error "fixeig with impossible inputs" |
248 | 260 | ||
261 | |||
262 | -- | Eigenvalues of a general real matrix, using LAPACK's /dgeev/ with jobz == \'N\'. | ||
263 | -- The eigenvalues are not sorted. | ||
264 | eigOnlyR :: Matrix Double -> Vector (Complex Double) | ||
265 | eigOnlyR = fixeig1 . eigOnlyAux dgeev "eigOnlyR" . fmat | ||
266 | |||
267 | |||
249 | ----------------------------------------------------------------------------- | 268 | ----------------------------------------------------------------------------- |
250 | 269 | ||
270 | eigSHAux f st m = unsafePerformIO $ do | ||
271 | l <- createVector r | ||
272 | v <- createMatrix ColumnMajor r r | ||
273 | app3 f mat m vec l mat v st | ||
274 | return (l,v) | ||
275 | where r = rows m | ||
276 | |||
251 | -- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. | 277 | -- | Eigenvalues and right eigenvectors of a symmetric real matrix, using LAPACK's /dsyev/. |
252 | -- The eigenvectors are the columns of v. | 278 | -- The eigenvectors are the columns of v. |
253 | -- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). | 279 | -- The eigenvalues are sorted in descending order (use 'eigS'' for ascending order). |
@@ -258,16 +284,7 @@ eigS m = (s', fliprl v) | |||
258 | 284 | ||
259 | -- | 'eigS' in ascending order | 285 | -- | 'eigS' in ascending order |
260 | eigS' :: Matrix Double -> (Vector Double, Matrix Double) | 286 | eigS' :: Matrix Double -> (Vector Double, Matrix Double) |
261 | eigS' m | 287 | eigS' = eigSHAux (dsyev 1) "eigS'" . fmat |
262 | | r == 1 = (fromList [flatten m `at` 0], singleton 1) | ||
263 | | otherwise = unsafePerformIO $ do | ||
264 | l <- createVector r | ||
265 | v <- createMatrix ColumnMajor r r | ||
266 | app3 dsyev mat m vec l mat v "eigS" | ||
267 | return (l,v) | ||
268 | where r = rows m | ||
269 | |||
270 | ----------------------------------------------------------------------------- | ||
271 | 288 | ||
272 | -- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. | 289 | -- | Eigenvalues and right eigenvectors of a hermitian complex matrix, using LAPACK's /zheev/. |
273 | -- The eigenvectors are the columns of v. | 290 | -- The eigenvectors are the columns of v. |
@@ -279,14 +296,20 @@ eigH m = (s', fliprl v) | |||
279 | 296 | ||
280 | -- | 'eigH' in ascending order | 297 | -- | 'eigH' in ascending order |
281 | eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) | 298 | eigH' :: Matrix (Complex Double) -> (Vector Double, Matrix (Complex Double)) |
282 | eigH' m | 299 | eigH' = eigSHAux (zheev 1) "eigH'" . fmat |
283 | | r == 1 = (fromList [realPart (flatten m `at` 0)], singleton 1) | 300 | |
284 | | otherwise = unsafePerformIO $ do | 301 | |
285 | l <- createVector r | 302 | -- | Eigenvalues of a symmetric real matrix, using LAPACK's /dsyev/ with jobz == \'N\'. |
286 | v <- createMatrix ColumnMajor r r | 303 | -- The eigenvalues are sorted in descending order. |
287 | app3 zheev mat m vec l mat v "eigH" | 304 | eigOnlyS :: Matrix Double -> Vector Double |
288 | return (l,v) | 305 | eigOnlyS = vrev . fst. eigSHAux (dsyev 0) "eigS'" . fmat |
289 | where r = rows m | 306 | |
307 | -- | Eigenvalues of a hermitian complex matrix, using LAPACK's /zheev/ with jobz == \'N\'. | ||
308 | -- The eigenvalues are sorted in descending order. | ||
309 | eigOnlyH :: Matrix (Complex Double) -> Vector Double | ||
310 | eigOnlyH = vrev . fst. eigSHAux (zheev 1) "eigH'" . fmat | ||
311 | |||
312 | vrev = flatten . flipud . reshape 1 | ||
290 | 313 | ||
291 | ----------------------------------------------------------------------------- | 314 | ----------------------------------------------------------------------------- |
292 | foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM | 315 | foreign import ccall "LAPACK/lapack-aux.h linearSolveR_l" dgesv :: TMMM |