{-
   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