|
| 1 | +import Data.Array |
| 2 | +import Data.Function (on) |
| 3 | +import Data.List (intercalate, maximumBy) |
| 4 | +import Data.Ratio |
| 5 | + |
| 6 | +type Matrix a = Array (Int, Int) a |
| 7 | + |
| 8 | +type Vector a = Array Int a |
| 9 | + |
| 10 | +swapRows :: Int -> Int -> Matrix a -> Matrix a |
| 11 | +swapRows r1 r2 m |
| 12 | + | r1 == r2 = m |
| 13 | + | otherwise = |
| 14 | + m // |
| 15 | + concat [[((r2, c), m ! (r1, c)), ((r1, c), m ! (r2, c))] | c <- [c1 .. cn]] |
| 16 | + where |
| 17 | + ((_, c1), (_, cn)) = bounds m |
| 18 | + |
| 19 | +subRows :: |
| 20 | + Fractional a |
| 21 | + => (Int, Int) -- pivot location |
| 22 | + -> (Int, Int) -- rows to cover |
| 23 | + -> (Int, Int) -- columns to cover |
| 24 | + -> Matrix a |
| 25 | + -> Matrix a |
| 26 | +subRows (r, c) (r1, rn) (c1, cn) m = |
| 27 | + accum |
| 28 | + (-) |
| 29 | + m |
| 30 | + [ ((i, j), m ! (i, c) * m ! (r, j) / m ! (r, c)) |
| 31 | + | i <- [r1 .. rn] |
| 32 | + , j <- [c1 .. cn] |
| 33 | + ] |
| 34 | + |
| 35 | +gaussianElimination :: (Fractional a, Ord a) => Matrix a -> Matrix a |
| 36 | +gaussianElimination mat = go (r1, c1) mat |
| 37 | + where |
| 38 | + ((r1, c1), (rn, cn)) = bounds mat |
| 39 | + go (r, c) m |
| 40 | + | c == cn = m |
| 41 | + | pivot == 0 = go (r, c + 1) m |
| 42 | + | otherwise = go (r + 1, c + 1) $ subRows (r, c) (r + 1, rn) (c, cn) m' |
| 43 | + where |
| 44 | + (target, pivot) = |
| 45 | + maximumBy (compare `on` abs . snd) [(k, m ! (k, c)) | k <- [r .. rn]] |
| 46 | + m' = swapRows r target m |
| 47 | + |
| 48 | +gaussJordan :: (Fractional a, Eq a) => Matrix a -> Matrix a |
| 49 | +gaussJordan mat = go (r1, c1) mat |
| 50 | + where |
| 51 | + ((r1, c1), (rn, cn)) = bounds mat |
| 52 | + go (r, c) m |
| 53 | + | c == cn = m |
| 54 | + | m ! (r, c) == 0 = go (r, c + 1) m |
| 55 | + | otherwise = go (r + 1, c + 1) $ subRows (r, c) (r1, r - 1) (c, cn) m' |
| 56 | + where |
| 57 | + m' = accum (/) m [((r, j), m ! (r, c)) | j <- [c .. cn]] |
| 58 | + |
| 59 | +backSubstitution :: (Fractional a) => Matrix a -> Vector a |
| 60 | +backSubstitution m = sol |
| 61 | + where |
| 62 | + ((r1, _), (rn, cn)) = bounds m |
| 63 | + sol = |
| 64 | + listArray (r1, rn) [(m ! (r, cn) - sum' r) / m ! (r, r) | r <- [r1 .. rn]] |
| 65 | + sum' r = sum [m ! (r, k) * sol ! k | k <- [r + 1 .. rn]] |
| 66 | + |
| 67 | +printM :: (Show a) => Matrix a -> String |
| 68 | +printM m = |
| 69 | + let ((r1, c1), (rn, cn)) = bounds m |
| 70 | + in unlines |
| 71 | + [ intercalate "\t" [show $ m ! (r, c) | c <- [c1 .. cn]] |
| 72 | + | r <- [r1 .. rn] |
| 73 | + ] |
| 74 | + |
| 75 | +printV :: (Show a) => Vector a -> String |
| 76 | +printV = unlines . map show . elems |
| 77 | + |
| 78 | +main :: IO () |
| 79 | +main = do |
| 80 | + let mat = [2, 3, 4, 6, 1, 2, 3, 4, 3, -4, 0, 10] :: [Ratio Int] |
| 81 | + m = listArray ((1, 1), (3, 4)) mat |
| 82 | + putStrLn "Original Matrix:" |
| 83 | + putStrLn $ printM m |
| 84 | + putStrLn "Echelon form" |
| 85 | + putStrLn $ printM $ gaussianElimination m |
| 86 | + putStrLn "Reduced echelon form" |
| 87 | + putStrLn $ printM $ gaussJordan $ gaussianElimination m |
| 88 | + putStrLn "Solution from back substitution" |
| 89 | + putStrLn $ printV $ backSubstitution $ gaussianElimination m |
0 commit comments