invlndet

4 files changed, 37 insertions, 6 deletions

diff --git a/CHANGES b/CHANGES
index cc05934..d9dddc2 100644
--- a/CHANGES
+++ b/CHANGES
@@ -1,3 +1,8 @@
+0.10.1.0
+========
+
+- invlndet (inverse and log of determinant)
+
 0.10.0.0
 ========
 
diff --git a/hmatrix.cabal b/hmatrix.cabal
index a187599..b8f23a1 100644
--- a/hmatrix.cabal
+++ b/hmatrix.cabal
@@ -1,5 +1,5 @@
 Name:               hmatrix
-Version:            0.10.0.2
+Version:            0.10.1.0
 License:            GPL
 License-file:       LICENSE
 Author:             Alberto Ruiz
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 (
     linearSolveLS,
     linearSolveSVD,
     inv, pinv,
-    det, rank, rcond,
+    det, invlndet,
+    rank, rcond,
 -- * Matrix factorizations
 -- ** Singular value decomposition
     svd,
@@ -166,6 +167,8 @@ vertical m = rows m >= cols m
 
 exactHermitian m = m `equal` ctrans m
 
+shSize m = "(" ++ show (rows m) ++"><"++ show (cols m)++")"
+
 --------------------------------------------------------------
 
 -- | Full singular value decomposition.
@@ -341,12 +344,25 @@ chol m | exactHermitian m = cholSH m
        | otherwise = error "chol requires positive definite complex hermitian or real symmetric matrix"
 
 
+-- | Joint computation of inverse and logarithm of determinant of a square matrix.
+invlndet :: (Floating t, Field t)
+         => Matrix t
+         -> (Matrix t, (t, t)) -- ^ (inverse, (log abs det, sign or phase of det)) 
+invlndet m | square m = (im,(ladm,sdm))
+           | otherwise = error $ "invlndet of nonsquare "++ shSize m ++ " matrix"
+  where
+    lp@(lup,perm) = luPacked m
+    s = signlp (rows m) perm
+    dg = toList $ takeDiag $ lup
+    ladm = sum $ map (log.abs) dg
+    sdm = s* product (map signum dg)
+    im = luSolve lp (ident (rows m))
 
 
--- | Determinant of a square matrix.
+-- | Determinant of a square matrix. To avoid possible overflow or underflow use 'invlndet'.
 det :: Field t => Matrix t -> t
 det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)
-      | otherwise = error "det of nonsquare matrix"
+      | otherwise = error $ "det of nonsquare "++ shSize m ++ " matrix"
     where (lup,perm) = luPacked m
           s = signlp (rows m) perm
 
@@ -357,10 +373,10 @@ det m | square m = {-# SCC "det" #-} s * (product $ toList $ takeDiag $ lup)
 lu :: Field t => Matrix t -> (Matrix t, Matrix t, Matrix t, t)
 lu = luFact . luPacked
 
--- | Inverse of a square matrix.
+-- | Inverse of a square matrix. See also 'invlndet'.
 inv :: Field t => Matrix t -> Matrix t
 inv m | square m = m `linearSolve` ident (rows m)
-      | otherwise = error "inv of nonsquare matrix"
+      | otherwise = error $ "inv of nonsquare "++ shSize m ++ " matrix"
 
 -- | Pseudoinverse of a general matrix.
 pinv :: Field t => Matrix t -> Matrix t
@@ -697,3 +713,4 @@ relativeError :: (Normed c t, Container c t) => c t -> c t -> Int
 relativeError x y = dig (norm (x `sub` y) / norm x)
     where norm = pnorm Infinity
           dig r = round $ -logBase 10 (realToFrac r :: Double)
+
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
             , 2, 8, i
             ]
 
+detTest2 = inv1 |~| inv2 && [det1] ~~ [det2]
+  where
+    m = complex (feye 6)
+    inv1 = inv m
+    det1 = det m
+    (inv2,(lda,sa)) = invlndet m
+    det2 = sa * exp lda
+
 --------------------------------------------------------------------
 
 polyEval cs x = foldr (\c ac->ac*x+c) 0 cs
@@ -493,6 +501,7 @@ runTests n = do
     _ <- runTestTT $ TestList
         [ utest "1E5 rots" rotTest
         , utest "det1" detTest1
+        , utest "invlndet" detTest2
         , utest "expm1" (expmTest1)
         , utest "expm2" (expmTest2)
         , utest "arith1" $ ((ones (100,100) * 5 + 2)/0.5 - 7)**2 |~| (49 :: RM)