Awful attempt to code gauss' method

qn-spghetti.hs

import Data.List
import Test.QuickCheck
import Control.Monad.IO.Class
import Criterion.Main

import Control.Parallel 
import Control.DeepSeq
import Control.Parallel.Strategies

data Matrix = Matrix [[Double]] 
    deriving (Show)

instance Arbitrary Matrix where
    arbitrary = sized matrix'
        where matrix' n = do
                v <- mapM (\_ -> vector (n+2)) [1..(n+1)]
                return $ Matrix $ map (map (+1)) v

gauss :: [[Double]] -> [Double]
gauss m = 
    let ls = buildMatrix [] m -- Строим матрицу упрощенных выражений, будет что-то вроде
                              -- [[a1, b],[a_1,a_2,b1],...,[.....]]
    in foldl (\r (x:xs) -> ( (last xs - sum (zipWith (*) xs r)) / x):r ) [] ls 
        -- решаем, для каждой строки с конца (они уже так упорядочены) нужно вычислить следующий 
        -- коффициент по формуле (b- sum (aI*rI)) / a1
        -- собственно это и делается, сначала берётся пустой массив решений и рассматриваем наш массив 
        -- как:  (x,xs) где x - старший коэффициент (вобще это 1, так что можно было бы и упростить)
        --       xs           - кроме последнего это остальные коэффициенты
        --       tail xs      - колонка решения
        --  zipWith (*) xs r  - находим сумму, важно что lenght xs гарантировано равно length r +1
        --  далее получили новый коэффициент и положили его в список найденный, и идём на следующий ряд
        --  в итоге получаем правильно отсортированный список решений
        --  NB. мне тут не нравится last и возможно его как-то можно убрать, но думать было лень
    where
        -- тут мы строим матрицу в правильном виде
        buildMatrix ms [] = ms  -- мы обработали все ряды
        buildMatrix ms m  =     -- обрабатываем кусок матрицы
            let (maxV, maxRow, other) = foldl' 
                            (\(mX,r,ls) xm@(x:xs)  -> if x > mX then (x,xm,r:ls) else (mX,r,xm:ls)) 
                                (head $ head m,head m,[]) 
                                (tail m) 
                -- нахождение максимального старшего коэффициента, ряда его содержащего и списка остальных
                -- сначала предполагаем, что это первый ряд (head $ head m {- 1 коэф 1-го ряда-}, head m, [])
                -- потом прогоняем по списку все ряды. если первый коэффициент другого ряда больше текущего
                -- максимального, то кладём это значение как максимум, запонимаем ряд, старый ряд кладём в список
                -- иначе просто кладём тестируемый ряд в список
                newR = map (/maxV) maxRow -- все коэффициенты максимального ряда делим на 1 коэффициент
            in buildMatrix (newR:ms) 
                        (map (\xs@(x:_) ->tail $ map (\(v',r') -> v'-x*r')  $ zip xs newR) other) 
                -- пераходим на следующую итерацию:
                -- 1). добавляем максимальный ряд в список обрадоботанных, 
                -- 2). остальные же надо обновить:
                --  нужно из каждого коэффициента ряда вычесть 1 коэф* тот же член максимального ряда:
                --  соотв map (...) other делает это
                --      к каждому ряду \xs@(x:_) {-посмотрели 1 элемент-} применяем вычитание
                --      для того, чтобы коэффициенты 1 ряда и остальных были "рядом" и не нужен был доступ по
                --      индексу делаем zip xs newR итого map (\(v',r') -> v' - x * r') обрабатывает ряд.
                --  Теперь первый коэффициент у нас заведомо 0, поэтому чтобы не тратить память и исключить
                --  ошибки арифметики с плавающей точкой выкивыдаем его tail $
                --  Итого мы получили список всех строк из которых вычеркнули первый коэффициент.

