diff options
author | Alberto Ruiz <aruiz@um.es> | 2010-12-19 13:05:31 +0000 |
---|---|---|
committer | Alberto Ruiz <aruiz@um.es> | 2010-12-19 13:05:31 +0000 |
commit | 0acb9b4e88b240d4a2f7aa7ed63b20bd7a332f4f (patch) | |
tree | daf5691df52776a510a8cf627a73b0ebbfa5fcaf /lib | |
parent | 20831d6521e54b42aa8410a1434d55c5b9bc004b (diff) |
invlndet
Diffstat (limited to 'lib')
-rw-r--r-- | lib/Numeric/LinearAlgebra/Algorithms.hs | 27 | ||||
-rw-r--r-- | lib/Numeric/LinearAlgebra/Tests.hs | 9 |
2 files changed, 31 insertions, 5 deletions
diff --git a/lib/Numeric/LinearAlgebra/Algorithms.hs b/lib/Numeric/LinearAlgebra/Algorithms.hs index dd93db2..83464f4 100644 --- a/lib/Numeric/LinearAlgebra/Algorithms.hs +++ b/lib/Numeric/LinearAlgebra/Algorithms.hs | |||
@@ -31,7 +31,8 @@ module Numeric.LinearAlgebra.Algorithms ( | |||
31 | linearSolveLS, | 31 | linearSolveLS, |
32 | linearSolveSVD, | 32 | linearSolveSVD, |
33 | inv, pinv, | 33 | inv, pinv, |
34 | det, rank, rcond, | 34 | det, invlndet, |
35 | rank, rcond, | ||
35 | -- * Matrix factorizations | 36 | -- * Matrix factorizations |
36 | -- ** Singular value decomposition | 37 | -- ** Singular value decomposition |
37 | svd, | 38 | svd, |
@@ -166,6 +167,8 @@ vertical m = rows m >= cols m | |||
166 | 167 | ||
167 | exactHermitian m = m `equal` ctrans m | 168 | exactHermitian m = m `equal` ctrans m |
168 | 169 | ||
170 | shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")" | ||
171 | |||
169 | -------------------------------------------------------------- | 172 | -------------------------------------------------------------- |
170 | 173 | ||
171 | -- | Full singular value decomposition. | 174 | -- | Full singular value decomposition. |
@@ -341,12 +344,25 @@ chol m | exactHermitian m = cholSH m | |||
341 | | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" | 344 | | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix" |
342 | 345 | ||
343 | 346 | ||
347 | -- | Joint computation of inverse and logarithm of determinant of a square matrix. | ||
348 | invlndet :: (Floating t, Field t) | ||
349 | => Matrix t | ||
350 | -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) | ||
351 | invlndet m | square m = (im,(ladm,sdm)) | ||
352 | | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix" | ||
353 | where | ||
354 | lp@(lup,perm) = luPacked m | ||
355 | s = signlp (rows m) perm | ||
356 | dg = toList $ takeDiag $ lup | ||
357 | ladm = sum $ map (log.abs) dg | ||
358 | sdm = s* product (map signum dg) | ||
359 | im = luSolve lp (ident (rows m)) | ||
344 | 360 | ||
345 | 361 | ||
346 | -- | Determinant of a square matrix. | 362 | -- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'. |
347 | det :: Field t => Matrix t -> t | 363 | det :: Field t => Matrix t -> t |
348 | det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup) | 364 | det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup) |
349 | | otherwise = error "det of nonsquare matrix" | 365 | | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix" |
350 | where (lup,perm) = luPacked m | 366 | where (lup,perm) = luPacked m |
351 | s = signlp (rows m) perm | 367 | s = signlp (rows m) perm |
352 | 368 | ||
@@ -357,10 +373,10 @@ det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup) | |||
357 | lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) | 373 | lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) |
358 | lu = luFact . luPacked | 374 | lu = luFact . luPacked |
359 | 375 | ||
360 | -- | Inverse of a square matrix. | 376 | -- | Inverse of a square matrix. See also 'invlndet'. |
361 | inv :: Field t => Matrix t -> Matrix t | 377 | inv :: Field t => Matrix t -> Matrix t |
362 | inv m | square m = m `linearSolve` ident (rows m) | 378 | inv m | square m = m `linearSolve` ident (rows m) |
363 | | otherwise = error "inv of nonsquare matrix" | 379 | | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix" |
364 | 380 | ||
365 | -- | Pseudoinverse of a general matrix. | 381 | -- | Pseudoinverse of a general matrix. |
366 | pinv :: Field t => Matrix t -> Matrix t | 382 | pinv :: Field t => Matrix t -> Matrix t |
@@ -697,3 +713,4 @@ relativeError :: (Normed c t, Container c t) => c t -> c t -> Int | |||
697 | relativeError x y = dig (norm (x `sub` y) / norm x) | 713 | relativeError x y = dig (norm (x `sub` y) / norm x) |
698 | where norm = pnorm Infinity | 714 | where norm = pnorm Infinity |
699 | dig r = round $ -logBase 10 (realToFrac r :: Double) | 715 | dig r = round $ -logBase 10 (realToFrac r :: Double) |
716 | |||
diff --git a/lib/Numeric/LinearAlgebra/Tests.hs b/lib/Numeric/LinearAlgebra/Tests.hs index 093e4ef..68e77f1 100644 --- a/lib/Numeric/LinearAlgebra/Tests.hs +++ b/lib/Numeric/LinearAlgebra/Tests.hs | |||
@@ -67,6 +67,14 @@ detTest1 = det m == 26 | |||
67 | , 2, 8, i | 67 | , 2, 8, i |
68 | ] | 68 | ] |
69 | 69 | ||
70 | detTest2 = inv1 |~| inv2 && [det1] ~~ [det2] | ||
71 | where | ||
72 | m = complex (feye 6) | ||
73 | inv1 = inv m | ||
74 | det1 = det m | ||
75 | (inv2,(lda,sa)) = invlndet m | ||
76 | det2 = sa * exp lda | ||
77 | |||
70 | -------------------------------------------------------------------- | 78 | -------------------------------------------------------------------- |
71 | 79 | ||
72 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs | 80 | polyEval cs x = foldr (\c ac->ac*x+c) 0 cs |
@@ -493,6 +501,7 @@ runTests n = do | |||
493 | _ <- runTestTT $ TestList | 501 | _ <- runTestTT $ TestList |
494 | [ utest "1E5 rots" rotTest | 502 | [ utest "1E5 rots" rotTest |
495 | , utest "det1" detTest1 | 503 | , utest "det1" detTest1 |
504 | , utest "invlndet" detTest2 | ||
496 | , utest "expm1" (expmTest1) | 505 | , utest "expm1" (expmTest1) |
497 | , utest "expm2" (expmTest2) | 506 | , utest "expm2" (expmTest2) |
498 | , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) | 507 | , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM) |