summaryrefslogtreecommitdiff
path: root/lib/Numeric/LinearAlgebra
diff options
context:
space:
mode:
Diffstat (limited to 'lib/Numeric/LinearAlgebra')
-rw-r--r--lib/Numeric/LinearAlgebra/Algorithms.hs27
-rw-r--r--lib/Numeric/LinearAlgebra/Tests.hs9
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
167exactHermitian m = m `equal` ctrans m 168exactHermitian m = m `equal` ctrans m
168 169
170shSize 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.
348invlndet :: (Floating t, Field t)
349 => Matrix t
350 -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det))
351invlndet 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'.
347det :: Field t => Matrix t -> t 363det :: Field t => Matrix t -> t
348det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup) 364det 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)
357lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t) 373lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
358lu = luFact . luPacked 374lu = luFact . luPacked
359 375
360-- | Inverse of a square matrix. 376-- | Inverse of a square matrix. See also 'invlndet'.
361inv :: Field t => Matrix t -> Matrix t 377inv :: Field t => Matrix t -> Matrix t
362inv m | square m = m `linearSolve` ident (rows m) 378inv 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.
366pinv :: Field t => Matrix t -> Matrix t 382pinv :: Field t => Matrix t -> Matrix t
@@ -697,3 +713,4 @@ relativeError :: (Normed c t, Container c t) => c t -> c t -> Int
697relativeError x y = dig (norm (x `sub` y) / norm x) 713relativeError 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
70detTest2 = 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
72polyEval cs x = foldr (\c ac->ac*x+c) 0 cs 80polyEval 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)