module Main where
import qualified Data.Matrix as M
import qualified Data.Vector as V
import Numeric.Field.Fraction (Fraction)
import Test.Hspec
import Test.Hspec.QuickCheck
import Test.QuickCheck
import Utils
main :: IO ()
main = hspec $ do
describe "luDecomp'" $ do
prop "works correctly" $ prop_luDecomp
return ()
prop_luDecomp :: Property
prop_luDecomp = forAll (resize 10 arbitrary) $ \(MatrixCase ms) ->
let m = M.fromLists ms :: M.Matrix (Fraction Integer)
Just (u,l,p,q,_,_) = M.luDecomp' m
in collect (minimum $ map (minimum . map abs) ms) $ collect (M.ncols m, M.nrows m) $
p * m * q == l * u && isUpperTriangle u && isLowerTriangle l && M.diagProd l == 1
isLowerTriangle :: (Num r, Eq r) => M.Matrix r -> Bool
isLowerTriangle m = all (\i -> all (\j -> m M.! (i, j) == 0) [i+1..M.ncols m]) [1..M.nrows m]
isUpperTriangle :: (Num r, Eq r) => M.Matrix r -> Bool
isUpperTriangle m = all (\i -> all (\j -> m M.! (i, j) == 0) [1..min (i-1) (M.ncols m)]) [1..M.nrows m]
rank :: (Eq r, Num r, Ord r, Fractional r) => M.Matrix r -> Int
rank mat =
let Just (u, _, _, _,_, _) = M.luDecomp' mat
in V.foldr (\a acc -> if a /= 0 then acc + 1 else acc) (0 :: Int) $ M.getDiag u