Added code in chapter

Author: jiegillet

chapters/matrix_methods/gaussian_elimination/code/haskell/gaussianElimination.hs

@@ -18,8 +18,9 @@ swapRows r1 r2 m
 subRows :: Fractional a =>
            (Int, Int) -> (Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
 subRows (r, c) (r1, rn) (c1, cn) m =
-  accum (-) m [ ((i, j), m!(i, c) * m!(r,j)/m!(r, c))
-              | i <- [r1..rn], j <- [c1..cn]
+  accum (-) m [ ((i, j), m!(i, c) * m!(r,j) / m!(r, c))
+              | i <- [r1..rn]
+              , j <- [c1..cn]
               ]
 
 toEchelon :: (Fractional a, Ord a) => Matrix a -> Matrix a

contents/gaussian_elimination/code/haskell/gaussianElimination.hs

@@ -0,0 +1,75 @@
+import Data.Array
+import Data.List (intercalate, maximumBy)
+import Data.Function (on)
+import Data.Ratio
+
+type Matrix a = Array (Int, Int) a
+
+type Vector a = Array Int a
+
+swapRows :: Int -> Int -> Matrix a -> Matrix a
+swapRows r1 r2 m
+  | r1 == r2  = m
+  | otherwise = m // concat [ [((r2, c), m!(r1, c)), ((r1, c), m!(r2, c))]
+                            | c <- [c1..cn]
+                            ]
+  where ((_, c1), (_, cn)) = bounds m
+
+subRows :: Fractional a =>
+           (Int, Int) -> (Int, Int) -> (Int, Int) -> Matrix a -> Matrix a
+subRows (r, c) (r1, rn) (c1, cn) m =
+  accum (-) m [ ((i, j), m!(i, c) * m!(r,j) / m!(r, c))
+              | i <- [r1..rn]
+              , j <- [c1..cn]
+              ]
+
+gaussianElimination :: (Fractional a, Ord a) => Matrix a -> Matrix a
+gaussianElimination mat = go (r1, c1) mat
+  where
+  ((r1, c1), (rn, cn)) = bounds mat
+  go (r, c) m
+    | c == cn    = m
+    | pivot == 0 = go (r, c + 1) m
+    | otherwise  = go (r + 1, c + 1) $ subRows (r, c) (r+1, rn) (c, cn) m'
+    where (target, pivot)
+            = maximumBy (compare `on` abs . snd) [ (k, m!(k, c)) | k <- [r..rn]]
+          m' = swapRows r target m
+
+gaussJordan :: (Fractional a, Eq a) => Matrix a -> Matrix a
+gaussJordan mat = go (r1, c1) mat
+  where
+  ((r1, c1), (rn, cn)) = bounds mat
+  go (r, c) m
+    | c == cn       = m
+    | m!(r, c) == 0 = go (r, c + 1) m
+    | otherwise     = go (r + 1, c + 1) $ subRows (r, c) (r1, r-1) (c, cn) m'
+    where m' = accum (/) m [((r, j), m!(r, c)) | j <- [c..cn] ]
+
+backSubstitution :: (Fractional a) => Matrix a -> Vector a
+backSubstitution m = sol
+  where
+  ((r1, c1), (rn, cn)) = bounds m
+  sol = listArray (r1,rn) [ (m!(r, cn) - sum' r) / m!(r, r) | r <- [r1..rn]]
+  sum' r = sum [ m!(r, k) * sol!k | k <- [r+1..rn]]
+
+printM :: (Show a) => Matrix a -> String
+printM m =
+  let  ((r1, c1), (rn, cn)) = bounds m
+  in unlines [ intercalate "\t" [ show $ m!(r, c) | c <- [c1..cn] ]
+             | r <- [r1..rn]
+             ]
+
+printV :: (Show a) => Vector a -> String
+printV = unlines . map show . elems
+
+main :: IO ()
+main = do
+  let m = listArray ((1,1),(3,4)) [2,3,4,6,1,2,3,4,3,-4,0,10] :: Matrix (Ratio Int)
+  putStrLn "Original Matrix:"
+  putStrLn $ printM m
+  putStrLn "Echelon form"
+  putStrLn $ printM $ gaussianElimination m
+  putStrLn "Reduced echelon form"
+  putStrLn $ printM $ gaussJordan $ gaussianElimination m
+  putStrLn "Solution from back substitution"
+  putStrLn $ printV $ backSubstitution $ gaussianElimination m

contents/gaussian_elimination/gaussian_elimination.md

@@ -149,7 +149,7 @@ $$
 \begin{array}{ccc|c}
 1 & 0 & 0 & \frac{18}{11} \\
 0 & 1 & 0 & \frac{-14}{11} \\
-0 & 0 & 1 & \frac{18}{11} 
+0 & 0 & 1 & \frac{18}{11}
 \end{array}
 \right]
 $$
@@ -358,6 +358,8 @@ In code, this looks like:
 [import:13-44, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import:41-78, lang:"rust"](code/rust/gaussian_elimination.rs)
+{% sample lang="hs" %}
+[import:10-36, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 Now, to be clear: this algorithm creates an upper-triangular matrix.
@@ -390,6 +392,8 @@ This code does not exist yet in C, so here's Julia code (sorry for the inconveni
 {% sample lang="rs" %}
 This code does not exist yet in rust, so here's Julia code (sorry for the inconvenience)
 [import:70-96, lang:"julia"](code/julia/gaussian_elimination.jl)
+{% sample lang="hs" %}
+[import:38-46, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 ## Back-substitution
@@ -418,6 +422,8 @@ In code, this involves keeping a rolling sum of all the values we substitute in
 [import:46-58, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import:79-94, lang:"rust"](code/rust/gaussian_elimination.rs)
+{% sample lang="hs" %}
+[import:48-53, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 ## Conclusions
@@ -438,6 +444,8 @@ As for what's next... Well, we are in for a treat! The above algorithm clearly h
 [import, lang:"c_cpp"](code/c/gaussian_elimination.c)
 {% sample lang="rs" %}
 [import, lang:"rust"](code/rust/gaussian_elimination.rs)
+{% sample lang="hs" %}
+[import, lang:"haskell"](code/haskell/gaussianElimination.hs)
 {% endmethod %}
 
 
@@ -463,4 +471,3 @@ $$
 \newcommand{\bfomega}{\boldsymbol{\omega}}
 \newcommand{\bftau}{\boldsymbol{\tau}}
 $$
-