Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- flatten [] = []
- flatten (x : xs) = x ++ (flatten xs)
- -- ### I/O ### ----------------------------------------------------------------
- -- reads first n lines from stdin and returns them in a list
- getLines :: (Ord t, Num t) => t -> IO [String]
- getLines n
- | n <= 0 = return []
- | otherwise = do
- x <- getLine
- xs <- getLines (n - 1)
- return (x:xs)
- -- transforms a list of lists to a printable matrix string
- vec_to_str [] = "\n"
- vec_to_str (x : xs) = (take 6 ((if x < 0 then "" else " ") ++ (show x) ++ [' ' | _ <- [1..]])) ++ "\t" ++ (vec_to_str xs)
- mat_to_str [] = ""
- mat_to_str (row : rows) = (vec_to_str row) ++ (mat_to_str rows)
- -- given a string representing vector b, finds vector x s.t. Ax = b
- get_solution_string str_b lower upper swaps =
- let b = map read (words str_b) :: [Double]
- y = forward_sub 1 lower [] (apply_swaps b swaps)
- x = backward_sub upper y
- in "Ax = " ++ (show b) ++ " => x = " ++ (show x) ++ "\n"
- -- given a list of strings of vectors, returns a string of solutions
- get_all_solutions strings lower upper swaps =
- flatten [get_solution_string b lower upper swaps | b <- strings]
- -- ######################################################################### --
- -- ### Linear Algebra ### -----------------------------------------------------
- mul_vv [] [] = 0
- mul_vv (x : xs) (y : ys) = x * y + (mul_vv xs ys)
- div_vc [] _ = []
- div_vc (x : xs) c = (x / c) : (div_vc xs c)
- div_mc [] _ = []
- div_mc (row : rows) c = (div_vc row c) : (div_mc rows c)
- outer_product [] _ = []
- outer_product (x : xs) ys = [x * y | y <- ys] : (outer_product xs ys)
- sub_vv [] [] = []
- sub_vv (x : xs) (y : ys) = (x - y) : (sub_vv xs ys)
- sub_mm [] [] = []
- sub_mm (row1 : mat1) (row2 : mat2) = (sub_vv row1 row2) : (sub_mm mat1 mat2)
- -- ######################################################################### --
- -- ### Mics ### ----------------------------------------------------------------
- -- splits the matrix returning the first column and the rest of the matrix
- split [] = ([], [])
- split ((x : xs) : rows) =
- let (col, mat) = split rows
- in ((x : col), (xs : mat))
- -- reverse of split
- merge [] [] = []
- merge (x : xs) (row : mat) = ((x : row) : (merge xs mat))
- -- swaps two elements of list at positions i1, j1
- swap i1 j1 list =
- if i1 == j1 then
- list
- else
- let i = minimum [i1, j1]
- j = maximum [i1, j1]
- first = take i list
- second = drop (i + 1) (take j list)
- third = drop (j + 1) list
- in first ++ [list !! j] ++ second ++ [list !! i] ++ third
- apply_swaps list [] = list
- apply_swaps list ((i, j) : swaps) = apply_swaps (swap i j list) swaps
- -- ######################################################################### --
- -- ### LUP decomposition ### ---------------------------------------------------
- -- returns the index of the maximum element in the list
- max_idx _ [] _ res = res
- max_idx i (x : xs) cur_max res =
- if (abs x) > cur_max then
- max_idx (i + 1) xs (abs x) i
- else
- max_idx (i + 1) xs cur_max res
- -- pivots the matrix bringing the biggest element to the top left corner.
- -- returns the new matrix and the permutation represented as a lsit of swaps.
- pivot mat =
- let (row, rest) = split mat
- i = max_idx 0 row (-1) (-1)
- in ((swap 0 i mat), (0, i))
- add_1 [] = []
- add_1 ((i, j) : rest) = (((i + 1), (j + 1)) : (add_1 rest))
- -- computes the LUP decomposition returning (L, U, P).
- lup_dec [[x]] = ([[1]], [[x]], [])
- lup_dec input_mat =
- let (((a : row) : rows), q) = pivot input_mat
- (col, mat) = split rows
- schur = sub_mm mat (div_mc (outer_product col row) a)
- (l1, u1, p1) = lup_dec schur
- (swap_i, swap_j) = head p1
- l_top = (1 : [0 | _ <- row])
- u_left = (a : [0 | _ <- col])
- l = (l_top : (merge (div_vc (apply_swaps col p1) a) l1))
- u = merge u_left (row : u1)
- p = (q : (add_1 p1))
- in (l, u, p)
- --in ((a : row) : rows, schur, [q])
- -- ######################################################################### --
- -- ### Determinant ### --------------------------------------------------------
- determinant [] = 1
- determinant ((x : row) : upper) =
- let (col, mat) = split upper
- in x * (determinant mat)
- -- ######################################################################### --
- -- ### Linear equations ### ---------------------------------------------------
- forward_sub _ [] xs _ = xs
- forward_sub i (row : lower) xs (b : bs) =
- let prefix = mul_vv [x | (x, j) <- (zip row [1..(i - 1)])] xs
- x = b - prefix
- in forward_sub (i + 1) lower (xs ++ [x]) bs
- backward_sub [[x]] [b] = [b / x]
- backward_sub ((x : row) : upper) (b : bs) =
- let (col, mat) = split upper
- xs = backward_sub mat bs
- prefix = mul_vv row xs
- in ((b - prefix) / x) : xs
- -- ######################################################################### --
- main :: IO()
- main = do
- -- reading and parsing the size of the matrix
- n_s <- getLine
- let n = read n_s :: Integer
- -- reading rows of the matrix as a list of strings
- rows <- getLines n
- -- transforming the list of strings to a double array of Integers
- let mat = [row | row_s <- rows, let row = (map read (words row_s) :: [Double])]
- -- calculating the LU-decomposition
- let (lower, upper, p) = lup_dec mat
- ---- calculating the determinant using the decomposed matrices
- let det = determinant upper * (if ((length p) `div` 2) == 0 then 1 else -1)
- putStrLn $ "Input matrix: \n" ++ (mat_to_str mat) ++ "\n\n"
- putStrLn $ "Unit-Lower Triangular: \n" ++ (mat_to_str lower) ++ "\n"
- putStrLn $ "Upper Triangular: \n" ++ (mat_to_str upper) ++ "\n"
- putStrLn $ "Permutations: " ++ (show p) ++ "\n"
- putStrLn $ "Determinant: " ++ (show det) ++ "\n"
- -- reading the amount of queries for linear equation solutions
- k_s <- getLine
- let k = read k_s :: Integer
- -- reading queries as a list of strings
- queries <- getLines k
- let solutions = get_all_solutions queries lower upper p
- putStrLn $ "\n\nSolutions: \n\n" ++ solutions