{-
Haskell Dynamics Engine, Copyright (C) 2007 Ruben Henner Zilibowitz
All rights reserved. Email: rubenz@cse.unsw.edu.au Web: www.cse.unsw.edu.au/~rubenz
This is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License that is include with this
package in the file LICENSE-GPL.TXT.
-}
-------------------------------------------------------------------------------
Maths Specification
Gaussian Elimination
------------------------------
Abstract
Input: Matrix A, vector b
Output: vector x st. Ax = b
------------------------------
Method
Input: Matrix A, vector b
1. Convert to upper triangular
For j = 1 to n-1
For i = j+1 to n
row i -> (-aij/ajj)*row j + row i (row of both a and b)
2. Do back substitution:
For i = n to 1
xi = bi - (ai(i+1)*x(i+1) + ... + ain*xn)
Output: vector x
------------------------------
More concrete version
Input: Matrix A = (aij), vector b = (bi)
1. For j = 1 to n-1
For i = j+1 to n
For k = 1 to n (equivalently, j to n)
aik = (-aij/ajj)*ajk + aik
bi = (-aij/ajj)*bj + bi
2. For i = n to 1
xi = bi - (ai(i+1)*x(i+1) + ... + ain*xn)
Output: vector x
-------------------------------------------------------------------------------
Haskell Specification
> module Gausselim(gauss,Vector,Matrix) where
Types:
Lists for vectors and lists of lists for matrices. For now, we just
use lists - can bother with finite sequences (FinSeq's) later.
> type Vector a = [a]
> type Matrix a = [[a]] -- list of rows
The following function is needed for a for loop in which the ith value
depends on the jth(j>i). The version written in terms of scan is chosen
here because it's easier to transform to
> accum_scanr1 :: (a->[b]->b) -> [a] -> [b]
> accum_scanr1 f zs = map head (scanr' g [] zs)
> where -- g :: a -> [b] -> [b]
> g x ys = (f x ys):ys
> -- scanr' is a form of scanr which returns the same number
> -- of results as the size of the initial list - it misses
> -- out the initial value - important for parallel things
> scanr' f a xs = init(scanr f a xs)
The tidied up Haskell specificaion:
> gauss :: RealFloat a => Matrix a -> Vector a -> Vector a
> gauss a b =
> let m = map2 join a b -- combine a and b into one matrix
> where join xs y = xs ++ [y]
> n = length b
> m2 = foldl deal_with_jth_row m [1..n-1] -- gaussian elimination
> -- for each row j=1..n-1, deal with the jth row
> x = accum_scanr1 solve_ith_eqn m2 -- back substitution
> -- for each eqn i=n,..1, solve that eqn for xi
> in x
> deal_with_jth_row :: RealFloat a => Matrix a -> Int -> Matrix a
> deal_with_jth_row m j = first_j_rows ++
> (map put_0_in_jth_pos other_rows)
> -- makes rows below row j have 0 in cols i<=j
> where first_j_rows = take j m
> other_rows = drop j m
> put_0_in_jth_pos row = map2 join row (m!!(j-1))
> where join x y = multiplier * y + x
> multiplier = -((row!!(j-1)) `safeDivide` (m!!(j-1)!!(j-1)))
> solve_ith_eqn :: RealFloat a => Vector a -> Vector a -> a
> -- calculates the value of xi, using the ith equation.
> -- row is the ith eqn, and zs are previously calculated x values
> solve_ith_eqn row zs =
> (foldl (+) (row!!n) (map2 f (drop i row)zs)) `safeDivide` (row!!(i-1))
> where n = (length row) - 1
> i = n - length zs
> f x y = -(x*y)
> map2 = zipWith
> isBadNum x = (isNaN x) || (isInfinite x)
> safeDivide x y | isBadNum z = error "error solving linear system"
> | otherwise = z
> where z = x / y