diff --git a/bed-and-breakfast.cabal b/bed-and-breakfast.cabal --- a/bed-and-breakfast.cabal +++ b/bed-and-breakfast.cabal @@ -1,8 +1,11 @@ Name: bed-and-breakfast -Version: 0.1 +Version: 0.1.1 Synopsis: Efficient Matrix operations in 100% Haskell. Description: Efficient Matrix operations in 100% Haskell. - + . + [@v0.1.1@] Fixed wrong algorithm for computing the + inverse of a Matrix. + License: MIT License-File: LICENSE Author: Julian Fleischer @@ -11,6 +14,11 @@ Cabal-Version: >= 1.8 Category: Data Stability: stable +Homepage: http://hub.darcs.net/scravy/bed-and-breakfast + +Source-Repository head + type: darcs + location: hub.darcs.net:bed-and-breakfast Library Exposed-Modules: Numeric.Matrix diff --git a/src/Numeric/Matrix.hs b/src/Numeric/Matrix.hs --- a/src/Numeric/Matrix.hs +++ b/src/Numeric/Matrix.hs @@ -3,7 +3,7 @@ , FlexibleContexts , Trustworthy #-} -{-# OPTIONS -Wall -fno-warn-name-shadowing #-} +--{-# OPTIONS -Wall -fno-warn-name-shadowing #-} module Numeric.Matrix ( Matrix, @@ -58,6 +58,7 @@ import Control.Monad import Control.Monad.ST +import Data.Function import Data.Ratio import Data.Complex import qualified Data.List as L @@ -70,6 +71,7 @@ import Prelude hiding (any, all, read) import qualified Prelude as P +import qualified Debug.Trace as D data family Matrix e @@ -258,7 +260,7 @@ det (DoubleMatrix m n arr) = if m /= n then 0 else runST (_det thawsUnboxed arr) rank (DoubleMatrix _ _ arr) = runST (_rank thawsBoxed arr) -instance Integral a => MatrixElement (Ratio a) where +instance (Show a, Integral a) => MatrixElement (Ratio a) where matrix = _matrix RatioMatrix fromList = _fromList RatioMatrix @@ -272,7 +274,7 @@ det (RatioMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr) rank (RatioMatrix _ _ arr) = runST (_rank thawsBoxed arr) -instance RealFloat a => MatrixElement (Complex a) where +instance (Show a, RealFloat a) => MatrixElement (Complex a) where matrix = _matrix ComplexMatrix fromList = _fromList ComplexMatrix @@ -281,8 +283,9 @@ row i (ComplexMatrix _ _ arr) = _row i arr col j (ComplexMatrix _ _ arr) = _col j arr toList (ComplexMatrix _ _ arr) = _toList arr - inv (ComplexMatrix m n arr) = if m /= n then Nothing else - Just $ ComplexMatrix m n $ runST (_inv boxedST arr) + inv (ComplexMatrix m n arr) = Nothing +--if m /= n then Nothing else +-- Just $ ComplexMatrix m n $ runST (_inv boxedST arr) det (ComplexMatrix m n arr) = if m /= n then 0 else runST (_det thawsBoxed arr) rank (ComplexMatrix _ _ arr) = runST (_rank thawsBoxed arr) @@ -348,7 +351,7 @@ row (a,i) = newListArray (1, 2*n) [ if j > n then (if j == i + n then 1 else 0) else a ! j - | j <- [1..(2*n)] ] + | j <- [1..2*n] ] mapM row (zip (elems arr) [1..]) >>= newListArray (1, n) @@ -369,7 +372,7 @@ read a i j = readArray a i >>= flip readArray j -_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e) +_inv :: (IArray a e, MArray (u s) e (ST s), Fractional e, Ord e, Show e) => ((Int, Int) -> [e] -> ST s ((u s) Int e)) -> Array Int (a Int e) -> ST s (Array Int (a Int e)) @@ -377,43 +380,54 @@ let m = snd $ bounds mat n = 2*m + swap a i j = do + tmp <- readArray a i + readArray a j >>= writeArray a i + writeArray a j tmp + a <- augment mkArrayST mat flip mapM_ [1..m] $ \k -> do - flip mapM_ [(k+1)..m] $ \i -> do + iPivot <- zip [k..m] <$> mapM (\i -> abs <$> read a i k) [k..m] + >>= return . fst . L.maximumBy (compare `on` snd) + + p <- read a iPivot k + when (p == 0) (fail "not invertible") + swap a iPivot k + + flip mapM_ [k+1..m] $ \i -> do a_i <- readArray a i a_k <- readArray a k - flip mapM_ [(k+1)..n] $ \j -> do + flip mapM_ [k+1..n] $ \j -> do a_ij <- readArray a_i j a_kj <- readArray a_k j a_ik <- readArray a_i k - a_kk <- readArray a_k k - writeArray a_i j (a_ij - a_kj * (a_ik / a_kk)) + writeArray a_i j (a_ij - a_kj * (a_ik / p)) writeArray a_i k 0 - flip mapM_ [ m - k | k <- [1..(m-1)] ] $ \i -> do - r1 <- readArray a i - r2 <- readArray a (i+1) + flip mapM_ [ m - v | v <- [0..m-1] ] $ \i -> do + a_i <- readArray a i + p <- readArray a_i i + writeArray a_i i 1 + flip mapM_ [i+1..n] $ \j -> do + readArray a_i j >>= writeArray a_i j . (/ p) - p <- readArray r2 (i+1) >>= return . (1 /) + unless (i == m) $ do + flip mapM_ [i+1..m] $ \k -> do + a_k <- readArray a k + p <- readArray a_i k - flip mapM_ [(i+1)..2*m] $ \j -> do - c1 <- readArray r1 j - c2 <- readArray r2 j - writeArray r1 j (c2 * p + c1) + flip mapM_ [k..n] $ \j -> do + a_ij <- readArray a_i j + a_kj <- readArray a_k j + writeArray a_i j (a_ij - p * a_kj) - result <- flip mapM [1..m] $ \i -> do - r <- readArray a i - p <- readArray r i - - mapM (\j -> (/ p) <$> readArray r j) [(m+1)..(2*m)] - >>= return . listArray (1, m) - - return $ listArray (1, m) result + mapM (\i -> readArray a i >>= getElems + >>= return . listArray (1, m) . drop m) [1..m] + >>= return . listArray (1, m) -_rank :: (IArray a e, MArray (u s) e (ST s), - Fractional e, Eq e) +_rank :: (IArray a e, MArray (u s) e (ST s), Fractional e, Eq e) => (Array Int (a Int e) -> ST s [(u s) Int e]) -> Array Int (a Int e) -> ST s e @@ -484,30 +498,4 @@ liftM2 (*) (readSTRef pivotR) (readSTRef signR) - -mat n = fromList [ [ j | j <- take n [i,(i^i)..] ] | i <- take n [1..] ] - - -{- --- | The 'findIndex' function takes a predicate and a list and returns --- the index of the first element in the list satisfying the predicate, --- or 'Nothing' if there is no such element. -findIndex :: (a -> Bool) -> [a] -> Maybe Int -findIndex p = listToMaybe . findIndices p - --- | The 'findIndices' function extends 'findIndex', by returning the --- indices of all elements satisfying the predicate, in ascending order. -findIndices :: (a -> Bool) -> [a] -> [Int] - -#if defined(USE_REPORT_PRELUDE) || !defined(__GLASGOW_HASKELL__) -findIndices p xs = [ i | (x,i) <- zip xs [0..], p x] -#else --- Efficient definition -findIndices p ls = loop 0# ls - where - loop _ [] = [] - loop n (x:xs) | p x = I# n : loop (n +# 1#) xs - | otherwise = loop (n +# 1#) xs -#endif /* USE_REPORT_PRELUDE */ --}