{-# language ScopedTypeVariables #-} ----------------------------------------------------------------------------- -- | -- Copyright : (C) 2016 Marco Zocca -- License : GPL-3 (see LICENSE) -- Maintainer : zocca.marco gmail -- Stability : provisional -- Portability : portable -- ----------------------------------------------------------------------------- module LibSpec where import Numeric.LinearAlgebra.Sparse -- import Numeric.LinearAlgebra.Class import Control.Monad (liftM, liftM2) import Control.Monad.Primitive import Data.Foldable (foldrM) import Data.Sparse.Common import Control.Monad (replicateM) import Control.Monad.State.Strict (execState) import qualified System.Random.MWC as MWC import qualified System.Random.MWC.Distributions as MWC import Test.Hspec -- import Test.Hspec.QuickCheck main :: IO () main = hspec spec -- niter = 5 spec :: Spec spec = do describe "Numeric.LinearAlgebra.Sparse : library" $ do -- prop "subtraction is cancellative" $ \(x :: SpVector Double) -> -- x ^-^ x `shouldBe` zero it "dot : inner product" $ tv0 `dot` tv0 `shouldBe` 61 it "transposeSM : sparse matrix transpose" $ transposeSM m1 `shouldBe` m1t it "matVec : matrix-vector product" $ nearZero ( normSq ((aa0 #> x0true) ^-^ b0 )) `shouldBe` True it "vecMat : vector-matrix product" $ nearZero ( normSq ((x0true <# aa0) ^-^ aa0tx0 ))`shouldBe` True it "matMat : matrix-matrix product" $ (m1 `matMat` m2) `shouldBe` m1m2 it "eye : identity matrix" $ infoSM (eye 10) `shouldBe` SMInfo 10 0.1 it "insertCol : insert a column in a SpMatrix" $ insertCol (eye 3) (fromListDenseSV 3 [2,2,2]) 0 `shouldBe` (fromListSM (3,3) [(0,0,2),(1,0,2),(1,1,1),(2,0,2),(2,2,1)]) it "insertRow : insert a row in a SpMatrix" $ insertRow (eye 3) (fromListDenseSV 3 [2,2,2]) 1 `shouldBe` (fromListSM (3,3) [(0,0,1), (1,0,2), (1,1,2), (1,2,2), (2,2,1)]) it "extractCol -> insertCol : identity" $ insertCol (eye 3) (extractCol (eye 3) 1) 1 `shouldBe` eye 3 it "extractRow -> insertRow : identity" $ insertRow (eye 3) (extractRow (eye 3) 1) 1 `shouldBe` eye 3 it "countSubdiagonalNZ : # of nonzero elements below the diagonal" $ countSubdiagonalNZSM m3 `shouldBe` 1 it "permutPairsSM : permutation matrices are orthogonal" $ do let pm0 = permutPairsSM 3 [(0,2), (1,2)] :: SpMatrix Double pm0 ##^ pm0 `shouldBe` eye 3 pm0 #^# pm0 `shouldBe` eye 3 it "isLowerTriSM : checks whether matrix is lower triangular" $ isLowerTriSM tm8' && isUpperTriSM tm8 `shouldBe` True it "modifyInspectN : early termination by iteration count" $ execState (modifyInspectN 2 (nearZero . diffSqL) (/2)) (1 :: Double) `shouldBe` 1/8 it "modifyInspectN : termination by value convergence" $ nearZero (execState (modifyInspectN (2^16) (nearZero . head) (/2)) (1 :: Double)) `shouldBe` True describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers" $ do -- it "TFQMR (2 x 2 dense)" $ -- normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True it "BCG (2 x 2 dense)" $ nearZero (normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true)) `shouldBe` True it "BiCGSTAB (2 x 2 dense)" $ nearZero (normSq (aa0 <\> b0 ^-^ x0true)) `shouldBe` True it "CGS (2 x 2 dense)" $ nearZero (normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true)) `shouldBe` True describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ do it "LU (unoptimized) (4 x 4 sparse)" $ checkLuSolve aa1 b1 `shouldBe` True describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do it "QR (4 x 4 sparse)" $ checkQr tm4 `shouldBe` True it "QR (3 x 3 dense)" $ checkQr tm2 `shouldBe` True describe "Numeric.LinearAlgebra.Sparse : LU decomposition" $ do it "LU (4 x 4 dense)" $ checkLu tm6 `shouldBe` True it "LU (10 x 10 sparse)" $ checkLu tm7 `shouldBe` True describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ it "chol (5 x 5 sparse)" $ checkChol tm7 `shouldBe` True describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration" $ do it "Arnoldi iteration (3 x 3 dense)" $ checkArnoldi aa2 3 `shouldBe` True it "Arnoldi iteration (5 x 5 sparse)" $ checkArnoldi tm7 5 `shouldBe` True {- QR-} checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool checkQr a = c1 && c2 where (q, r) = qr a c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a) c2 = isOrthogonalSM q {- LU -} checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool checkLu a = lup == a where (l, u) = lu a lup = l #~# u {- Cholesky -} checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where l = chol a {- direct linear solver -} checkLuSolve :: (Epsilon a, Real a, Floating a) => SpMatrix a -> SpVector a -> Bool checkLuSolve amat rhs = nearZero (normSq ( (lmat #> (umat #> xlu)) ^-^ rhs )) where (lmat, umat) = lu amat xlu = luSolve lmat umat rhs {- Arnoldi iteration -} checkArnoldi :: (Epsilon a, Floating a, Eq a) => SpMatrix a -> Int -> Bool checkArnoldi aa kn = nearZero $ normFrobenius $ (aa #~# qvprev) ^-^ (qv #~# hh) where (qv, hh) = arnoldi aa kn (m, n) = dim qv qvprev = extractSubmatrix qv (0, m - 1) (0, n - 2) {- example 0 : 2x2 linear system [1 2] [2] = [8] [3 4] [3] [18] [1 3] [2] = [11] [2 4] [3] [16] -} aa0 :: SpMatrix Double aa0 = fromListDenseSM 2 [1,3,2,4] -- b0, x0 : r.h.s and initial solution resp. b0, x0, x0true :: SpVector Double b0 = mkSpVectorD 2 [8,18] x0 = mkSpVectorD 2 [0.3,1.4] -- x0true : true solution x0true = mkSpVectorD 2 [2,3] aa0tx0 = mkSpVectorD 2 [11,16] {- 4x4 system -} aa1 :: SpMatrix Double aa1 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12] x1, b1 :: SpVector Double x1 = mkSpVectorD 4 [1,2,3,4] b1 = mkSpVectorD 4 [30,56,60,101] {- 3x3 system -} aa2 :: SpMatrix Double aa2 = sparsifySM $ fromListDenseSM 3 [2, -1, 0, -1, 2, -1, 0, -1, 2] x2, b2 :: SpVector Double x2 = mkSpVectorD 3 [3,2,3] b2 = mkSpVectorD 3 [4,-2,4] aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double -- -- {- example 1 : random linear system -} -- -- dense -- solveRandom n = do -- aa0 <- randMat n -- let aa = aa0 ^+^ eye n -- xtrue <- randVec n -- -- x0 <- randVec n -- let b = aa #> xtrue -- dx = aa <\> b ^-^ xtrue -- return $ normSq dx -- -- let xhatB = _xBicgstab (bicgstab aa b x0 x0) -- -- xhatC = _x (cgs aa b x0 x0) -- -- return (aa, x, x0, b, xhatB, xhatC) -- -- sparse -- solveSpRandom :: Int -> Int -> IO Double -- solveSpRandom n nsp = do -- aa0 <- randSpMat n nsp -- let aa = aa0 ^+^ eye n -- xtrue <- randSpVec n nsp -- let b = (aa ^+^ eye n) #> xtrue -- dx = aa <\> b ^-^ xtrue -- return $ normSq dx -- solveRandomBanded n bw mu sig = do -- let ndiags = 2*bw -- bands <- replicateM (ndiags + 1) (randArray n mu sig) -- xtrue <- randVec n -- b <- randVec n -- let -- diags = [-bw .. bw - 1] randDiagMat :: PrimMonad m => Rows -> Double -> Double -> Int -> m (SpMatrix Double) randDiagMat n mu sig i = do x <- randArray n mu sig return $ mkSubDiagonal n i x go (m:ms) mat = m ^+^ go ms mat go [] mat = mat plusM x y = return $ x ^+^ y -- {- matMat [1, 2] [5, 6] = [19, 22] [3, 4] [7, 8] [43, 50] -} m1 = fromListDenseSM 2 [1,3,2,4] m2 = fromListDenseSM 2 [5, 7, 6, 8] m1m2 = fromListDenseSM 2 [19, 43, 22, 50] -- transposeSM m1t = fromListDenseSM 2 [1,2,3,4] -- {- countSubdiagonalNZ -} m3 = fromListSM (3,3) [(0,2,3),(2,0,4),(1,1,3)] {- mkSubDiagonal -} {- eigenvalues -} aa3 = fromListDenseSM 3 [1,1,3,2,2,2,3,1,1] :: SpMatrix Double b3 = mkSpVectorD 3 [1,1,1] :: SpVector Double -- aa4 : eigenvalues 1 (mult.=2) and -1 aa4 = fromListDenseSM 3 [3,2,-2,2,2,-1,6,5,-4] :: SpMatrix Double b4 = fromListDenseSV 3 [-3,-3,-3] :: SpVector Double -- test data tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double tm0 = fromListSM (2,2) [(0,0,pi), (1,0,sqrt 2), (0,1, exp 1), (1,1,sqrt 5)] tv0, tv1 :: SpVector Double tv0 = mkSpVectorD 2 [5, 6] tv1 = fromListSV 2 [(0,1)] -- wikipedia test matrix for Givens rotation tm1 = sparsifySM $ fromListDenseSM 3 [6,5,0,5,1,4,0,4,3] tm1g1 = givens tm1 1 0 tm1a2 = tm1g1 ## tm1 tm1g2 = givens tm1a2 2 1 tm1a3 = tm1g2 ## tm1a2 tm1q = transposeSM (tm1g2 ## tm1g1) -- wp test matrix for QR decomposition via Givens rotation tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41] tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9] tm3g1 = fromListDenseSM 3 [1, 0,0, 0,c,-s, 0, s, c] where c= 0.4961 s = 0.8682 -- tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12] tm5 = fromListDenseSM 3 [2, -4, -4, -1, 6, -2, -2, 3, 8] :: SpMatrix Double tm6 = fromListDenseSM 4 [1,3,4,2,2,5,2,10,3,6,8,11,4,7,9,12] :: SpMatrix Double tm7 :: SpMatrix Double tm7 = a ^+^ b ^+^ c where n = 5 a = mkSubDiagonal n 1 $ replicate n (-1) b = mkSubDiagonal n 0 $ replicate n 2 c = mkSubDiagonal n (-1) $ replicate n (-1) -- -- run N iterations -- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB -- runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where -- x0 = mkSpVectorD nd $ replicate nd 0.9 -- nd = dim r0 -- r0 = b ^-^ (aa #> x0) -- p0 = r0 -- bicgsInit = BICGSTAB x0 r0 p0 -- -- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS -- runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where -- x0 = mkSpVectorD nd $ replicate nd 0.1 -- nd = dim r0 -- r0 = b ^-^ (aa #> x0) -- residual of initial guess solution -- p0 = r0 -- u0 = r0 -- cgsInit = CGS x0 r0 p0 u0 -- solveRandomN ndim nsp niter = do -- aa0 <- randSpMat ndim (nsp ^ 2) -- let aa = aa0 ^+^ eye ndim -- xtrue <- randSpVec ndim nsp -- let b = aa #> xtrue -- xhatB = head $ runNBiC niter aa b -- xhatC = head $ runNCGS niter aa b -- -- printDenseSM aa -- return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue)) tm8 :: SpMatrix Double tm8 = fromListSM (2,2) [(0,0,1), (0,1,1), (1,1,1)] tm8' :: SpMatrix Double tm8' = fromListSM (2,2) [(0,0,1), (1,0,1), (1,1,1)]