We want to program the gauss-elimination to calculate a basis (linear algebra) as exercise for ourselves. It is not homework.
I thought first of [[Int]]
as structure for our matrix. I thought then that we can sort the lists lexicographically. But then we must calculate with the matrix. And there is the problem. Can someone give us some hints.
Consider using matrices from the hmatrix package. Among its modules you can find both a fast implementation of a matrix and a lot of linear algebra algorithms. Browsing their sources might help you with your doubts.
Here's a simple example of adding one row to another by splitting the matrix into rows.
import Numeric.Container
import Data.Packed.Matrix
addRow :: Container Vector t => Int -> Int -> Matrix t -> Matrix t
addRow from to m = let rows = toRows m in
fromRows $ take to rows ++
[(rows !! from) `add` (rows !! to)] ++
drop (to + 1) rows
Another example, this time by using matrix multiplication.
addRow :: (Product e, Container Vector e) =>
Int -> Int -> Matrix e -> Matrix e
addRow from to m = m `add` (e <> m)
where
nrows = rows m
e = buildMatrix nrows nrows
(\(r,c) -> if (r,c) /= (to,from) then 0 else 1)
Cf. Container, Vector, Product.
It will be easier if you use [[Rational]]
instead of [[Int]]
since you get nice division.
You probably want to start by implementing the elementary row operations.
swap :: Int -> Int -> [[Rational]] -> [[Rational]
swap r1 r2 m = --a matrix with r1 and r2 swapped
scale :: Int -> Rational -> [[Rational]] -> [[Rational]]
scale r c m = --a matrix with row r multiplied by c
addrow :: Int -> Int -> Rational -> [[Rational]] -> [[Rational]]
addrow r1 r2 c m = --a matrix with (c * r1) added to r2
In order to actually do guassian elimination, you need a way to decide what multiple of one row to add to another to get a zero. So given two rows..
5 4 3 2 1
7 6 5 4 3
We want to add c times row 1 to row 2 so that the 7 becomes a zero. So 7 + c * 5 = 0
and c = -7/5
. So in order to solve for c all we need are the first elements of each row. Here's a function that finds c:
whatc :: Rational -> Rational -> Rational
whatc _ 0 = 0
whatc a b = - a / b
Also, as others have said, using lists to represent your matrix will give you worse performance. But if you're just trying to understand the algorithm, lists should be fine.