module Main where

import Numeric.LinearAlgebra as LA
import Graphics.EasyPlot

data TriangType = Up | Low | None deriving Show

a n = build (n, n) $ gen n :: Matrix Double

gen n x y | x == (fromIntegral $ n - 1) = 1
          | x == y = 1
          | x < y = -1
          | x > y = 0

luGet a = case lu a of
    (l, u , p, s) -> (p, ((l, Low), (u, Up)))

solveTriang (a, t) f = case t of
    Up -> let rs = toRows a in
        calcu (cols a) (reverse rs) f 
    Low -> let rs = toRows a in
        calcl (cols a) rs f 
    None -> error "Not triang"

calcu nrs rs f = calc' rs 0 [] where 
    calc' [] n us = us
    calc' (r:rs) n us = calc' rs (n + 1) (ui : us) where
        ui = (fi - (subVector 0 n $ vector us) <.> (subVector (nrs - n) n r)) / ri where 
                  fi = f `atIndex` (nrs - n - 1)
                  ri = (r `atIndex` (nrs - n - 1))

calcl nrs rs f = calc' rs 0 [] where 
    calc' [] n us = reverse us
    calc' (r:rs) n us = calc' rs (n + 1) (ui : us) where
        ui = (fi - ((subVector 0 n $ vector $ reverse us) <.> (subVector 0 n r))) / ri where
                  fi = f `atIndex` n
                  ri =  r `atIndex` n

-- Для вычисления числа обусловенности
norm1 = maximum . map (sum . map abs) . toLists
norm2 = norm1 . tr
norm3 m =  sqrt $ maximum $  map magnitude $ toList $ eigenvalues $ m LA.<> (tr m)

mu norm m = norm m * (norm $ inv m)

collectDots norm 1 = [(0, mu norm (a 1))]
collectDots norm n = (fromIntegral n, mu norm (a n)) : collectDots norm (n - 1)

main :: IO ()
main = do
    -- Пример решения системы n = 5
    let m = a 5
    let f = vector $ take 5 $ repeat 1
    let v = vector $ solveTriang (fst $ snd $ luGet m) f
    let u = solveTriang (snd $ snd $ luGet m) v
    --print u :: Vector R

    let dots = collectDots norm1 100
    plot X11 $ Data2D [Title "norm1", Style Lines] [] dots
    print dots