instead of (!) (or readArray in stateful code) helped a lot. So then I
went to optimize my gaussian elimination function and found just the
opposite. unsafeRead is slower than readArray. This struck me as very
The functions gaussElimSafe and gaussElimUnsafe are different in their
outputs (*). I guess you made a mistake during conversion from one
version to the other.
BTW, the attached file contains an outline of a pure functional and
indexless implementation. It calculates the inverse of a quadratic
matrix in time O(n^3), so maybe you can get rid of all this exclamation
marks and IO stuff...
BR,
(*) Try:
newStdGen >>= \ g -> flip mapM_ [gaussElimSafe,gaussElimUnsafe] $ \ f ->
makeMatrix g >>= \ m -> printMatrix m >> f m >> printMatrix m
--
-- Mirko Rahn -- Tel +49-721 608 7504 --
--- http://liinwww.ira.uka.de/~rahn/ ---
import Control.Monad (mplus,liftM)
type Matrix a = [[a]]
inverse :: (Monad m, Fractional a) => Matrix a -> m (Matrix a)
inverse m
| null m = fail "empty matrix"
| any (/=w) ws = fail "matrix not rectangular"
| h /= w = fail "matrix not quadractic"
| otherwise = liftM (map (drop h)) . clean . gauss $ e
where h = length m
w:ws = map length m
e = zipWith (++) m $ munit h
munit n =
let k = n-1
rep = flip replicate 0
in [0..k] >>= \ i -> [rep i ++ [1] ++ rep (k-i)]
gauss m =
flip (maybe m) (pivot id m) $ \ (f,row:rows) ->
map f . (row:) . map (0:) . gauss . map tail $ rows
pivot adjust m
| all null m = fail "no pivot at all"
| otherwise = mplus (pivot1 m >>= return . (,) adjust)
(pivot ((0:) . adjust) (map tail m))
pivot1 (row@(0:_):rows) = pivot1 rows >>= \ (r:rs) -> return (r:row:rs)
pivot1 (row@(p:_):rows) = return . (first:) . map add $ rows
where first = map (/p) row
add r@(x:_) = sub r x first
pivot1 _ = fail "no pivot in first colum"
sub v c = zipWith (-) v . map (*c)
clean (row@(1:rt):rows) =
liftM (map (0:)) (clean . map tail $ rows) >>= \ crows ->
return $ (foldl ( \ v (a,w) -> sub v a w ) row $ zip rt crows) : crows
clean ( (_: _):_ ) = fail "matrix singular"
clean m = return m
_______________________________________________
Glasgow-haskell-users mailing list
Glasgow-haskell-users@haskell.org
http://www.haskell.org/mailman/listinfo/glasgow-haskell-users