-- 1 we should evaluate newR before next step
gauss1 :: [[Double]] -> [Double]
gauss1 m = 
    let ls = buildMatrix [] m 
    in foldl (\r (x:xs) -> ( (last xs - sum (zipWith (*) xs r)) / x):r ) [] ls 
    where
        buildMatrix ms [] = ms
        buildMatrix ms m  =
            let (maxV, maxRow, other) = foldl' 
                            (\(mX,r,ls) xm@(x:xs)  -> if x > mX then (x,xm,r:ls) else (mX,r,xm:ls)) 
                                (head $ head m,head m,[]) 
                                (tail m) 
                newR = force $ map (/maxV) maxRow
            in  buildMatrix (newR:ms) 
                        ((map (\xs@(x:_) ->tail $ map (\(v',r') -> v'-x*r')  $ zip xs newR) other) `using` (parList rpar))

main = do
    (Matrix x) <- testT 100
    defaultMain
       [ bench "default" $ nf gauss  x
       , bench "1"       $ nf gauss1 x
       ]

testT n = do
    (a:_) <- sample' (resize n  (arbitrary::Gen Matrix))
    return a


--g_prop :: [[Double]] -> Bool
g_prop :: Matrix -> Bool
g_prop (Matrix m) = 
    let r = gauss m
    in all (\x -> sum (zipWith (*) r x) ~= last x) m


test (Matrix m) = 
    let r = gauss m
    in putStrLn $ unlines $ map (show . (\x -> abs $ sum (zipWith (*) r x) - last x)) m

(~=) :: Double -> Double -> Bool
x ~= y = abs (x - y) < 1e-8 
infixl 4 ~=

testM :: [[Double]]
testM = [ [1,2,3,4]
        , [0,1,0,0]
        ]
testM1 :: [[Double]]
testM1 = [ [3,2,7,14]
         , [8,3,19,22]
         , [3,4,22,18]
         ]

Spaghetti.hs

gauss  :: [[Double]] -> [Double]
gauss  m = let (m',is) = gauss0 m
           in order is $ gauss2 $ gauss1 m

gauss0 :: [[Double]] -> ([[Double]],[Int])
gauss0 m = let (     mt,is) = foldr choose (trns m,straight) straight
           in  (trns mt,is)
    where choose k  (mt,is) = let row = m !! (k-1)
                                  q   = (+) 1 $ fromJust $ elemIndex (maximum $ init $ drop (k-1) row) row
                              in (swap' k q mt, swap' k q is)
          straight = [1..sizeM m]

gauss1 :: [[Double]] ->  [[Double]]
gauss1 m = foldr eliminate m [1..sizeM m]
    where eliminate k m' = let oldRow = m'     !! (k-1)
                               mkk    = oldRow !! (k-1)
                               newRow = map (/mkk) oldRow
                               sub curRow = let mik = curRow !! (k-1)
                                            in  map (\(j,mij) -> mij - (newRow !! (j-1))*mik) $ zip [1..] curRow
                           in  map (\(i,row) -> if i < k+1
                                                then if i == k
                                                     then newRow
                                                     else row
                                                else sub  row) $ zip [1..] m'

gauss2 :: [[Double]] -> [Double]
gauss2   = gauss2' . ([],) . reverse . (map $ dropWhile ((<0.05) . abs))
    where  gauss2' (r,[   ]  ) = r
           gauss2' (r,[x,f]:m) = let r' = f/x
                               in gauss2' (r':r, map (reverse . (\(g:x':t) -> (g - x'*r'):t) . reverse) m)
           gauss2' z = unsafePerformIO (print z) `seq` undefined

--swap'  :: Int -> Int -> [a] -> [a]
swap'  :: Show a => Int -> Int -> [a] -> [a]
swap' i j xs = if   i == j
               then xs
               else
                   let     (l,a:t) = splitAt (i'   -1) xs
                   in  let (m,b:r) = splitAt (j'-i'-1) t
                       in   l ++ (b:m) ++ (a:r)
    where i' = min i j
          j' = max i j

order  :: [Int] -> [a] -> [a]
order is xs = let m = fromList $ zip [1..] xs
              in  map (m!) is