sparse-linear-algebra 0.2.2.0 → 0.2.9
raw patch · 21 files changed
+3788/−1754 lines, 21 filesdep +exceptionsdep +transformersdep +vector-algorithmsdep ~QuickCheckdep ~base
Dependencies added: exceptions, transformers, vector-algorithms, vector-space, writer-cps-transformers
Dependency ranges changed: QuickCheck, base
Files
- README.md +66/−59
- app/Main.hs +0/−4
- sparse-linear-algebra.cabal +35/−19
- src/Control/Exception/Common.hs +78/−0
- src/Control/Iterative.hs +333/−0
- src/Data/Sparse/Common.hs +139/−50
- src/Data/Sparse/Internal/CSR.hs +0/−34
- src/Data/Sparse/Internal/IntM.hs +127/−0
- src/Data/Sparse/Internal/IntMap2.hs +58/−80
- src/Data/Sparse/Internal/Utils.hs +301/−0
- src/Data/Sparse/PPrint.hs +58/−0
- src/Data/Sparse/SpMatrix.hs +200/−85
- src/Data/Sparse/SpVector.hs +138/−61
- src/Data/Sparse/Types.hs +6/−0
- src/Data/Sparse/Utils.hs +28/−0
- src/Numeric/Eps.hs +25/−9
- src/Numeric/LinearAlgebra/Class.hs +352/−112
- src/Numeric/LinearAlgebra/EigenSolvers/Experimental.hs +51/−0
- src/Numeric/LinearAlgebra/LinearSolvers/Experimental.hs +59/−0
- src/Numeric/LinearAlgebra/Sparse.hs +1017/−1049
- test/LibSpec.hs +717/−192
README.md view
@@ -2,13 +2,13 @@ Numerical computation in native Haskell -TravisCI : [](https://travis-ci.org/ocramz/sparse-linear-algebra)+[](https://hackage.haskell.org/package/sparse-linear-algebra) [](https://travis-ci.org/ocramz/sparse-linear-algebra) This library provides common numerical analysis functionality, without requiring any external bindings. It is not optimized for performance (yet), but it serves as an experimental platform for scientific computation in a purely functional setting. Contents : -* Iterative linear solvers (`linSolve`)+* Iterative linear solvers (`<\>`) * Generalized Minimal Residual (GMRES) (non-Hermitian systems) @@ -18,12 +18,12 @@ * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems) - * Transpose-Free Quasi-Minimal Residual (TFQMR)+ * Moore-Penrose pseudoinverse (`pinv`) (rectangular systems) * Direct linear solvers - * LU-based (`luSolve`)-+ * LU-based (`luSolve`); forward and backward substitution (`triLowerSolve`, `triUpperSolve`)+ * Matrix factorization algorithms * QR (`qr`)@@ -32,10 +32,10 @@ * Cholesky (`chol`) -* Eigenvalue algorithms- * Arnoldi iteration (`arnoldi`) +* Eigenvalue algorithms+ * QR (`eigsQR`) * Rayleigh quotient iteration (`eigRayleigh`)@@ -45,6 +45,19 @@ * Predicates : Matrix orthogonality test (A^T A ~= I) ++### Under development++* Matrix factorization algorithms++ * Golub-Kahan-Lanczos bidiagonalization+ + * Singular value decomposition (SVD)++* Iterative linear solvers++ * Transpose-Free Quasi-Minimal Residual (TFQMR)+ --------- ## Examples@@ -53,123 +66,116 @@ ### Creation of sparse data -The `fromListSM` function creates a sparse matrix from an array of its entries we use :+The `fromListSM` function creates a sparse matrix from a collection of its entries in (row, column, value) format. This is its type signature: fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a -e.g.+and, in case you have a running GHCi session (the terminal is denoted from now on by `λ>`), you can try something like this: - > amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)]+ λ> amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)] :: SpMatrix Double -and similarly+Similarly, `fromListSV` is used to create sparse vectors: fromListSV :: Int -> [(Int, a)] -> SpVector a+ -can be used to create sparse vectors. Alternatively, the user can copy the contents of a list to a (dense) SpVector using fromListDenseSV :: Int -> [a] -> SpVector a + ### Displaying sparse data Both sparse vectors and matrices can be pretty-printed using `prd`: - > prd amat+ λ> prd amat ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 ) - [2,0,0]- [4,3,2]- [0,0,5]+ 2.0 _ _ + 4.0 3.0 2.0+ _ _ 5.0 -The zeros are just added at printing time; sparse vectors and matrices should only contain non-zero entries.+Note: sparse data should only contain non-zero entries not to waste memory and computation. ### Matrix operations -There are a few common matrix factorizations available; in the following example we compute the LU factorization of a matrix and verify it with the matrix-matrix product `##` :+There are a few common matrix factorizations available; in the following example we compute the LU factorization of matrix `amat` and verify it with the matrix-matrix product `##` of its factors : - > (l, u) = lu amat- > prd $ l ## u+ λ> (l, u) <- lu amat+ λ> prd $ l ## u ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 ) - [2.0,0.0,0.0]- [4.0,3.0,2.0]- [0.0,0.0,5.0]+ 2.0 _ _ + 4.0 3.0 2.0+ _ _ 5.0 Notice that the result is _dense_, i.e. certain entries are numerically zero but have been inserted into the result along with all the others (thus taking up memory!). To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined in `Numeric.Eps`) depends on the numerical type used (e.g. it is 10^-6 for `Float`s and 10^-12 for `Double`s). - > prd $ l #~# u+ λ> prd $ l #~# u ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 ) - [2.0,0.0,0.0]- [4.0,3.0,2.0]- [0.0,0.0,5.0]+ 2.0 _ _ + 4.0 3.0 2.0+ _ _ 5.0 -A matrix is transposed using `transposeSM`.+A matrix is transposed using the `transpose` function. -Sometimes we need to compute matrix-matrix transpose products, which is why the library offers the infix operators `#^#` (M^T N) and `##^` (M N^T):+Sometimes we need to compute matrix-matrix transpose products, which is why the library offers the infix operators `#^#` (i.e. matrix transpose * matrix) and `##^` (matrix * matrix transpose): - > amat' = amat #^# amat- > prd amat'+ λ> amat' = amat #^# amat+ λ> prd amat' ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 ) - [20.0,12.0,8.0]- [12.0,9.0,6.0]- [8.0,6.0,29.0]-- > l = chol amat'- > prd $ l ##^ l+ 20.0 12.0 8.0+ 12.0 9.0 6.0+ 8.0 6.0 29.0+ + λ> l <- chol amat'+ λ> prd $ l ##^ l ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 ) - [20.000000000000004,12.0,8.0]- [12.0,9.0,10.8]- [8.0,10.8,29.0]+ 20.000000000000004 12.0 8.0+ 12.0 9.0 10.8+ 8.0 10.8 29.0 -In the above example we have also shown the Cholesky decomposition (M = L L^T where L is a lower-triangular matrix), which is only possible for symmetric positive-definite matrices.+In the last example we have also shown the Cholesky decomposition (M = L L^T where L is a lower-triangular matrix), which is only defined for symmetric positive-definite matrices. ### Linear systems -Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `linSolve` function, or alternatively `<\>` (which uses GMRES as default solver method) :+Large sparse linear systems are best solved with iterative methods. `sparse-linear-algebra` provides a selection of these via the `<\>` (inspired by Matlab's "backslash" function. Here we use GMRES as default solver method) : - > b = fromListDenseSV 3 [3,2,5]- > x = amat <\> b- > prd x+ λ> b = fromListDenseSV 3 [3,2,5] :: SpVector Double+ λ> x <- amat <\> b+ λ> prd x ( 3 elements ) , 3 NZ ( sparsity 1.0 ) - [1.4999999999999998,-1.9999999999999998,0.9999999999999998]+ 1.4999999999999998 -1.9999999999999998 0.9999999999999998 The result can be verified by computing the matrix-vector action `amat #> x`, which should (ideally) be very close to the right-hand side `b` : - > prd $ amat #> x+ λ> prd $ amat #> x ( 3 elements ) , 3 NZ ( sparsity 1.0 ) - [2.9999999999999996,1.9999999999999996,4.999999999999999]+ 2.9999999999999996 1.9999999999999996 4.999999999999999 The library also provides a forward-backward substitution solver (`luSolve`) based on a triangular factorization of the system matrix (usually LU). This should be the preferred for solving smaller, dense systems. Using the data defined above we can cross-verify the two solution methods: - > x' = luSolve l u b- > prd x'+ λ> x' <- luSolve l u b+ λ> prd x' ( 3 elements ) , 3 NZ ( sparsity 1.0 ) - [1.5,-2.0,1.0]+ 1.5 -2.0 1.0 ----------- -This is also an experiment in principled scientific programming : -* set the stage by declaring typeclasses and some useful generic operations (normed linear vector spaces, i.e. finite-dimensional spaces equipped with an inner product that induces a distance function), -* define appropriate data structures, and how they relate to those properties (sparse vectors and matrices, defined internally via `Data.IntMap`, are made instances of the VectorSpace and Additive classes respectively). This allows to decouple the algorithms from the actual implementation of the backend,--* implement the algorithms, following 1:1 the textbook [1, 2] -- ## License GPL3, see LICENSE@@ -179,6 +185,7 @@ Inspired by * `linear` : https://hackage.haskell.org/package/linear+* `vector-space` : https://hackage.haskell.org/package/vector-space * `sparse-lin-alg` : https://github.com/laughedelic/sparse-lin-alg ## References
− app/Main.hs
@@ -1,4 +0,0 @@-module Main where--main :: IO ()-main = putStrLn $ "2 + 3 = " ++ show (2 + 3)
sparse-linear-algebra.cabal view
@@ -1,6 +1,6 @@ name: sparse-linear-algebra-version: 0.2.2.0-synopsis: Numerical computation in native Haskell +version: 0.2.9+synopsis: Sparse linear algebra algorithms and datastructures for scientific computation, in native Haskell. Iterative linear solvers, matrix factorizations, linear eigensolvers and related utilities. description: /Overview/ .@@ -30,30 +30,43 @@ hs-source-dirs: src exposed-modules: Numeric.LinearAlgebra.Sparse Numeric.LinearAlgebra.Class+ Numeric.Eps Data.Sparse.SpVector Data.Sparse.SpMatrix Data.Sparse.Common+ Control.Exception.Common other-modules: Data.Sparse.Internal.IntMap2- Data.Sparse.Internal.CSR+ Data.Sparse.Internal.IntM+ Data.Sparse.Internal.Utils Data.Sparse.Utils+ Data.Sparse.PPrint Data.Sparse.Types- Numeric.Eps- build-depends: QuickCheck- , base >= 4.7 && < 5+ Control.Iterative+ Numeric.LinearAlgebra.LinearSolvers.Experimental+ Numeric.LinearAlgebra.EigenSolvers.Experimental+ build-depends: base >= 4.7 && < 5+ , vector-space , containers- , hspec- , primitive >= 0.6.1.0- , mtl >= 2.2.1- , mwc-random , vector+ , vector-algorithms+ , exceptions+ , mtl >= 2.2.1 + , transformers+ , writer-cps-transformers+ -- , monad-log == 0.1.1.0+ -- , monad-par+ , QuickCheck >= 2.8.2+ , hspec+ -- , primitive >= 0.6.1.0+ -- , mwc-random -executable sparse-linear-algebra- default-language: Haskell2010- ghc-options: -threaded -rtsopts -with-rtsopts=-N- hs-source-dirs: app- main-is: Main.hs- build-depends: base+-- executable sparse-linear-algebra+-- default-language: Haskell2010+-- ghc-options: -threaded -rtsopts -with-rtsopts=-N+-- hs-source-dirs: app+-- main-is: Main.hs+-- build-depends: base test-suite spec@@ -63,16 +76,19 @@ hs-source-dirs: test other-modules: LibSpec main-is: Spec.hs- build-depends: base+ build-depends: QuickCheck >= 2.8.2+ , base , containers+ , criterion , hspec , mtl >= 2.2.1+ , exceptions , mwc-random , primitive >= 0.6.1.0 , sparse-linear-algebra- , criterion- -- , QuickCheck+ , vector-space source-repository head type: git location: https://github.com/ocramz/sparse-linear-algebra+
+ src/Control/Exception/Common.hs view
@@ -0,0 +1,78 @@+module Control.Exception.Common where++-- import Control.Exception.Safe (Exception, MonadThrow, SomeException, throwM)+import Control.Exception+import Control.Monad.Catch (MonadThrow (..))+import Data.Typeable -- (TypeRep, Typeable, typeRep)++import Data.Sparse.Utils++++data PartialFunctionError = EmptyList String deriving (Eq, Typeable)+instance Show PartialFunctionError where+ show (EmptyList s) = unwords [s, ": empty list"]+instance Exception PartialFunctionError+++-- | Input error++data InputError = NonNegError String Int deriving (Eq, Typeable)+instance Show InputError where+ show (NonNegError s i) = unwords [s, ": parameter must be nonnegative, instead I got", show i]+instance Exception InputError++++-- | Out of bounds index error+data OutOfBoundsIndexError i = OOBIxError String i+ | OOBIxsError String [i]+ | OOBNoCompatRows String (i,i) deriving (Eq, Typeable)+instance Show i => Show (OutOfBoundsIndexError i) where+ show (OOBIxError e i) = unwords [e, ": index", show i,"out of bounds"]+ show (OOBIxsError e ixs) = unwords [e, ":, indices", show ixs, "out of bounds"]+ show (OOBNoCompatRows e ij) = unwords [e, ": no compatible rows for indices", show ij]+instance (Show i, Typeable i) => Exception (OutOfBoundsIndexError i)++checkIxBound :: MonadThrow m => String -> Int -> UB -> m a -> m a+checkIxBound e i n ff+ | not (inBounds0 n i) = throwM (OOBIxError e i)+ | otherwise = ff++ ++-- | Operand size mismatch errors+data OperandSizeMismatch = DotSizeMismatch Int Int+ | NonTriangularException String+ | MatVecSizeMismatchException String (Int, Int) Int deriving (Eq, Typeable)+instance Show OperandSizeMismatch where+ show (DotSizeMismatch na nb) = unwords ["<.> : Incompatible dimensions : ", show na, show nb]+ show (NonTriangularException s )= unwords [s, ": Matrix must be triangular"]+ show (MatVecSizeMismatchException s dm o) = unwords [s, ": Matrix-vector dimensions are incompatible: Matrix is",show dm,", whereas vector is",show o]+instance Exception OperandSizeMismatch+++++-- | Matrix exceptions+data MatrixException i = HugeConditionNumber String i+ | NeedsPivoting String String deriving (Eq, Typeable)+instance Show i => Show (MatrixException i) where+ show (HugeConditionNumber s x) = unwords [s, ": Rank-deficient system: condition number", show x]+ show (NeedsPivoting s1 s2) = unwords [s1, ":", s2, "is close to 0. Permute the rows to obtain a nonzero diagonal"]+instance (Show i, Typeable i) => Exception (MatrixException i)+++++++-- | Numerical iteration errors+data IterationException a = NotConvergedE String Int a + | DivergingE String Int a a+ deriving (Eq, Typeable)+instance Show a => Show (IterationException a) where+ show (NotConvergedE s niters x) = unwords [s, ": Could not converge within",show niters, "iterations; final state:", show x]+ show (DivergingE s niters x0 x1) = unwords [s, ": Diverging at iteration #", show niters, "; latest state:",show x0, "; previous state:",show x1]+ -- show (NotConvergedToExpectedE s niters x x0) = unwords [s, ": Could not converge within",show niters, "iterations; final state:", show x, "expected state:", show x0]+instance (Show a, Typeable a) => Exception (IterationException a)
+ src/Control/Iterative.hs view
@@ -0,0 +1,333 @@+{-# language FlexibleContexts, GeneralizedNewtypeDeriving, DeriveFunctor #-}+-----------------------------------------------------------------------------+-- |+-- Module : Control.Iterative+-- Copyright : (c) Marco Zocca 2017+-- License : GPL-style (see the file LICENSE)+--+-- Maintainer : zocca marco gmail+-- Stability : experimental+-- Portability : portable+--+-- Combinators and helper functions for iterative algorithms, with support for monitoring and exceptions.+--+-----------------------------------------------------------------------------+module Control.Iterative where++import Control.Exception.Common+import Numeric.LinearAlgebra.Class+import Numeric.Eps++import Control.Monad.Catch+import Data.Typeable++import Control.Monad (when)+-- import Control.Monad.Trans.Reader+import Control.Monad.State.Strict+-- import Control.Monad.Trans.Writer.CPS+import Control.Monad.Trans.Class (lift)+import qualified Control.Monad.Trans.State.Strict as MTS -- (runStateT)+import Data.Foldable (foldrM)++import Data.VectorSpace+++data ConvergenceStatus a = BufferNotReady+ | Converging+ | Converged a+ | Diverging a a+ | NotConverged+ deriving (Eq, Show)++ +data IterationConfig a b =+ IterConf { numIterationsMax :: Int,+ printDebugInfo :: Bool,+ iterationView :: a -> b, + printDebugIO :: b -> IO ()}+instance Show (IterationConfig a b) where+ show (IterConf n qd _ _) = unwords ["Max. # of iterations:",show n,", print debug information:", show qd]+ +++-- * Control primitives for bounded iteration with convergence check++-- -- | transform state until a condition is met+modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s+modifyUntil q f = modifyUntilM q (pure . f)+ +modifyUntilM :: MonadState s m => (s -> Bool) -> (s -> m s) -> m s+modifyUntilM q f = do+ x <- get+ y <- f x+ put y+ if q y then return y+ else modifyUntilM q f ++++++++-- | `untilConvergedG0` is a special case of `untilConvergedG` that assesses convergence based on the L2 distance to a known solution `xKnown`+untilConvergedG0 ::+ (Normed v, MonadThrow m, MonadIO m, Typeable (Magnitude v), Typeable s, Show s) =>+ String+ -> IterationConfig s v+ -> v -- ^ Known value+ -> (s -> s)+ -> s+ -> m s+untilConvergedG0 fname config xKnown f x0 = + modifyInspectGuarded fname config norm2Diff nearZero qdiverg qfin f x0+ where+ qfin s = nearZero $ norm2 (xKnown ^-^ s)+ +++-- | This function makes some default choices on the `modifyInspectGuarded` machinery: convergence is assessed using the squared L2 distance between consecutive states, and divergence is detected when this function is increasing between pairs of measurements.+untilConvergedG :: (Normed v, MonadThrow m, MonadIO m, Typeable (Magnitude v), Typeable s, Show s) =>+ String+ -> IterationConfig s v+ -> (v -> Bool)+ -> (s -> s) + -> s + -> m s+untilConvergedG fname config =+ modifyInspectGuarded fname config norm2Diff nearZero qdiverg+++-- | ", monadic version+untilConvergedGM ::+ (Normed v, MonadThrow m, MonadIO m, Typeable (Magnitude v), Typeable s, Show s) =>+ String+ -> IterationConfig s v+ -> (v -> Bool)+ -> (s -> m s) + -> s+ -> m s+untilConvergedGM fname config =+ modifyInspectGuardedM fname config norm2Diff nearZero qdiverg++++++-- | `modifyInspectGuarded` is a high-order abstraction of a numerical iterative process. It accumulates a rolling window of 3 states and compares a summary `q` of the latest 2 with that of the previous two in order to assess divergence (e.g. if `q latest2 > q prev2` then it). The process ends when either we hit an iteration budget or relative convergence is verified. The function then assesses the final state with a predicate `qfinal` (e.g. against a known solution; if this is not known, the user can just supply `const True`)+modifyInspectGuarded ::+ (MonadThrow m, MonadIO m, Typeable s, Typeable a, Show s, Show a) =>+ String -- ^ Calling function name+ -> IterationConfig s v -- ^ Configuration+ -> ([v] -> a) -- ^ State summary array projection+ -> (a -> Bool) -- ^ Convergence criterion+ -> (a -> a -> Bool) -- ^ Divergence criterion+ -> (v -> Bool) -- ^ Final state acceptance criterion+ -> (s -> s) -- ^ State evolution+ -> s -- ^ Initial state+ -> m s -- ^ Final state+modifyInspectGuarded fname config sf qc qd qfin f x0 =+ modifyInspectGuardedM fname config sf qc qd qfin (pure . f) x0++ +++-- | ", monadic version+modifyInspectGuardedM ::+ (MonadThrow m, MonadIO m, Typeable s, Show s, Typeable a, Show a) =>+ String+ -> IterationConfig s v+ -> ([v] -> a)+ -> (a -> Bool)+ -> (a -> a -> Bool)+ -> (v -> Bool)+ -> (s -> m s)+ -> s+ -> m s+modifyInspectGuardedM fname config sf qconverg qdiverg qfinal f x0 + | nitermax > 0 = MTS.execStateT (go 0 []) x0+ | otherwise = throwM (NonNegError fname nitermax)+ where+ lwindow = 3+ nitermax = numIterationsMax config+ pf = iterationView config+ checkConvergStatus y i ll+ | length ll < lwindow = BufferNotReady+ | qdiverg qi qt && not (qconverg qi) = Diverging qi qt + | qconverg qi || qfinal (pf y) = Converged qi+ | i == nitermax - 1 = NotConverged + | otherwise = Converging+ where llf = pf <$> ll+ qi = sf $ init llf -- summary of latest 2 states+ qt = sf $ tail llf -- " " previous 2 states+ go i ll = do+ x <- MTS.get+ y <- lift $ f x+ when (printDebugInfo config) $ liftIO $ do+ putStrLn $ unwords ["Iteration", show i]+ printDebugIO config (pf y) + case checkConvergStatus y i ll of+ BufferNotReady -> do + MTS.put y+ let ll' = y : ll -- cons current state to buffer+ go (i + 1) ll'+ Converged qi -> MTS.put y+ Diverging qi qt -> do+ MTS.put y+ throwM (DivergingE fname i qi qt)+ Converging -> do+ MTS.put y+ let ll' = init (y : ll) -- rolling state window+ go (i + 1) ll'+ NotConverged -> do+ MTS.put y+ throwM (NotConvergedE fname nitermax y)+ +++-- | Some useful combinators+++-- | Apply a function over a range of integer indices, zip the result with it and filter out the almost-zero entries+onRangeSparse :: Epsilon b => (Int -> b) -> [Int] -> [(Int, b)]+onRangeSparse f ixs = foldr ins [] ixs where+ ins x xr | isNz (f x) = (x, f x) : xr+ | otherwise = xr++-- | ", monadic version+onRangeSparseM :: (Epsilon b, Foldable t, Monad m) =>+ (a -> m b) -> t a -> m [(a, b)]+onRangeSparseM f ixs = unfoldZipM mf f ixs where+ mf x = isNz <$> f x+ +++unfoldZipM0 :: (Foldable t, Monad m) =>+ (a -> Bool) -> (a -> b) -> t a -> m [(a, b)]+unfoldZipM0 q f = unfoldZipM (pure . q) (pure . f)+++unfoldZipM :: (Foldable t, Monad m) =>+ (a -> m Bool) -> (a -> m b) -> t a -> m [(a, b)]+unfoldZipM q f ixs = foldrM insf [] ixs where+ insf x xr = do+ qx <- q x+ if qx+ then do+ y <- f x+ pure $ (x, y) : xr+ else pure xr++-- | A combinator I don't know how to call+combx :: Functor f => (a -> b) -> (t -> f a) -> t -> f b+combx g f x = g <$> f x +++ + ++++-- | Helpers+++-- meanl :: (Foldable t, Fractional a) => t a -> a+-- meanl xx = 1/fromIntegral (length xx) * sum xx++-- norm2l :: (Foldable t, Functor t, Floating a) => t a -> a+-- norm2l xx = sqrt $ sum (fmap (**2) xx)++-- | Squared difference of a 2-element list.+-- | NB: unsafe !+diffSqL :: Floating a => [a] -> a+diffSqL xx = (x1 - x2)**2 where [x1, x2] = [head xx, xx!!1]+++-- | Relative tolerance :+-- relTol a b := ||a - b|| / (1 + min (||norm2 a||, ||norm2 b||))+relTol :: Normed v => v -> v -> Magnitude v+relTol a b = norm2 (a ^-^ b) / m where+ m = 1 + min (norm2 a) (norm2 b)+++qdiverg :: Ord a => a -> a -> Bool+qdiverg = (>)++norm2Diff [s1, s0] = norm2 (s1 ^-^ s0)+norm2Diff _ = 1/0+++++++-- test data++data S = S {unS1 :: Double, unS2 :: String} deriving (Eq, Show)+liftS1 f (S x i) = S (f x) i+s0 = S 1 "blah"+ic1 = IterConf 2 True unS1 print++++++-- playground++-- instance MonadThrow m => MonadThrow (WriterT w m) where+-- throwM = lift . throwM++-- -- | iter0 also accepts a configuration, e.g. for optional printing of debug info+-- -- iter0 :: MonadIO m =>+-- -- Int -> (s -> m s) -> (s -> String) -> IterationConfig s -> s -> m s+-- iter0 nmax f sf config x0 = flip runReaderT config $ MTS.execStateT (go (0 :: Int)) x0+-- where+-- go i = do+-- x <- get+-- c <- lift $ asks printDebugInfo -- neat+-- y <- lift . lift $ f x -- not neat+-- when c $ liftIO $ putStrLn $ sf y +-- put y+-- unless (i >= nmax) (go $ i + 1)+ ++-- -- | iter1 prints output at every iteration until the loop terminates OR is interrupted by an exception, whichever happens first+-- -- iter1 :: (MonadThrow m, MonadIO m, Typeable t, Show t) =>+-- -- (t -> m t) -> (t -> String) -> (t -> Bool) -> (t -> Bool) -> t -> m t+-- iter1 f wf qe qx x0 = execStateT (go 0) x0 where+-- go i = do+-- x <- get+-- y <- lift $ f x+-- _ <- liftIO $ wf y+-- when (qx y) $ throwM (NotConvergedE "bla" (i+1) y)+-- put y+-- unless (qe y) $ go (i + 1) +++-- -- | iter2 concatenates output with WriterT but does NOT `tell` any output if an exception is raised before the end of the loop+-- iter2 :: (MonadThrow m, Monoid w, Typeable t, Show t) => (t -> m t)+-- -> (t -> w) -> (t -> Bool) -> (t -> Bool) -> t -> m (t, w)+-- iter2 f wf qe qx x0 = runWriterT $ execStateT (go 0) x0 where+-- go i = do+-- x <- get+-- y <- lift . lift $ f x+-- lift $ tell $ wf y+-- when (qx y) $ throwM (NotConvergedE "bla" (i+1) y)+-- put y+-- unless (qe y) $ go (i + 1) +++-- -- test :: IO (Int, [String])+-- test :: IO Int+-- -- test :: IO ()+-- test = do+-- (yt, w ) <- iter2 f wf qe qexc x0+-- putStrLn w+-- return yt+-- -- iter1 f wf qe qexc x0+-- where+-- f = pure . (+ 1)+-- wf v = unwords ["state =", show v]+-- qe = (== 5)+-- qexc = (== 3)+-- x0 = 0 :: Int
src/Data/Sparse/Common.hs view
@@ -1,3 +1,6 @@+{-# language TypeFamilies #-}+{-# language TypeOperators, GADTs #-}+{-# language FlexibleInstances, FlexibleContexts #-} ----------------------------------------------------------------------------- -- | -- Copyright : (C) 2016 Marco Zocca@@ -18,34 +21,48 @@ extractDiagDense, extractSubRow, extractSubCol, extractSubRow_RK, extractSubCol_RK,- matVec, (#>), vecMat, (<#),- fromCols,- prd) where+ fromCols, fromColsL, toCols) where +-- import Control.Exception+-- import Control.Exception.Common+-- import Control.Monad.Catch+ import Data.Sparse.Utils as X+import Data.Sparse.PPrint as X import Data.Sparse.Types as X import Data.Sparse.Internal.IntMap2 -- as X+import qualified Data.Sparse.Internal.IntM as I+import Data.Sparse.Internal.IntM (IntM(..)) import Data.Sparse.SpMatrix as X import Data.Sparse.SpVector as X+-- import Data.Sparse.Internal.CSR as X import Numeric.Eps as X import Numeric.LinearAlgebra.Class as X import qualified Data.IntMap as IM +import GHC.Exts+import Data.Complex++-- import Control.Applicative+-- import Data.Traversable+ import Data.Maybe (fromMaybe, maybe) import qualified Data.Vector as V +import Data.VectorSpace+ -- withBoundsSM m ij e f -- | isValidIxSM m ij = f m ij -- | otherwise = error e --- | modify the size of a SpVector. Do not use directly+-- | Modify the size of a SpVector. Do not use directly resizeSV :: Int -> SpVector a -> SpVector a resizeSV d2 (SV _ sv) = SV d2 sv---mapKeysSV fk (SV d sv) = SV d $ IM.mapKeys fk sv+-- | Remap the keys of a SpVector. Do not use directly+mapKeysSV :: (IM.Key -> IM.Key) -> SpVector a -> SpVector a+mapKeysSV fk (SV d sv) = SV d $ I.mapKeys fk sv -- * Insert row/column vector in matrix @@ -53,10 +70,10 @@ insertRowWith :: (IxCol -> IxCol) -> SpMatrix a -> SpVector a -> IM.Key -> SpMatrix a insertRowWith fj (SM (m,n) im) (SV d sv) i | not (inBounds0 m i) = error "insertRowSM : index out of bounds"- | n >= d = SM (m,n) $ IM.insert i (insertOrUnion i sv' im) im+ | n >= d = SM (m,n) $ I.insert i (insertOrUnion i sv' im) im | otherwise = error $ "insertRowSM : incompatible dimensions " ++ show (n, d)- where sv' = IM.mapKeys fj sv- insertOrUnion i sv im = maybe sv (IM.union sv) (IM.lookup i im)+ where sv' = I.mapKeys fj sv+ insertOrUnion i' sv' im' = maybe sv' (I.union sv') (I.lookup i' im') -- | Insert row insertRow :: SpMatrix a -> SpVector a -> IM.Key -> SpMatrix a@@ -71,7 +88,7 @@ (m, n) = dim smm mv = dim sv vl = toListSV sv- insIM2 im2 ((i,x):xs) j = insIM2 (insertSpMatrix (fi i) j x im2) xs j+ insIM2 im2 ((i,x):xs) j' = insIM2 (insertSpMatrix (fi i) j' x im2) xs j' insIM2 im2 [] _ = im2 -- | Insert column@@ -107,24 +124,26 @@ -- | promote a SV to SM svToSM :: SpVector a -> SpMatrix a-svToSM (SV n d) = SM (n, 1) $ IM.singleton 0 d+svToSM (SV n d) = SM (n, 1) $ I.singleton 0 d -- |Demote (n x 1) or (1 x n) SpMatrix to SpVector toSV :: SpMatrix a -> SpVector a-toSV (SM (m,n) im) = SV d (ff im) where- ff | m > n = IM.map g -- column case- | otherwise = g- g = snd . head . IM.toList+toSV (SM (m, n) im) = SV d im' where+ im' | m < n = snd . head . toList $ im+ | otherwise = fmap g im+ g = snd . head . toList d | m==1 && n==1 = 1 | m==1 && n>1 = n | n==1 && m>1 = m- | otherwise = error $ "toSV : incompatible matrix dimension " ++ show (m,n)+ | otherwise = error $ "toSV : incompatible matrix dimension " ++ show (m,n) ++ -- | Lookup a row in a SpMatrix; returns an SpVector with the row, if this is non-empty lookupRowSM :: SpMatrix a -> IxRow -> Maybe (SpVector a)-lookupRowSM sm i = SV (ncols sm) <$> IM.lookup i (dat sm) +lookupRowSM sm i = SV (ncols sm) <$> I.lookup i (dat sm) -- * Extract a SpVector from an SpMatrix@@ -148,7 +167,7 @@ Num a => (Int -> (IxRow, IxCol)) -> SpMatrix a -> SpVector a extractVectorDenseWith f mm = fromListDenseSV n $ foldr ins [] ll where ll = [0 .. n - 1]- (m, n) = dim mm+ (_, n) = dim mm ins i acc = mm @@ f i : acc -- | Extract ith row (dense)@@ -211,29 +230,87 @@ :: (Normed f1, Num b, Functor f) => f (f1 b) -> f1 b -> f b -} +instance LinearVectorSpace (SpVector Double) where+ type MatrixType (SpVector Double) = SpMatrix Double+ (#>) = matVecSD+ (<#) = vecMatSD +instance LinearVectorSpace (SpVector (Complex Double)) where+ type MatrixType (SpVector (Complex Double)) = SpMatrix (Complex Double)+ (#>) = matVecSD+ (<#) = vecMatSD --- |Matrix-on-vector-matVec, (#>) :: Num a => SpMatrix a -> SpVector a -> SpVector a-matVec (SM (nr, nc) mdata) (SV n sv)- | nc == n = SV nr $ fmap (`dot` sv) mdata- | otherwise = error $ "matVec : mismatching dimensions " ++ show (nc, n)+ -(#>) = matVec+matVecSD :: InnerSpace (IntM t) =>+ SpMatrix t -> SpVector t -> SpVector (Scalar (IntM t))+matVecSD (SM (nr, nc) mdata) (SV n sv)+ | nc == n = SV nr $ fmap (`dot` sv) mdata+ | otherwise = error $ "matVec : mismatched dimensions " ++ show (nc, n) -- |Vector-on-matrix (FIXME : transposes matrix: more costly than `matVec`, I think)-vecMat, (<#) :: Num a => SpVector a -> SpMatrix a -> SpVector a -vecMat (SV n sv) (SM (nr, nc) mdata)+vecMatSD :: InnerSpace (IntM t) =>+ SpVector t -> SpMatrix t -> SpVector (Scalar (IntM t))+vecMatSD (SV n sv) (SM (nr, nc) mdata) | n == nr = SV nc $ fmap (`dot` sv) (transposeIM2 mdata) | otherwise = error $ "vecMat : mismatching dimensions " ++ show (n, nr) -(<#) = vecMat +-- -- generalized matVec : we require a function `rowsf` that produces a functor of elements of a Hilbert space (the rows of `m`)+-- matVecG :: (Hilbert v, Functor f, f (Scalar v) ~ v) => (m -> f v) -> m -> v -> v+-- matVecG rowsf m v = fmap (`dot` v) (rowsf m) +-- matVecGA+-- :: (Hilbert v, Traversable t, t (Scalar v) ~ v) =>+-- (m -> t v) -> m -> v -> v+-- matVecGA rowsf m v = traverse (<.> v) (rowsf m) +-- -- -- Really, a matrix is just notation for a linear map between two finite-dimensional Hilbert spaces, i.e.+-- matVec :: (Hilbert u, Hilbert v) => (u -> v) -> u -> v+-- which is a specialization of a function application operator like ($) :: (a -> b) -> a -> b+++++-- -- -- from `vector-space`++-- data a -* b where+-- Dot :: VectorSpace b => b -> (b -* Scalar b)+-- (:&&) :: (a -* c) -> (a -* d) -> (a -* (c, d)) -- a,c,d should be constrained++-- apply :: Hilbert a => (a -* b) -> (a -> b)+-- apply (Dot b) = dot b+-- apply (f :&& g) = apply f &&& apply g+-- where (u &&& v) a = (u a, v a) -- (&&&) from Control.Arrow++-- -- type a :~ b = Scalar a ~ Scalar b+++++++++-- | Pack a list of SpVectors into an SpMatrix+fromColsL :: [SpVector a] -> SpMatrix a+fromColsL = fromCols . V.fromList++-- | Unpack an SpMatrix into a list of SpVectors+toCols :: SpMatrix a -> [SpVector a]+toCols aa = map (extractCol aa) [0 .. n-1] where+ (m,n) = dim aa++++++++ -- | Pack a V.Vector of SpVectors as columns of an SpMatrix fromCols :: V.Vector (SpVector a) -> SpMatrix a fromCols qv = V.ifoldl' ins (zeroSM m n) qv where@@ -244,29 +321,36 @@ +++++ -- * Pretty printing -showNonZero :: (Show a, Num a, Eq a) => a -> String-showNonZero x = if x == 0 then " " else show x+showNz :: (Epsilon a, Show a) => a -> String+showNz x | nearZero x = " _ "+ | otherwise = show x +-- toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]+toDenseRow sm irow =+ fmap (\icol -> sm @@ (irow,icol)) [0..ncol-1] where (_, ncol) = dim sm -toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a]-toDenseRow (SM (_,ncol) im) irow =- fmap (\icol -> im `lookupWD_IM` (irow,icol)) [0..ncol-1] -toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String++-- toDenseRowClip :: (Show a, Num a) => SpMatrix a -> IM.Key -> Int -> String toDenseRowClip sm irow ncomax- | ncols sm > ncomax = unwords (map show h) ++ " ... " ++ show t- | otherwise = show dr+ | nco > ncomax = unwords (map showNz h) ++ " ... " ++ showNz t+ | otherwise = unwords $ showNz <$> dr where dr = toDenseRow sm irow h = take (ncomax - 2) dr t = last dr+ (_, nco) = dim sm -newline :: IO ()-newline = putStrLn "" -printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()+-- printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()+printDenseSM :: (ScIx c ~ (Int, Int), FDSize c ~ (Int, Int), SpContainer c a, Show a, Epsilon a) => c a -> IO () printDenseSM sm = do newline putStrLn $ sizeStr sm@@ -274,22 +358,24 @@ printDenseSM' sm 5 5 newline where - printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()- printDenseSM' sm'@(SM (nr,_) _) nromax ncomax = mapM_ putStrLn rr_' where+ -- printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()+ printDenseSM' sm' nromax ncomax = mapM_ putStrLn rr_' where+ (nr, _) = dim sm' rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]- rr_' | nrows sm > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]+ rr_' | nr > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_] | otherwise = rr_ -toDenseListClip :: (Show a, Num a) => SpVector a -> Int -> String++toDenseListClip :: (Show a, Epsilon a) => SpVector a -> Int -> String toDenseListClip sv ncomax- | dim sv > ncomax = unwords (map show h) ++ " ... " ++ show t- | otherwise = show dr+ | dim sv > ncomax = unwords (map showNz h) ++ " ... " ++ showNz t+ | otherwise = unwords $ showNz <$> dr where dr = toDenseListSV sv h = take (ncomax - 2) dr t = last dr -printDenseSV :: (Show t, Num t) => SpVector t -> IO ()+printDenseSV :: (Show t, Epsilon t) => SpVector t -> IO () printDenseSV sv = do newline putStrLn $ sizeStrSV sv@@ -302,19 +388,22 @@ | otherwise = rr_ -- ** Pretty printer typeclass-class PrintDense a where- prd :: a -> IO () -instance (Show a, Num a) => PrintDense (SpVector a) where++instance (Show a, Num a, Epsilon a) => PrintDense (SpVector a) where prd = printDenseSV -instance (Show a, Num a) => PrintDense (SpMatrix a) where+instance (Show a, Num a, Epsilon a) => PrintDense (SpMatrix a) where prd = printDenseSM +++-- instance (Elt a, Show a) => PrintDense (CsrMatrix a) where+-- prd = printDenseSM
− src/Data/Sparse/Internal/CSR.hs
@@ -1,34 +0,0 @@-module Data.Sparse.Internal.CSR where--import qualified Data.Vector as V-import qualified Data.Vector.Unboxed as VU-import qualified Data.Vector.Unboxed.Mutable as VM--{-| Compressed Row Storage specification :-- http://netlib.org/utk/people/JackDongarra/etemplates/node373.html-- The compressed row storage (CRS) format puts the subsequent nonzeros of the matrix- rows in contiguous memory locations. Assuming we have a nonsymmetric sparse matrix- $A$, we create three vectors: one for floating point numbers (val) and the other- two for integers (col_ind, row_ptr).-- The val vector stores the values of the nonzero elements of the matrix $A$ as- they are traversed in a row-wise fashion.- - The col_ind vector stores the column indexes of the elements in the val vector,- that is, if val(k)=a_{i,j}, then col_ind(k)=j$.-- The row_ptr vector stores the locations in the val vector that start a row;- that is, if val(k)=a_{i,j}, then row_ptr(i) <= k < row_ptr(i+1)---}--data CsrMatrix a =- CsrMatrix { csrVal :: VU.Vector a,- csrColInd :: VU.Vector Int,- csrRowPtr :: VU.Vector Int,- csrNnz :: {-# UNPACK #-} !Int,- csrNrows :: {-# UNPACK #-} !Int,- csrNcols :: {-# UNPACK #-} !Int } deriving Eq-
+ src/Data/Sparse/Internal/IntM.hs view
@@ -0,0 +1,127 @@+{-# language GeneralizedNewtypeDeriving, DeriveFunctor, DeriveFoldable, CPP, TypeFamilies, FlexibleInstances #-}+module Data.Sparse.Internal.IntM where++import Data.Sparse.Utils+import Numeric.LinearAlgebra.Class++import GHC.Exts+import Data.Complex+import Data.VectorSpace+import qualified Data.IntMap as IM++++++-- | A synonym for IntMap +newtype IntM a = IntM {unIM :: IM.IntMap a} deriving (Eq, Show, Functor, Foldable)+++empty = IntM IM.empty++size (IntM x) = IM.size x++singleton i x = IntM $ IM.singleton i x++filterWithKey f im = IntM $ IM.filterWithKey f (unIM im)++insert k x (IntM im) = IntM $ IM.insert k x im++filterI f (IntM im) = IntM $ IM.filter f im++lookup x (IntM im) = IM.lookup x im++lookupLT x (IntM im) = IM.lookupLT x im++foldlWithKey f z (IntM im) = IM.foldlWithKey f z im++foldlWithKey' f z (IntM im) = IM.foldlWithKey' f z im++mapWithKey f (IntM im) = IntM $ IM.mapWithKey f im++keys (IntM im) = IM.keys im++mapKeys f (IntM im) = IntM $ IM.mapKeys f im++union (IntM a) (IntM b) = IntM $ IM.union a b++findMin (IntM im) = IM.findMin im+findMax (IntM im) = IM.findMax im++(IntM im) ! i = im IM.! i+++instance IsList (IntM a) where+ type Item (IntM a) = (Int, a)+ fromList = IntM . IM.fromList+ toList = IM.toList . unIM++++instance Set IntM where+ liftU2 f (IntM a) (IntM b) = IntM $ IM.unionWith f a b+ liftI2 f (IntM a) (IntM b) = IntM $ IM.intersectionWith f a b++instance Num a => AdditiveGroup (IntM a) where+ zeroV = IntM IM.empty+ {-# INLINE zeroV #-}+ (^+^) = liftU2 (+)+ {-# INLINE (^+^) #-}+ (^-^) = liftU2 (-)+ {-# INLINE (^-^) #-}+ negateV = fmap negate+ {-# INLINE negateV #-}+++++-- | ParamInstance can be used with all types that are instances of Set (which are by construction also instances of Functor)+#define ParamInstance(f, t) \+ instance VectorSpace (f t) where {type (Scalar (f (t))) = (t); n *^ im = fmap (* n) im};\+ instance VectorSpace (f (Complex t)) where {type (Scalar (f (Complex t))) = Complex (t); n *^ im = fmap (* n) im};\+ instance InnerSpace (f t) where {a <.> b = sum $ liftI2 (*) a b};\+ instance InnerSpace (f (Complex t)) where {a <.> b = sum $ liftI2 (*) (conjugate <$> a) b};\+ -- instance Normed (f t) where {type RealScalar (f t) = t ; type Magnitude (f t) = t ; norm1 a = sum (abs <$> a) ; norm2Sq a = sum $ liftI2 (*) a a; normP p v = sum u**(1/p) where u = fmap (**p) v; normalize = normzPR ; normalize2 = normz2R}; \+ -- instance Normed (f (Complex t)) where {type RealScalar (f (Complex t)) = t; type Magnitude (f (Complex t)) = t; norm1 a = realPart $ sum (abs <$> a); norm2Sq a = realPart $ sum $ liftI2 (*) (conjugate <$> a) a; normP p v = realPart $ sum u**(1/(p :+ 0)) where u = fmap (**(p :+ 0)) v; normalize = normzPC; normalize2 = normz2C }+++instance Normed (IntM Double) where+ type RealScalar (IntM Double) = Double+ type Magnitude (IntM Double) = Double+ norm1 a = sum (abs <$> a)+ norm2Sq a = sum $ liftI2 (*) a a+ normP p v = sum u**(1/p) where u = fmap (**p) v+ normalize p v = v ./ normP p v + normalize2 v = v ./ norm2 v + +instance Normed (IntM (Complex Double)) where+ type RealScalar (IntM (Complex Double)) = Double+ type Magnitude (IntM (Complex Double)) = Double+ norm1 a = realPart $ sum (abs <$> a)+ norm2Sq a = realPart $ sum $ liftI2 (*) (conjugate <$> a) a+ normP p v = realPart $ sum u**(1/(p :+ 0)) where u = fmap (**(p :+ 0)) v+ normalize p v = v ./ toC (normP p v)+ normalize2 v = v ./ toC (norm2 v)+++++++-- | IntMap instances+#define IntMapInstance(t) \+ ParamInstance( IntM, t )++IntMapInstance(Double)+-- IntMapInstance(Float)++++++-- -- | list to IntMap+-- mkIm :: [Double] -> IM.IntMap Double+mkIm xs = fromList $ denseIxArray xs :: IntM Double++-- mkImC :: [Complex Double] -> IM.IntMap (Complex Double)+mkImC xs = fromList $ denseIxArray xs :: IntM (Complex Double)
src/Data/Sparse/Internal/IntMap2.hs view
@@ -1,43 +1,14 @@+{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-} module Data.Sparse.Internal.IntMap2 where -import Numeric.LinearAlgebra.Class+import qualified Data.Sparse.Internal.IntM as I +-- import Numeric.LinearAlgebra.Class import qualified Data.IntMap.Strict as IM import Data.Sparse.Types import Data.Maybe----instance Set IM.IntMap where- liftU2 = IM.unionWith- {-# INLINE liftU2 #-}- liftI2 = IM.intersectionWith- {-# INLINE liftI2 #-}--instance Additive IM.IntMap where- zero = IM.empty- {-# INLINE zero #-}- (^+^) = liftU2 (+)- {-# INLINE (^+^) #-}---instance VectorSpace IM.IntMap where- n .* im = IM.map (* n) im- -instance Hilbert IM.IntMap where- a `dot` b = sum $ liftI2 (*) a b- --instance Normed IM.IntMap where- norm p v | p==1 = norm1 v- | p==2 = norm2 v- | otherwise = normP p v------ set-like brackets-+import GHC.Exts @@ -47,18 +18,20 @@ -- * Insertion -- | Insert an element-insertIM2 ::- IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)-insertIM2 i j x imm = IM.insert i ro imm where- ro = maybe (IM.singleton j x) (IM.insert j x) (IM.lookup i imm)+-- insertIM2 ::+-- IM.Key -> IM.Key -> a -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+insertIM2+ :: IM.Key -> IM.Key -> a -> I.IntM (I.IntM a) -> I.IntM (I.IntM a)+insertIM2 i j x imm = I.insert i ro imm where+ ro = maybe (I.singleton j x) (I.insert j x) (I.lookup i imm) {-# inline insertIM2 #-} -- * Lookup -- |Lookup a key-lookupIM2 ::- IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a-lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j+-- lookupIM2 ::+-- IM.Key -> IM.Key -> IM.IntMap (IM.IntMap a) -> Maybe a+lookupIM2 i j imm = I.lookup i imm >>= I.lookup j {-# inline lookupIM2 #-} -- | Lookup with default 0@@ -69,9 +42,9 @@ -- |Populate an IM2 from a list of (row index, column index, value) -fromListIM2 ::- Foldable t =>- t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+-- fromListIM2 ::+-- Foldable t =>+-- t (IM.Key, IM.Key, a) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) fromListIM2 iix sm = foldl ins sm iix where ins t (i,j,x) = insertIM2 i j x t @@ -79,32 +52,32 @@ -- * Folding -- |Indexed left fold over an IM2, with general accumulator-ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b-ifoldlIM2' f empty mm = IM.foldlWithKey' accRow empty mm where- accRow acc i r = IM.foldlWithKey' (accElem i) acc r+-- ifoldlIM2' :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b+ifoldlIM2' f empty mm = I.foldlWithKey' accRow empty mm where+ accRow acc i r = I.foldlWithKey' (accElem i) acc r accElem i acc j x = f i j x acc {-# inline ifoldlIM2' #-} -- |Indexed left fold over an IM2-ifoldlIM2 ::- (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) ->- IM.IntMap (IM.IntMap t) -> - IM.IntMap a-ifoldlIM2 f m = IM.foldlWithKey' accRow IM.empty m where- accRow acc i row = IM.foldlWithKey' (accElem i) acc row+-- ifoldlIM2 ::+-- (IM.Key -> IM.Key -> t -> IM.IntMap a -> IM.IntMap a) ->+-- IM.IntMap (IM.IntMap t) -> +-- IM.IntMap a+ifoldlIM2 f m = I.foldlWithKey' accRow I.empty m where+ accRow acc i row = I.foldlWithKey' (accElem i) acc row accElem i acc j x = f i j x acc {-# inline ifoldlIM2 #-} -- |Left fold over an IM2, with general accumulator-foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b-foldlIM2 f empty mm = IM.foldl accRow empty mm where- accRow acc r = IM.foldl accElem acc r+-- foldlIM2 :: (a -> b -> b) -> b -> IM.IntMap (IM.IntMap a) -> b+foldlIM2 f empty mm = foldl accRow empty mm where+ accRow acc r = foldl accElem acc r accElem acc x = f x acc {-# inline foldlIM2 #-} -- | Inner indices become outer ones and vice versa. No loss of information because both inner and outer IntMaps are nubbed.-transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+-- transposeIM2 :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) transposeIM2 = ifoldlIM2 (flip insertIM2) {-# inline transposeIM2 #-} @@ -120,25 +93,25 @@ -- * Filtering -- |Map over outer IM and filter all inner IM's-ifilterIM2 ::- (IM.Key -> IM.Key -> a -> Bool) ->- IM.IntMap (IM.IntMap a) ->- IM.IntMap (IM.IntMap a)+-- ifilterIM2 ::+-- (IM.Key -> IM.Key -> a -> Bool) ->+-- IM.IntMap (IM.IntMap a) ->+-- IM.IntMap (IM.IntMap a) ifilterIM2 f =- IM.mapWithKey (\irow row -> IM.filterWithKey (f irow) row) + I.mapWithKey (\irow row -> I.filterWithKey (f irow) row) {-# inline ifilterIM2 #-} -- |Specialized filtering : keep only sub-diagonal elements-filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+-- filterSubdiag :: IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a) filterSubdiag = ifilterIM2 (\i j _ -> i>j) -countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int+-- countSubdiagonalNZ :: IM.IntMap (IM.IntMap a) -> Int countSubdiagonalNZ im =- IM.size $ IM.filter (not . IM.null) (filterSubdiag im)+ I.size $ I.filterI (not . null) (filterSubdiag im) -- |List of (row, col) indices of (nonzero) subdiagonal elements-subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]-subdiagIndices im = concatMap rpairs $ IM.toList (IM.map IM.keys im') where+-- subdiagIndices :: IM.IntMap (IM.IntMap a) -> [(IM.Key, IM.Key)]+subdiagIndices im = concatMap rpairs $ toList (I.keys <$> im') where im' = filterSubdiag im rpairs :: (a, [b]) -> [(a, b)]@@ -157,37 +130,42 @@ -- * Mapping -- |Map over IM2-mapIM2 :: (a -> b) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap b)-mapIM2 = IM.map . IM.map -- imapIM2 (\_ _ x -> f x)+mapIM2 :: (a -> b) -> I.IntM (I.IntM a) -> I.IntM (I.IntM b)+mapIM2 = fmap . fmap -- |Indexed map over IM2-imapIM2 ::- (IM.Key -> IM.Key -> a -> b) ->- IM.IntMap (IM.IntMap a) ->- IM.IntMap (IM.IntMap b)-imapIM2 f im = IM.mapWithKey ff im where- ff j x = IM.mapWithKey (`f` j) x+-- imapIM2 ::+-- (IM.Key -> IM.Key -> a -> b) ->+-- IM.IntMap (IM.IntMap a) ->+-- IM.IntMap (IM.IntMap b)+imapIM2 f im = I.mapWithKey ff im where+ ff j x = I.mapWithKey (`f` j) x -- |Mapping keys-mapKeysIM2 ::- (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)-mapKeysIM2 fi fj im = IM.map adjCols adjRows where- adjRows = IM.mapKeys fi im- adjCols = IM.mapKeys fj +-- mapKeysIM2 ::+-- (IM.Key -> IM.Key) -> (IM.Key -> IM.Key) -> IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+mapKeysIM2 fi fj im = adjCols <$> adjRows where+ adjRows = I.mapKeys fi im+ adjCols = I.mapKeys fj -- map over a single `column` -mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b)+-- mapColumnIM2 :: (b -> b) -> IM.IntMap (IM.IntMap b) -> Int -> IM.IntMap (IM.IntMap b) mapColumnIM2 f im jj = imapIM2 (\i j x -> if j == jj then f x else x) im+++++-- | utilities
+ src/Data/Sparse/Internal/Utils.hs view
@@ -0,0 +1,301 @@+{-# language ExistentialQuantification, TypeFamilies, FlexibleInstances, MultiParamTypeClasses #-}+module Data.Sparse.Internal.Utils where++import Control.Monad (unless)+import Control.Monad.State+import Control.Monad.ST++import Data.Ord (comparing)++import qualified Data.Vector as V +-- import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Mutable as VM+import qualified Data.Vector.Algorithms.Merge as VA++import Data.Ix+import Data.Maybe++import Data.Sparse.Types+import Numeric.LinearAlgebra.Class++++++-- | Given a number of rows(resp. columns) `n` and a _sorted_ Vector of Integers in increasing order (containing the row(col) indices of nonzero entries), return the cumulative vector of nonzero entries of length `n + 1` (the "row(col) pointer" of the CSR(CSC) format). NB: Fused count-and-accumulate+-- E.g.:+-- > csrPtrV (==) 4 (V.fromList [1,1,2,3])+-- [0,0,2,3,4]+csrPtrV :: (a -> Int -> Bool) -> Int -> V.Vector a -> V.Vector Int+csrPtrV eqf n xs = V.create createf where+ createf :: ST s (VM.MVector s Int)+ createf = do+ let c = 0+ vm <- VM.new (n + 1)+ VM.write vm 0 0 -- write `0` at position 0+ let loop v ll i count | i == n = return ()+ | otherwise = do+ let lp = V.length $ V.takeWhile (`eqf` i) ll+ count' = count + lp+ VM.write v (i + 1) count'+ loop v (V.drop lp ll) (succ i) count'+ loop vm xs 0 c+ return vm+++-- csrPtrV' eqf n xs = V.create createf where+-- createf :: ST s (VM.MVector s Int)+-- createf = do+-- let c = 0+-- vm <- VM.new (n + 1)+-- VM.write vm 0 0+-- numLoop fw 0 where+-- fw ix = let lp = V.length $ V.takeWhile (`eqf` ix) ll+-- count+-- return vm++numLoop :: Monad m => (Int -> m a) -> Int -> m ()+numLoop fm n = go 0 where+ go i | i == n = return ()+ | otherwise = do+ _ <- fm i+ go (succ i)++numLoopST' :: Monad m => (Int -> s -> m a) -> Int -> (s -> s) -> s -> m ()+numLoopST' fm n fs s0 = go 0 s0 where+ go i s | i == n = return ()+ | otherwise = do+ _ <- fm i s+ go (succ i) (fs s)++numLoopST'' ::+ MonadState s m => (Int -> s -> m a) -> Int -> (s -> s) -> m ()+numLoopST'' fm n fs = go 0 where+ go i = unless (i == n) $ do+ s <- get+ _ <- fm i s -- ignore result of `fm`+ put $ fs s+ go (succ i)++ + +++++ +-- | O(N) : Intersection between sorted vectors, in-place updates+intersectWith :: Ord o =>+ (a -> o) -> (a -> a -> c) -> V.Vector a -> V.Vector a -> V.Vector c+intersectWith f = intersectWithCompare (comparing f)+++intersectWithCompare ::+ (a1 -> a2 -> Ordering) ->+ (a1 -> a2 -> a) ->+ V.Vector a1 ->+ V.Vector a2 ->+ V.Vector a+intersectWithCompare fcomp g u_ v_ = V.create $ do+ let n = min (V.length u_) (V.length v_)+ vm <- VM.new n+ let go u_ v_ i vm | V.null u_ || V.null v_ || i == n = return (vm, i)+ | otherwise = do+ let (u,us) = (V.head u_, V.tail u_)+ (v,vs) = (V.head v_, V.tail v_)+ case fcomp u v of EQ -> do VM.write vm i (g u v)+ go us vs (i + 1) vm+ LT -> go us v_ i vm+ GT -> go u_ vs i vm+ (vm', i') <- go u_ v_ 0 vm+ let vm'' = VM.take i' vm'+ return vm''++++++++-- | O(N) : Union between sorted vectors, in-place updates+unionWith :: Ord o =>+ (t -> o) -> (t -> t -> a) -> t -> V.Vector t -> V.Vector t -> V.Vector a+unionWith f = unionWithCompare (comparing f)+++unionWithCompare ::+ (t -> t -> Ordering) -> (t -> t -> a) -> t -> V.Vector t -> V.Vector t -> V.Vector a+unionWithCompare fcomp g z u_ v_ = V.create $ do+ let n = (V.length u_) + (V.length v_)+ vm <- VM.new n+ let go u_ v_ i vm+ | (V.null u_ && V.null v_) || i==n = return (vm, i)+ | V.null u_ = do+ VM.write vm i (g z (V.head v_))+ go u_ (V.tail v_) (i+1) vm+ | V.null v_ = do+ VM.write vm i (g (V.head u_) z)+ go (V.tail u_) v_ (i+1) vm+ | otherwise = do+ let (u,us) = (V.head u_, V.tail u_)+ (v,vs) = (V.head v_, V.tail v_)+ case fcomp u v of EQ -> do VM.write vm i (g u v)+ go us vs (i + 1) vm+ LT -> do VM.write vm i (g u z)+ go us v_ (i + 1) vm+ GT -> do VM.write vm i (g z v)+ go u_ vs (i + 1) vm+ (vm', nfin) <- go u_ v_ 0 vm+ let vm'' = VM.take nfin vm'+ return vm''++++++-- * Sorting+sortWith :: Ord o => (t -> o) -> V.Vector t -> V.Vector t+sortWith f v = V.modify (VA.sortBy (comparing f)) v++sortWith3 :: Ord o =>+ ((a, b, c) -> o) ->+ V.Vector a ->+ V.Vector b ->+ V.Vector c ->+ V.Vector (a, b, c)+sortWith3 f x y z = sortWith f $ V.zip3 x y z+++sortByFst3 :: Ord a => V.Vector a -> V.Vector b -> V.Vector c -> V.Vector (a, b, c)+sortByFst3 = sortWith3 fst3++sortBySnd3 :: Ord b => V.Vector a -> V.Vector b -> V.Vector c -> V.Vector (a, b, c)+sortBySnd3 = sortWith3 snd3+++++++++++++-- * Utilities++-- ** 3-tuples+fst3 :: (a, b, c) -> a+fst3 (a, _, _) = a+snd3 :: (a, b, c) -> b+snd3 (_, b, _) = b+third3 :: (a, b, c) -> c+third3 (_, _, c) = c+++tail3 :: (t, t1, t2) -> (t1, t2)+tail3 (_,j,x) = (j,x)++++mapFst3 :: (a -> b) -> (a, y, z) -> (b, y, z)+mapFst3 f (a, b, c) = (f a, b, c)++mapSnd3 :: (a -> b) -> (x, a, z) -> (x, b, z)+mapSnd3 f (a, b, c) = (a, f b, c)++mapThird3 :: (a -> b) -> (x, y, a) -> (x, y, b)+mapThird3 f (a, b, c) = (a, b, f c)++++lift2 :: (a -> b) -> (b -> b -> c) -> a -> a -> c+lift2 p f t1 t2 = f (p t1) (p t2)+++++++-- | Stream fusion based version of the above, from [1]+-- [1] : https://www.schoolofhaskell.com/user/edwardk/revisiting-matrix-multiplication/part-3++data Stream m a = forall s . Stream (s -> m (Step s a)) s ++data Step s a = Yield a s | Skip s | Done ++data MergeState sa sb i a+ = MergeL sa sb i a+ | MergeR sa sb i a+ | MergeLeftEnded sb+ | MergeRightEnded sa+ | MergeStart sa sb++mergeStreamsWith :: (Ord i, Monad m) =>+ (a -> a -> Maybe a) -> Stream m (i, a) -> Stream m (i, a) -> Stream m (i, a)+mergeStreamsWith f (Stream stepa sa0) (Stream stepb sb0)+ = Stream step (MergeStart sa0 sb0) where+ step (MergeStart sa sb) = do+ r <- stepa sa+ return $ case r of+ Yield (i, a) sa' -> Skip (MergeL sa' sb i a)+ Skip sa' -> Skip (MergeStart sa' sb)+ Done -> Skip (MergeLeftEnded sb)+ step (MergeL sa sb i a) = do+ r <- stepb sb+ return $ case r of+ Yield (j, b) sb' -> case compare i j of+ LT -> Yield (i, a) (MergeR sa sb' j b)+ EQ -> case f a b of+ Just c -> Yield (i, c) (MergeStart sa sb')+ Nothing -> Skip (MergeStart sa sb')+ GT -> Yield (j, b) (MergeL sa sb' i a)+ Skip sb' -> Skip (MergeL sa sb' i a)+ Done -> Yield (i, a) (MergeRightEnded sa)+ step (MergeR sa sb j b) = do+ r <- stepa sa+ return $ case r of+ Yield (i, a) sa' -> case compare i j of+ LT -> Yield (i, a) (MergeR sa' sb j b)+ EQ -> case f a b of+ Just c -> Yield (i, c) (MergeStart sa' sb)+ Nothing -> Skip (MergeStart sa' sb)+ GT -> Yield (j, b) (MergeL sa' sb i a)+ Skip sa' -> Skip (MergeR sa' sb j b)+ Done -> Yield (j, b) (MergeLeftEnded sb)+ step (MergeLeftEnded sb) = do+ r <- stepb sb+ return $ case r of+ Yield (j, b) sb' -> Yield (j, b) (MergeLeftEnded sb')+ Skip sb' -> Skip (MergeLeftEnded sb')+ Done -> Done+ step (MergeRightEnded sa) = do+ r <- stepa sa+ return $ case r of+ Yield (i, a) sa' -> Yield (i, a) (MergeRightEnded sa')+ Skip sa' -> Skip (MergeRightEnded sa')+ Done -> Done+ {-# INLINE [0] step #-}+{-# INLINE [1] mergeStreamsWith #-} +++-- test data++-- m0 = V.fromList [O (0,0,1), O(0,1,2)]+-- m1 = V.fromList [O (0,0,1), O(0,2,3)]+++isOrderedV :: Ord a => V.Vector a -> Bool+isOrderedV l = V.all (== True) $ V.zipWith (<) l (V.tail l)+++++++++++
+ src/Data/Sparse/PPrint.hs view
@@ -0,0 +1,58 @@+-----------------------------------------------------------------------------+-- |+-- Copyright : (C) 2016 Marco Zocca+-- License : GPL-3 (see LICENSE)+-- Maintainer : zocca.marco gmail+-- Stability : provisional+-- Portability : portable+--+-----------------------------------------------------------------------------+module Data.Sparse.PPrint where++import Data.Complex+import Text.Printf++import Numeric.Eps++++class PrintDense a where+ prd :: a -> IO ()++newline :: IO ()+newline = putStrLn "" +++++printfDouble :: (PrintfArg t, PrintfType t1) => PPrintOptions -> t -> t1+printfDouble opts x = printf pstr x where+ pstr = concat ["%" , show ni, ".", show nd, "f"]+ nd = pprintDecimals opts+ ni = pprintDigits opts - nd++data PPrintOptions =+ PPrintOptions {+ pprintDigits :: Int,+ pprintDecimals :: Int } deriving (Eq, Show)++pprintDefaults = PPrintOptions 5 2+++++-- | Cleaner way to display Complex values++newtype C a = C {unC :: Complex a} deriving Eq+instance (Num a, Epsilon a, Ord a, Show a) => Show (C a) where+ show (C (r :+ i)) = unwords [show r, oi] where+ oi | isNz i = unwords [s, show i' ++ "j"]+ | otherwise = []+ s | signum i >= 0 = "+"+ | otherwise = "-"+ i' = abs i++-- c0, c1 :: C Double+-- c0 = C $ 1 :+ (-2) +-- c1 = C $ pi :+ 0+
src/Data/Sparse/SpMatrix.hs view
@@ -1,4 +1,5 @@-{-# language FlexibleInstances, MultiParamTypeClasses, TypeFamilies #-}+{-# language FlexibleInstances, FlexibleContexts, MultiParamTypeClasses, TypeFamilies #-}+{-# language DeriveFunctor, DeriveFoldable #-} ----------------------------------------------------------------------------- -- | -- Copyright : (C) 2016 Marco Zocca@@ -10,18 +11,28 @@ ----------------------------------------------------------------------------- module Data.Sparse.SpMatrix where +import Control.Exception.Common+import Data.Sparse.SpVector import Data.Sparse.Utils import Data.Sparse.Types import Numeric.Eps import Numeric.LinearAlgebra.Class +import Data.Sparse.Internal.IntM (IntM (..))+import qualified Data.Sparse.Internal.IntM as I import Data.Sparse.Internal.IntMap2 -import qualified Data.IntMap as IM +import GHC.Exts++import qualified Data.IntMap.Strict as IM++import Data.Complex+import Data.Foldable (foldl') import Data.Maybe +import Data.VectorSpace -- *@@ -40,37 +51,52 @@ -- * Sparse Matrix -data SpMatrix a = SM {smDim :: (Rows, Cols),- smData :: IM.IntMap (IM.IntMap a)} deriving Eq+data SpMatrix a = SM {smDim :: {-# UNPACK #-} !(Rows, Cols),+ smData :: !(IntM (IntM a))}+ deriving (Eq, Functor, Foldable) -sizeStr :: SpMatrix a -> String++sizeStr :: (FDSize f ~ (a1, a2), Sparse f a, Show a2, Show a1) => f a -> String sizeStr sm =- unwords ["(",show (nrows sm),"rows,",show (ncols sm),"columns ) ,",show nz,"NZ ( sparsity",show sy,")"] where- (SMInfo nz sy) = infoSM sm+ unwords ["(",show nr,"rows,",show nc,"columns ) ,",show nz,"NZ ( sparsity",show sy,")"] where+ (nr, nc) = dim sm+ nz = nnz sm+ sy = spy sm :: Double instance Show a => Show (SpMatrix a) where- show sm@(SM _ x) = "SM " ++ sizeStr sm ++ " "++ show (IM.toList x)+ show sm@(SM _ x) = unwords ["SM",sizeStr sm,show (toList $ toList <$> x)]+ -- show sm@(SM _ x) = show x -instance Functor SpMatrix where- fmap f (SM d md) = SM d ((fmap . fmap) f md)+-- instance Functor SpMatrix where+-- fmap f (SM d md) = SM d ((fmap . fmap) f md) ++ + instance Set SpMatrix where liftU2 f2 (SM n1 x1) (SM n2 x2) = SM (maxTup n1 n2) ((liftU2.liftU2) f2 x1 x2) liftI2 f2 (SM n1 x1) (SM n2 x2) = SM (minTup n1 n2) ((liftI2.liftI2) f2 x1 x2) --- | 'SpMatrix'es form a ring, in that they can be added and possess a zero element -instance Additive SpMatrix where- zero = SM (0,0) IM.empty+-- | 'SpMatrix'es form an additive group, in that they can have an invertible associtative operation (matrix sum)+instance Num a => AdditiveGroup (SpMatrix a) where+ zeroV = SM (0,0) I.empty (^+^) = liftU2 (+)+ negateV = fmap negate+ (^-^) = liftU2 (-) ++++ -- | 'SpMatrix'es are maps between finite-dimensional spaces instance FiniteDim SpMatrix where type FDSize SpMatrix = (Rows, Cols) dim = smDim instance HasData SpMatrix a where- type HDData SpMatrix a = IM.IntMap (IM.IntMap a)+ type HDData SpMatrix a = IntM (IntM a)+ nnz = nzSM dat = smData instance Sparse SpMatrix a where@@ -87,11 +113,16 @@ +++++ -- ** Creation -- | `zeroSM m n` : Empty SpMatrix of size (m, n) zeroSM :: Rows -> Cols -> SpMatrix a-zeroSM m n = SM (m,n) IM.empty+zeroSM m n = SM (m,n) I.empty -- *** Diagonal matrix@@ -173,7 +204,7 @@ -- | Add to existing SpMatrix using data from list (row, col, value) fromListSM' :: Foldable t => t (IxRow, IxCol, a) -> SpMatrix a -> SpMatrix a-fromListSM' iix sm = foldl ins sm iix where+fromListSM' iix sm = foldl' ins sm iix where ins t (i,j,x) = insertSpMatrix i j x t -- | Create new SpMatrix using data from list (row, col, value)@@ -181,6 +212,17 @@ fromListSM (m,n) iix = fromListSM' iix (zeroSM m n) +mkSpMR :: Foldable t =>+ (Int, Int) -> t (IxRow, IxCol, Double) -> SpMatrix Double+mkSpMR d ixv = fromListSM d ixv :: SpMatrix Double++mkSpMC :: Foldable t =>+ (Int, Int) -> t (IxRow, IxCol, Complex Double) -> SpMatrix (Complex Double)+mkSpMC d ixv = fromListSM d ixv :: SpMatrix (Complex Double)++++ -- | Create new SpMatrix assuming contiguous, 0-based indexing of elements fromListDenseSM :: Int -> [a] -> SpMatrix a fromListDenseSM m ll = fromListSM (m, n) $ denseIxArray2 m ll where@@ -189,6 +231,12 @@ -- ** toList +-- | Populate list with SpMatrix contents+toListSM :: SpMatrix t -> [(IxRow, IxCol, t)]+toListSM = ifoldlSM buildf [] where+ buildf i j x y = (i, j, x) : y ++ -- | Populate list with SpMatrix contents and populate missing entries with 0 toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)] toDenseListSM m =@@ -203,15 +251,13 @@ -- -- ** Lookup lookupSM :: SpMatrix a -> IxRow -> IxCol -> Maybe a-lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j+lookupSM (SM _ im) i j = I.lookup i im >>= I.lookup j -- | Looks up an element in the matrix with a default (if the element is not found, zero is returned) @@ -254,6 +300,7 @@ + -- | Extract a submatrix given the specified index bounds, rebalancing keys with the two supplied functions extractSubmatrixSM :: (IM.Key -> IM.Key) -> -- row index function@@ -265,7 +312,7 @@ | q = SM (m', n') imm' | otherwise = error $ "extractSubmatrixSM : invalid index " ++ show (i1, i2) ++ ", " ++ show (j1, j2) where imm' = mapKeysIM2 fi gi $ -- rebalance keys- IM.filter (not . IM.null) $ -- remove all-null rows+ I.filterI (not . null) $ -- remove all-null rows ifilterIM2 ff im -- keep `submatrix` ff i j _ = i1 <= i && i <= i2 &&@@ -346,10 +393,10 @@ -- |Is the matrix diagonal? isDiagonalSM :: SpMatrix a -> Bool-isDiagonalSM m = IM.size d == nrows m where- d = IM.filterWithKey ff (immSM m)- ff irow row = IM.size row == 1 &&- IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1+isDiagonalSM m = I.size d == nrows m where+ d = I.filterWithKey ff (immSM m)+ ff irow row = I.size row == 1 &&+ I.size (I.filterWithKey (\j _ -> j == irow) row) == 1 -- | Is the matrix lower/upper triangular? isLowerTriSM, isUpperTriSM :: Eq a => SpMatrix a -> Bool@@ -359,11 +406,10 @@ isUpperTriSM m = m == lm where lm = ifilterSM (\i j _ -> i <= j) m --- |Is the matrix orthogonal? i.e. Q^t ## Q == I-isOrthogonalSM :: (Eq a, Epsilon a) => SpMatrix a -> Bool+-- -- -- |Is the matrix orthogonal? i.e. Q^t ## Q == I+-- -- isOrthogonalSM :: (Eq a, Epsilon a) => SpMatrix a -> Bool isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where- rsm = roundZeroOneSM $ transposeSM sm ## sm-+ rsm = roundZeroOneSM $ transpose sm ## sm @@ -378,7 +424,7 @@ -- ** Matrix data and metadata -- | Data in internal representation (do not export)-immSM :: SpMatrix t -> IM.IntMap (IM.IntMap t)+-- immSM :: SpMatrix t -> IM.IntMap (IM.IntMap t) immSM (SM _ imm) = imm -- | (Number of rows, Number of columns)@@ -404,7 +450,7 @@ infoSM s = SMInfo (nzSM s) (spySM s) nzSM :: SpMatrix a -> Int-nzSM s = sum $ fmap IM.size (immSM s)+nzSM s = sum $ fmap I.size (immSM s) spySM :: Fractional b => SpMatrix a -> b spySM s = fromIntegral (nzSM s) / fromIntegral (nelSM s)@@ -417,7 +463,7 @@ nzRow s i | inBounds0 (nrows s) i = nzRowU s i | otherwise = error "nzRow : index out of bounds" where nzRowU :: SpMatrix a -> IM.Key -> Int- nzRowU s i = maybe 0 IM.size (IM.lookup i $ immSM s)+ nzRowU s i = maybe 0 I.size (I.lookup i $ immSM s) @@ -432,12 +478,12 @@ bwBoundsSM :: SpMatrix a -> (Int, Int) bwBoundsSM s = -- b- (snd $ IM.findMin b,- snd $ IM.findMax b)+ (snd $ I.findMin b,+ snd $ I.findMax b) where ss = immSM s- fmi = fst . IM.findMin- fma = fst . IM.findMax+ fmi = fst . I.findMin+ fma = fst . I.findMax b = fmap (\x -> fma x - fmi x + 1:: Int) ss @@ -474,11 +520,7 @@ --- encode :: (Int, Int) -> (Rows, Cols) -> Int--- encode (nr,_) (i,j) = i + (j * nr) --- decode :: (Int, Int) -> Int -> (Rows, Cols)--- decode (nr, _) ci = (r, c) where (c,r ) = quotRem ci nr @@ -499,12 +541,12 @@ -- | Vertical stacking vertStackSM, (-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a-vertStackSM mm1 mm2 = SM (m, n) $ IM.union u1 u2 where+vertStackSM mm1 mm2 = SM (m, n) $ I.union u1 u2 where nro1 = nrows mm1 m = nro1 + nrows mm2 n = max (ncols mm1) (ncols mm2) u1 = immSM mm1- u2 = IM.mapKeys (+ nro1) (immSM mm2)+ u2 = I.mapKeys (+ nro1) (immSM mm2) (-=-) = vertStackSM @@ -517,12 +559,14 @@ ------+-- | Assembles a square matrix from a list of square matrices, arranging these along the main diagonal+fromBlocksDiag :: [SpMatrix a] -> SpMatrix a+fromBlocksDiag mml = fromListSM (n, n) lstot where+ dims = map nrows mml+ n = sum dims+ shifts = init $ scanl (+) 0 dims+ lstot = concat $ zipWith shiftDims shifts $ map toListSM mml --lsts+ shiftDims s = map (\(i,j,x) -> (i + s, j + s, x)) @@ -533,16 +577,20 @@ -- ** Misc. SpMatrix operations +-- | Indexed filter over SpMatrix+{-# INLINE ifilterSM #-} ifilterSM :: (IM.Key -> IM.Key -> a -> Bool) -> SpMatrix a -> SpMatrix a ifilterSM f (SM d im) = SM d $ ifilterIM2 f im -- | Left fold over SpMatrix+{-# INLINE foldlSM #-} foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b foldlSM f n (SM _ m)= foldlIM2 f n m -- | Indexed left fold over SpMatrix+{-# INLINE ifoldlSM #-} ifoldlSM :: (IM.Key -> IM.Key -> a -> b -> b) -> b -> SpMatrix a -> b ifoldlSM f n (SM _ m) = ifoldlIM2' f n m @@ -578,7 +626,7 @@ -- ** Sparsify : remove almost-0 elements (|x| < eps) sparsifyIM2 ::- Epsilon a => IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+ Epsilon a => I.IntM (I.IntM a) -> I.IntM (I.IntM a) sparsifyIM2 = ifilterIM2 (\_ _ x -> isNz x) -- | Sparsify an SpMatrix@@ -598,16 +646,38 @@ +++-- | Modify (row, column) keys, leaving data intact. Be careful when using this!+-- modifyKeysSM' :: (IxRow -> IxRow) -> (IxCol -> IxCol) -> SpMatrix a -> SpMatrix a+modifyKeysSM' :: (IxRow -> a) -> (IxCol -> b) -> SpMatrix c -> [(a, b, c)]+modifyKeysSM' fi fj mm = zip3 (fi <$> ii) (fj <$> jj) xx where+ (ii, jj, xx) = unzip3 $ toListSM mm+++modifyKeysSM :: (IxRow -> IxRow) -> (IxCol -> IxCol) -> SpMatrix a -> SpMatrix a+modifyKeysSM fi fj mm = fromListSM (dim mm) $ zip3 (fi <$> ii) (fj <$> jj) xx where+ (ii, jj, xx) = unzip3 $ toListSM mm++++++++++ -- * Primitive algebra operations -- ** Matrix row swap -- | Swap two rows of a SpMatrix (bounds not checked) swapRows :: IxRow -> IxRow -> SpMatrix a -> SpMatrix a-swapRows i1 i2 (SM d im) = SM d $ IM.insert i1 ro2 im' where- ro1 = im IM.! i1- ro2 = im IM.! i2- im' = IM.insert i2 ro1 im+swapRows i1 i2 (SM d im) = SM d $ I.insert i1 ro2 im' where+ ro1 = im I.! i1+ ro2 = im I.! i2+ im' = I.insert i2 ro1 im -- | Swap two rows of a SpMatrix (bounds checked) swapRowsSafe :: IxRow -> IxRow -> SpMatrix a -> SpMatrix a@@ -627,6 +697,9 @@ transposeSM :: SpMatrix a -> SpMatrix a transposeSM (SM (m, n) im) = SM (n, m) (transposeIM2 im) +-- | Hermitian conjugate+hermitianConj :: Num a => SpMatrix (Complex a) -> SpMatrix (Complex a)+hermitianConj m = conjugate <$> transposeSM m @@ -635,7 +708,7 @@ -- ** Multiply matrix by a scalar matScale :: Num a => a -> SpMatrix a -> SpMatrix a-matScale a = fmap (*a)+matScale a = fmap (* a) @@ -643,19 +716,23 @@ --- ** Frobenius norm-normFrobenius :: Floating a => SpMatrix a -> a-normFrobenius m = sqrt $ foldlSM (+) 0 m' where- m' | nrows m > ncols m = transposeSM m ## m- | otherwise = m ## transposeSM m +-- ** Trace +trace :: Num b => SpMatrix b -> b+trace m = foldlSM (+) 0 $ extractDiag m +-- ** Frobenius norm +normFrobeniusSM :: (MatrixRing (SpMatrix a), Floating a) => SpMatrix a -> a+normFrobeniusSM m = sqrt $ trace (m ##^ m) +normFrobeniusSMC ::+ (MatrixRing (SpMatrix (Complex a)), RealFloat a) => SpMatrix (Complex a) -> a+normFrobeniusSMC m = sqrt $ realPart $ trace (m ##^ m) @@ -665,39 +742,61 @@ +-- ** Matrix-matrix product +instance MatrixRing (SpMatrix Double) where+ type MatrixNorm (SpMatrix Double) = Double+ (##) = matMat_ AB+ (##^) = matMat_ ABt+ transpose = transposeSM+ normFrobenius = normFrobeniusSM --- ** Matrix-matrix product+instance MatrixRing (SpMatrix (Complex Double)) where+ type MatrixNorm (SpMatrix (Complex Double)) = Double+ (##) = matMat_ AB+ (##^) = matMat_ ABt+ transpose = hermitianConj+ normFrobenius = normFrobeniusSMC -matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a-matMat m1 m2- | c1 == r2 = matMatU m1 m2- | otherwise = error $ "matMat : incompatible matrix sizes" ++ show (d1, d2) where- d1@(r1, c1) = dim m1- d2@(r2, c2) = dim m2- matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a- matMatU m1 m2 =- SM (nrows m1, ncols m2) im where- im = fmap (\vm1 -> (`dot` vm1) <$> transposeIM2 (immSM m2)) (immSM m1)- -(##) = matMat+-- | Internal implementation+data MatProd_ = AB | ABt deriving (Eq, Show) --- matMat m1 m2 =--- withDim2 m1 m2--- (\(r1,c1) (r2,c2) _ _ -> c1 == r2)--- matMatU--- "matMat : incompatible matrix sizes"--- (\m1 m2 -> unwords [show (dim m1), show (dim m2)])+{-# INLINE matMat_ #-}+matMat_ pt mm1 mm2 =+ case pt of AB -> matMatCheck (matMatUnsafeWith transposeIM2) mm1 mm2+ ABt -> matMatCheck (matMatUnsafeWith id) mm1 (trDim mm2)+ where+ trDim (SM (a, b) x) = SM (b, a) x+ matMatCheck mmf m1 m2+ | c1 == r2 = mmf m1 m2+ | otherwise = error $ "matMat : incompatible matrix sizes" ++ show (d1, d2)+ where+ d1@(_, c1) = dim m1+ d2@(r2, _) = dim m2 +-- | Matrix product without dimension checks+{-# INLINE matMatUnsafeWith #-}+-- matMatUnsafeWith :: Num a =>+-- (IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)) ->+-- SpMatrix a ->+-- SpMatrix a ->+-- SpMatrix a+matMatUnsafeWith ff2 m1 m2 = SM (nrows m1, ncols m2) (overRows2 <$> immSM m1) where+ overRows2 vm1 = (`dott` vm1) <$> ff2 (immSM m2)+ dott x y = sum $ liftI2 (*) x y -- NB !! no complex conjugation +++ -- ** Matrix-matrix product, sparsified -- | Removes all elements `x` for which `| x | <= eps`)-matMatSparsified, (#~#) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a-matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2+matMatSparsified, (#~#) :: (MatrixRing (SpMatrix a), Epsilon a) =>+ SpMatrix a -> SpMatrix a -> SpMatrix a+matMatSparsified m1 m2 = sparsifySM $ m1 ## m2 (#~#) = matMatSparsified @@ -707,13 +806,15 @@ -- *** Sparsified matrix products of two matrices -- | A^T B-(#^#) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a-a #^# b = transposeSM a #~# b+(#~#^) :: (MatrixRing (SpMatrix a), Epsilon a) =>+ SpMatrix a -> SpMatrix a -> SpMatrix a+a #~^# b = transpose a #~# b -- | A B^T-(##^) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a-a ##^ b = a #~# transposeSM b+(#~^#) :: (MatrixRing (SpMatrix a), Epsilon a) =>+ SpMatrix a -> SpMatrix a -> SpMatrix a+a #~#^ b = a #~# transpose b @@ -730,9 +831,23 @@ -- ** Partial inner product -- | Contract row `i` of A with column `j` of B up to an index `n`, i.e. summing over repeated indices: -- Aij Bjk , for j in [0 .. n] -contractSub :: Num a => SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> Int -> a+contractSub :: Elt a => SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> Int -> a contractSub a b i j n | ncols a == nrows b && isValidIxSM a (i,j) &&- n <= ncols a = sum $ map (\i' -> a@@!(i,i')*b@@!(i',j)) [0 .. n]+ n <= ncols a = sum $ map (\i' -> (a@@!(i,i'))*b@@!(i',j)) [0 .. n] | otherwise = error "contractSub : n must be <= i"++++++++-- -- * Misc.utilities++-- encode :: (Int, Int) -> (Rows, Cols) -> Int+-- encode (nr,_) (i,j) = i + (j * nr)++-- decode :: (Int, Int) -> Int -> (Rows, Cols)+-- decode (nr, _) ci = (r, c) where (c,r ) = quotRem ci nr
src/Data/Sparse/SpVector.hs view
@@ -1,4 +1,7 @@+{-# LANGUAGE FlexibleContexts #-} {-# language TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-}+{-# language CPP #-}+{-# language GeneralizedNewtypeDeriving, DeriveFunctor #-} ----------------------------------------------------------------------------- -- | -- Copyright : (C) 2016 Marco Zocca@@ -10,32 +13,47 @@ ----------------------------------------------------------------------------- module Data.Sparse.SpVector where +import Control.Exception+import Control.Monad.Catch (MonadThrow (..))+import Control.Exception.Common++import GHC.Exts+ import Data.Sparse.Utils import Data.Sparse.Types-import Data.Sparse.Internal.IntMap2+import Data.Sparse.Internal.IntM import Numeric.Eps import Numeric.LinearAlgebra.Class -+import Data.Complex import Data.Maybe -import qualified Data.IntMap as IM+import qualified Data.IntMap.Strict as IM import qualified Data.Foldable as F import qualified Data.Vector as V +import Data.VectorSpace hiding (magnitude)+++++ -- * Sparse Vector -data SpVector a = SV { svDim :: Int ,- svData :: IM.IntMap a} deriving Eq+data SpVector a = SV { svDim :: {-# UNPACK #-} !Int ,+ svData :: !(IntM a)} deriving Eq +instance Show a => Show (SpVector a) where+ show (SV d x) = "SV (" ++ show d ++ ") "++ show (toList x)+ -- | SpVector sparsity spySV :: Fractional b => SpVector a -> b-spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (dim s)+spySV s = fromIntegral (size (dat s)) / fromIntegral (dim s) -- | Number of nonzeros nzSV :: SpVector a -> Int-nzSV sv = IM.size (dat sv)+nzSV sv = size (dat sv) sizeStrSV :: SpVector a -> String@@ -53,88 +71,127 @@ instance Foldable SpVector where foldr f d v = F.foldr f d (svData v) -instance Additive SpVector where- zero = SV 0 IM.empty- (^+^) = liftU2 (+) ++ -- | 'SpVector's form a vector space because they can be multiplied by a scalar-instance VectorSpace SpVector where- n .* v = scale n v + -- | 'SpVector's are finite-dimensional vectors instance FiniteDim SpVector where type FDSize SpVector = Int dim = svDim instance HasData SpVector a where- type HDData SpVector a = IM.IntMap a+ type HDData SpVector a = IntM a dat = svData+ nnz (SV _ x) = length x instance Sparse SpVector a where spy = spySV -- | 'SpVector's are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0)-instance Num a => SpContainer SpVector a where+instance Elt a => SpContainer SpVector a where type ScIx SpVector = Int scInsert = insertSpVector scLookup v i = lookupSV i v+ scToList = toListSV v @@ i = lookupDenseSV i v +-- instance (Elt e, RealFloat e) => SparseVector SpVector e where+-- type SpvIx SpVector = Int+-- svFromList = fromListSV+-- svFromListDense = fromListDenseSV+-- svConcat = foldr concatSV zero --- | 'SpVector's form a Hilbert space, in that we can define an inner product over them-instance Hilbert SpVector where- a `dot` b | dim a == dim b = dot (dat a) (dat b)- | otherwise =- error $ "dot : sizes must coincide, instead we got " ++- show (dim a, dim b)+-- instance SparseVector SpVector (Complex Double) where --- | Since 'SpVector's form a Hilbert space, we can define a norm for them -instance Normed SpVector where- norm p (SV _ v) = norm p v +#define SpVectorInstance(t) \+ instance AdditiveGroup (SpVector (t)) where { zeroV = SV 0 empty; (^+^) = liftU2 (+); negateV = fmap negate };\+ instance AdditiveGroup (SpVector (Complex t)) where { zeroV = SV 0 empty; (^+^) = liftU2 (+); negateV = fmap negate };\+ instance VectorSpace (SpVector t) where { type (Scalar (SpVector t)) = t; n *^ v = scale n v};\+ instance VectorSpace (SpVector (Complex t)) where { type (Scalar (SpVector (Complex t))) = Complex t; n *^ v = scale n v};\+ instance InnerSpace (SpVector (t)) where { (<.>) = dotS };\+ instance InnerSpace (SpVector (Complex (t))) where { (<.>) = dotS };\+ instance Normed (SpVector (t)) where {type RealScalar (SpVector (t)) = t; type Magnitude (SpVector (t)) = t; norm1 (SV _ v) = norm1 v; norm2Sq (SV _ v) = norm2Sq v ; normP p (SV _ v) = normP p v; normalize p (SV n v) = SV n (normalize p v); normalize2 (SV n v) = SV n (normalize2 v)};\+ instance Normed (SpVector (Complex t)) where {type RealScalar (SpVector (Complex t)) = t; type Magnitude (SpVector (Complex t)) = t; norm1 (SV _ v) = norm1 v; norm2Sq (SV _ v) = norm2Sq v ; normP p (SV _ v) = normP p v; normalize p (SV n v) = SV n (normalize p v); normalize2 (SV n v) = SV n (normalize2 v)} +SpVectorInstance(Double)+-- SpVectorInstance(Float) +dotS :: InnerSpace (IntM t) => SpVector t -> SpVector t -> Scalar (IntM t)+(SV m a) `dotS` (SV n b)+ | n == m = a <.> b+ | otherwise = error $ unwords ["<.> : Incompatible dimensions:", show m, show n]++-- dotSSafe :: (MonadThrow m, InnerSpace (IM.IntMap t)) =>+-- SpVector t -> SpVector t -> m (Scalar (IM.IntMap t))+dotSSafe :: (InnerSpace (IntM t), MonadThrow m) =>+ SpVector t -> SpVector t -> m (Scalar (IntM t))+dotSSafe (SV m a) (SV n b)+ | n == m = return $ a <.> b+ | otherwise = throwM (DotSizeMismatch m n)++++++++++ -- ** Creation -- | Empty sparse vector (length n, no entries) zeroSV :: Int -> SpVector a-zeroSV n = SV n IM.empty+zeroSV n = SV n empty -- | Singleton sparse vector (length 1) singletonSV :: a -> SpVector a-singletonSV x = SV 1 (IM.singleton 0 x)+singletonSV x = SV 1 (singleton 0 x) -- | Canonical basis vector in R^n ei :: Num a => Int -> IM.Key -> SpVector a-ei n i = SV n (IM.insert (i - 1) 1 IM.empty)+ei n i = SV n (insert (i - 1) 1 empty) --- | create a sparse vector from an association list while discarding all zero entries+-- | Sparse vector from an association list while discarding all zero entries mkSpVector :: Epsilon a => Int -> IM.IntMap a -> SpVector a-mkSpVector d im = SV d $ IM.filterWithKey (\k v -> isNz v && inBounds0 d k) im+mkSpVector d im = SV d $ IntM $ IM.filterWithKey (\k v -> isNz v && inBounds0 d k) im --- | ", from logically dense array (consecutive indices)-mkSpVectorD :: Epsilon a => Int -> [a] -> SpVector a-mkSpVectorD d ll = mkSpVector d (IM.fromList $ denseIxArray (take d ll)) -- ", don't filter zero elements mkSpVector1 :: Int -> IM.IntMap a -> SpVector a-mkSpVector1 d ll = SV d $ IM.filterWithKey (\ k _ -> inBounds0 d k) ll+mkSpVector1 d ll = SV d $ IntM $ IM.filterWithKey (\ k _ -> inBounds0 d k) ll ++-- | Dense real SpVector (monomorphic Double)+mkSpVR :: Int -> [Double] -> SpVector Double+mkSpVR d ll = SV d $ mkIm ll++-- | Dense complex SpVector (monomorphic Double)+mkSpVC :: Int -> [Complex Double] -> SpVector (Complex Double)+mkSpVC d ll = SV d $ mkImC ll++++ -- | Create new sparse vector, assumin 0-based, contiguous indexing fromListDenseSV :: Int -> [a] -> SpVector a-fromListDenseSV d ll = SV d (IM.fromList $ denseIxArray (take d ll))+fromListDenseSV d ll = SV d (fromList $ denseIxArray (take d ll)) -- | Map a function over a range of indices and filter the result (indices and values) to fit in a `n`-long SpVector@@ -155,7 +212,7 @@ -- | one-hot encoding : `oneHotSV n k` produces a SpVector of length n having 1 at the k-th position oneHotSVU :: Num a => Int -> IxRow -> SpVector a-oneHotSVU n k = SV n (IM.singleton k 1)+oneHotSVU n k = SV n (singleton k 1) oneHotSV :: Num a => Int -> IxRow -> SpVector a oneHotSV n k |inBounds0 n k = oneHotSVU n k@@ -164,11 +221,11 @@ -- | DENSE vector of `1`s onesSV :: Num a => Int -> SpVector a-onesSV d = SV d $ IM.fromList $ denseIxArray $ replicate d 1+onesSV d = SV d $ fromList $ denseIxArray $ replicate d 1 -- | DENSE vector of `0`s zerosSV :: Num a => Int -> SpVector a-zerosSV d = SV d $ IM.fromList $ denseIxArray $ replicate d 0+zerosSV d = SV d $ fromList $ denseIxArray $ replicate d 0 @@ -196,36 +253,40 @@ -- ** Element insertion --- |insert element `x` at index `i` in a preexisting SpVector-insertSpVector :: Int -> a -> SpVector a -> SpVector a-insertSpVector i x (SV d xim)- | inBounds0 d i = SV d (IM.insert i x xim)- | otherwise = error "insertSpVector : index out of bounds"+-- |insert element `x` at index `i` in a preexisting SpVector; discards out-of-bounds entries+insertSpVector :: IM.Key -> a -> SpVector a -> SpVector a+insertSpVector i x (SV d xim) | inBounds0 d i = SV d (insert i x xim) +insertSpVectorSafe :: MonadThrow m => Int -> a -> SpVector a -> m (SpVector a)+insertSpVectorSafe i x (SV d xim)+ | inBounds0 d i = return $ SV d (insert i x xim)+ | otherwise = throwM (OOBIxError "insertSpVector" i) +++ -- ** fromList-fromListSV :: Int -> [(Int, a)] -> SpVector a-fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix ))+fromListSV :: Foldable t => Int -> t (Int, a) -> SpVector a+fromListSV d iix = SV d $ foldr insf empty iix where+ insf (i, x) xacc | inBounds0 d i = insert i x xacc+ | otherwise = xacc --- fromListSV' d iix = SV d (F.foldl' insf IM.empty iix') where--- insf im (i, x) = IM.insert i x im--- iix' = filter (inBounds0 d . fst) iix -- filtering forces list as only instance -- ** toList toListSV :: SpVector a -> [(IM.Key, a)]-toListSV sv = IM.toList (dat sv)+toListSV sv = toList (dat sv) -- |To dense list (default = 0) toDenseListSV :: Num b => SpVector b -> [b]-toDenseListSV (SV d im) = fmap (\i -> IM.findWithDefault 0 i im) [0 .. d-1]+toDenseListSV (SV d (IntM im)) = fmap (\i -> IM.findWithDefault 0 i im) [0 .. d-1] -- | Indexed fold over SpVector ifoldSV :: (IM.Key -> a -> b -> b) -> b -> SpVector a -> b-ifoldSV f e (SV d im) = IM.foldWithKey f e im+ifoldSV f e (SV _ (IntM im)) = IM.foldrWithKey f e im @@ -234,19 +295,18 @@ -instance Show a => Show (SpVector a) where- show (SV d x) = "SV (" ++ show d ++ ") "++ show (IM.toList x) + -- ** Lookup -- | Lookup an index in a SpVector lookupSV :: IM.Key -> SpVector a -> Maybe a-lookupSV i (SV _ im) = IM.lookup i im+lookupSV i (SV _ (IntM im)) = IM.lookup i im -- | Lookup an index, return a default value if lookup fails lookupDefaultSV :: a -> IM.Key -> SpVector a -> a-lookupDefaultSV def i (SV _ im) = IM.findWithDefault def i im+lookupDefaultSV def i (SV _ (IntM im)) = IM.findWithDefault def i im -- |Lookup an index in a SpVector, returns 0 if lookup fails lookupDenseSV :: Num a => IM.Key -> SpVector a -> a@@ -260,52 +320,69 @@ -- ** Sub-vectors -- | Tail elements tailSV :: SpVector a -> SpVector a-tailSV (SV n sv) = SV (n-1) ta where+tailSV (SV n (IntM sv)) = SV (n-1) $ IntM ta where ta = IM.mapKeys (\i -> i - 1) $ IM.delete 0 sv -- | Head element headSV :: Num a => SpVector a -> a-headSV sv = fromMaybe 0 (IM.lookup 0 (dat sv))+headSV (SV _ (IntM im)) = fromMaybe 0 (IM.lookup 0 im) -- | Keep the first n components of the SpVector (like `take` for lists) takeSV, dropSV :: Int -> SpVector a -> SpVector a-takeSV n (SV _ sv) = SV n $ IM.filterWithKey (\i _ -> i < n) sv+takeSV n (SV _ sv) = SV n $ filterWithKey (\i _ -> i < n) sv -- | Discard the first n components of the SpVector and rebalance the keys (like `drop` for lists)-dropSV n (SV n0 sv) = SV (n0 - n) $ IM.mapKeys (subtract n) $ IM.filterWithKey (\i _ -> i >= n) sv+dropSV n (SV n0 (IntM sv)) = SV (n0 - n) $ IntM $ IM.mapKeys (subtract n) $ IM.filterWithKey (\i _ -> i >= n) sv +-- | Keep a range of entries +rangeSV :: (IM.Key, IM.Key) -> SpVector a -> SpVector a+rangeSV (rmin, rmax) (SV n (IntM sv))+ | len > 0 && len <= n = SV len $ IntM sv'+ | otherwise = error $ unwords ["rangeSV : invalid bounds", show (rmin, rmax) ] where+ len = rmax - rmin+ sv' = IM.mapKeys (subtract rmin) $ IM.filterWithKey (\i _ -> i >= rmin && i <= rmax) sv +++ -- | Concatenate two sparse vectors concatSV :: SpVector a -> SpVector a -> SpVector a-concatSV (SV n1 s1) (SV n2 s2) = SV (n1+n2) (IM.union s1 s2') where+concatSV (SV n1 (IntM s1)) (SV n2 (IntM s2)) = SV (n1+n2) $ IntM (IM.union s1 s2') where s2' = IM.mapKeys (+ n1) s2 -- | Filter filterSV :: (a -> Bool) -> SpVector a -> SpVector a-filterSV q sv = SV (dim sv) (IM.filter q (dat sv)) +filterSV q sv = SV (dim sv) $ IntM (IM.filter q (unIM $ dat sv)) -- | Indexed filter ifilterSV :: (Int -> a -> Bool) -> SpVector a -> SpVector a-ifilterSV q sv = SV (dim sv) (IM.filterWithKey q (dat sv))+ifilterSV q sv = SV (dim sv) (filterWithKey q (dat sv)) +-- * Sparsify : remove almost-0 elements (|x| < eps)+-- | Sparsify an SpVector+sparsifySV :: Epsilon a => SpVector a -> SpVector a+sparsifySV = filterSV isNz ++ -- * Orthogonal vector -- | Generate an arbitrary (not random) vector `u` such that `v dot u = 0`-orthogonalSV :: Fractional a => SpVector a -> SpVector a+orthogonalSV :: (Scalar (SpVector t) ~ t, InnerSpace (SpVector t), Fractional t) =>+ SpVector t -> SpVector t orthogonalSV v = u where (h, t) = (headSV v, tailSV v) n = dim v
src/Data/Sparse/Types.hs view
@@ -5,3 +5,9 @@ type IxRow = Int type IxCol = Int++-- * Lexicographic ordering types++type LexIx = Int++data LexOrd = RowsFirst | ColsFirst deriving (Eq, Show)
src/Data/Sparse/Utils.hs view
@@ -33,6 +33,11 @@ n = ln `div` m +-- indexed :: [b] -> [(Int, b)]+-- indexed x = zip [0 .. length x - 1] x+++ -- folds -- | foldr over the results of a fmap@@ -84,3 +89,26 @@ tail' :: V.Vector a -> Maybe (V.Vector a) tail' = harness V.null V.tail+++++-- | a cons-based moving-window datatype of length at least 3++data W3 a = W3 a a a [a] deriving (Eq, Show)+initW3 :: a -> a -> a -> [a] -> W3 a+initW3 = W3++pushW3 :: a -> W3 a -> W3 a+pushW3 i (W3 i0 i1 i2 is) = W3 i i0 i1 (i2 : initSafe is)++fstW3 (W3 i _ _ _) = i+sndW3 (W3 _ i _ _) = i+thirdW3 (W3 _ _ i _) = i++withInitW3, withTailW3 :: W3 t -> (t -> t -> t1) -> t1+withInitW3 (W3 a b _ _) f = f a b+withTailW3 (W3 _ b c _) f = f b c++initSafe (x:xs) = x : initSafe xs+initSafe [] = []
src/Numeric/Eps.hs view
@@ -1,3 +1,4 @@+{-# language FlexibleInstances #-} ----------------------------------------------------------------------------- -- | -- Copyright : (C) 2016 Marco Zocca, 2012-2015 Edward Kmett@@ -9,8 +10,10 @@ -- Testing for values "near" zero ----------------------------------------------------------------------------- module Numeric.Eps- ( Epsilon(..), nearZero, isNz, roundZero, roundOne, roundZeroOne+ ( Epsilon(..), isNz, roundZero, roundOne, roundZeroOne,+ nearOne ) where+import Data.Complex import Foreign.C.Types (CFloat, CDouble) -- | Provides a test to see if a quantity is near zero.@@ -26,7 +29,7 @@ -- -- >>> nearZero (1e-7 :: Float) -- True-class Num a => Epsilon a where+class (Floating a, Num a) => Epsilon a where -- | Determine if a quantity is near zero. nearZero :: a -> Bool @@ -47,24 +50,37 @@ nearZero a = abs a <= 1e-12 +-- | Complex types +instance Epsilon (Complex Float) where+ nearZero a = magnitude a <= 1e-6 +instance Epsilon (Complex Double) where+ nearZero a = magnitude a <= 1e-12++instance Epsilon (Complex CFloat) where+ nearZero a = magnitude a <= 1e-6++instance Epsilon (Complex CDouble) where+ nearZero a = magnitude a <= 1e-12+++ -- * Rounding operations -- | Rounding rule-almostZero, almostOne, isNz :: Epsilon a => a -> Bool-almostZero = nearZero-almostOne x = nearZero (1 - x)-isNz x = not (almostZero x)+nearOne, isNz :: Epsilon a => a -> Bool+nearOne x = nearZero (1 - x)+isNz x = not (nearZero x) withDefault :: (t -> Bool) -> t -> t -> t withDefault q d x | q x = d | otherwise = x roundZero, roundOne, roundZeroOne :: Epsilon a => a -> a-roundZero = withDefault almostZero (fromIntegral 0)-roundOne = withDefault almostOne (fromIntegral 1)+roundZero = withDefault nearZero (fromIntegral (0 :: Int))+roundOne = withDefault nearOne (fromIntegral (1 :: Int)) with2Defaults :: (t -> Bool) -> (t -> Bool) -> t -> t -> t -> t with2Defaults q1 q2 d1 d2 x | q1 x = d1@@ -72,4 +88,4 @@ | otherwise = x -- | Round to respectively 0 or 1-roundZeroOne = with2Defaults almostZero almostOne (fromIntegral 0) (fromIntegral 1)+roundZeroOne = with2Defaults nearZero nearOne (fromIntegral (0 :: Int)) (fromIntegral (1 :: Int))
src/Numeric/LinearAlgebra/Class.hs view
@@ -1,111 +1,187 @@-{-# LANGUAGE TypeFamilies, MultiParamTypeClasses, KindSignatures #-}+{-# language TypeFamilies, MultiParamTypeClasses, KindSignatures, FlexibleContexts, FlexibleInstances, ConstraintKinds #-}+{-# language AllowAmbiguousTypes #-}+{-# language CPP #-}+-----------------------------------------------------------------------------+-- |+-- Module : Numeric.LinearAlgebra.Class+-- Copyright : (c) Marco Zocca 2017+-- License : GPL-style (see the file LICENSE)+--+-- Maintainer : zocca marco gmail+-- Stability : experimental+-- Portability : portable+--+-- Typeclasses for linear algebra and related concepts+--+----------------------------------------------------------------------------- module Numeric.LinearAlgebra.Class where --- * Additive ring -class Functor f => Additive f where- -- | Ring zero element- zero :: Num a => f a- - -- | Ring +- (^+^) :: Num a => f a -> f a -> f a+-- import Control.Applicative+import Data.Complex+-- import Data.Ratio+-- import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble) +-- import Control.Exception+-- import Control.Exception.Common+import Control.Monad.Catch+import Control.Monad.IO.Class - one :: Num a => f a+-- import Data.Typeable (Typeable) - (^*^) :: Num a => f a -> f a -> f a+import qualified Data.Vector as V (Vector) +import Data.VectorSpace hiding (magnitude) +import Data.Sparse.Types+import Numeric.Eps --- | negate the values in a functor-negated :: (Num a, Functor f) => f a -> f a-negated = fmap negate --- | subtract two Additive objects-(^-^) :: (Additive f, Num a) => f a -> f a -> f a-x ^-^ y = x ^+^ negated y +-- * Matrix and vector elements (possibly Complex)+class (Eq e , Fractional e, Floating e, Num (EltMag e), Ord (EltMag e)) => Elt e where+ type EltMag e :: *+ conj :: e -> e+ conj = id+ mag :: e -> EltMag e +instance Elt Double where {type EltMag Double = Double ; mag = id}+instance Elt Float where {type EltMag Float = Float; mag = id}+instance RealFloat e => Elt (Complex e) where+ type EltMag (Complex e) = e+ conj = conjugate+ mag = magnitude+ ++ -- * Vector space-class Additive f => VectorSpace f where- -- | multiplication by a scalar- (.*) :: Num a => a -> f a -> f a- --- |linear interpolation-lerp :: (VectorSpace f, Num a) => a -> f a -> f a -> f a+-- | Scale a vector +(.*) :: VectorSpace v => Scalar v -> v -> v+(.*) = (*^)++-- | Scale a vector by the reciprocal of a number (e.g. for normalization)+(./) :: (VectorSpace v, Fractional (Scalar v)) => v -> Scalar v -> v+v ./ n = recip n .* v++-- | Convex combination of two vectors (NB: 0 <= `a` <= 1). +lerp :: (VectorSpace e, Num (Scalar e)) => Scalar e -> e -> e -> e lerp a u v = a .* u ^+^ ((1-a) .* v) +-- linearCombination :: (VectorSpace v , Foldable t) => t (Scalar v, v) -> v+-- linearCombination = foldr (\(a, x) (b, y) -> (a .* x) ^+^ (b .* y)) ++-- linComb a v = a .* v++++ -- * Hilbert space (inner product)-class VectorSpace f => Hilbert f where- -- | inner product- dot :: Num a => f a -> f a -> a +-- | Inner product+dot :: InnerSpace v => v -> v -> Scalar v+dot = (<.>)+ ++++ -- ** Hilbert-space distance function--- |`hilbertDistSq x y = || x - y ||^2`-hilbertDistSq :: (Hilbert f, Num a) => f a -> f a -> a-hilbertDistSq x y = dot t t where++-- |`hilbertDistSq x y = || x - y ||^2` computes the squared L2 distance between two vectors+hilbertDistSq :: InnerSpace v => v -> v -> Scalar v+hilbertDistSq x y = t <.> t where t = x ^-^ y - --- * Normed vector space-class Hilbert f => Normed f where- norm :: (Floating a, Eq a) => a -> f a -> a --- ** Norms and related results --- | Squared 2-norm-normSq :: (Hilbert f, Num a) => f a -> a-normSq v = v `dot` v+-- * Normed vector spaces +class (InnerSpace v, Num (RealScalar v), Eq (RealScalar v), Epsilon (Magnitude v), Show (Magnitude v), Ord (Magnitude v)) => Normed v where+ type Magnitude v :: *+ type RealScalar v :: *+ norm1 :: v -> Magnitude v -- ^ L1 norm+ norm2Sq :: v -> Magnitude v -- ^ Euclidean (L2) norm+ normP :: RealScalar v -> v -> Magnitude v -- ^ Lp norm (p > 0)+ normalize :: RealScalar v -> v -> v -- ^ Normalize w.r.t. Lp norm+ normalize2 :: v -> v -- ^ Normalize w.r.t. L2 norm+ normalize2' :: Floating (Scalar v) => v -> v -- ^ Normalize w.r.t. norm2' instead of norm2+ normalize2' x = x ./ norm2' x+ norm2 :: Floating (Magnitude v) => v -> Magnitude v -- ^ Euclidean norm+ norm2 x = sqrt (norm2Sq x)+ norm2' :: Floating (Scalar v) => v -> Scalar v -- ^ Euclidean norm; returns a Complex (norm :+ 0) for containers of complex values+ norm2' x = sqrt $ x <.> x+ norm :: Floating (Magnitude v) => RealScalar v -> v -> Magnitude v -- ^ Lp norm (p > 0)+ norm p v+ | p == 1 = norm1 v+ | p == 2 = norm2 v+ | otherwise = normP p v --- |L1 norm-norm1 :: (Foldable t, Num a, Functor t) => t a -> a-norm1 v = sum (fmap abs v) --- |Euclidean norm-norm2 :: (Hilbert f, Floating a) => f a -> a-norm2 v = sqrt (normSq v) --- |Lp norm (p > 0)-normP :: (Foldable t, Functor t, Floating a) => a -> t a -> a-normP p v = sum u**(1/p) where- u = fmap (**p) v+-- | Infinity-norm (Real)+normInftyR :: (Foldable t, Ord a) => t a -> a+normInftyR x = maximum x --- |Infinity-norm-normInfty :: (Foldable t, Ord a) => t a -> a-normInfty = maximum+-- | Infinity-norm (Complex)+normInftyC :: (Foldable t, RealFloat a, Functor t) => t (Complex a) -> a+normInftyC x = maximum (magnitude <$> x) --- |Normalize w.r.t. p-norm (p finite)-normalize :: (Normed f, Floating a, Eq a) => a -> f a -> f a-normalize p v = (1 / norm p v) .* v +instance Normed Double where+ type Magnitude Double = Double+ type RealScalar Double = Double+ norm1 = abs+ norm2Sq = abs+ normP _ = abs+ normalize _ _ = 1+ normalize2 _ = 1+ norm2 = abs+ norm2' = abs +instance Normed (Complex Double) where+ type Magnitude (Complex Double) = Double+ type RealScalar (Complex Double) = Double+ norm1 (r :+ i) = abs r + abs i+ norm2Sq = (**2) . magnitude+ normP p (r :+ i) = (r**p + i**p)**(1/p)+ normalize p c = c ./ normP p c+ normalize2 c = c ./ magnitude c+ norm2 = magnitude+ norm2' = magnitude+ --- |Lp inner product (p > 0)+ +++++++-- | Lp inner product (p > 0) dotLp :: (Set t, Foldable t, Floating a) => a -> t a -> t a -> a dotLp p v1 v2 = sum u**(1/p) where f a b = (a*b)**p u = liftI2 f v1 v2 --- |Reciprocal+-- | Reciprocal reciprocal :: (Functor f, Fractional b) => f b -> f b reciprocal = fmap recip @@ -120,121 +196,285 @@ ++++-- * Matrix ring++-- | A matrix ring is any collection of matrices over some ring R that form a ring under matrix addition and matrix multiplication++class (AdditiveGroup m, Epsilon (MatrixNorm m)) => MatrixRing m where+ type MatrixNorm m :: *+ (##) :: m -> m -> m+ (##^) :: m -> m -> m -- ^ A B^T+ (#^#) :: m -> m -> m -- ^ A^T B+ a #^# b = transpose a ## b+ transpose :: m -> m+ normFrobenius :: m -> MatrixNorm m+++++-- a "sparse matrix ring" ?++-- class MatrixRing m a => SparseMatrixRing m a where+-- (#~#) :: Epsilon a => Matrix m a -> Matrix m a -> Matrix m a++++++++-- * Linear vector space++class (VectorSpace v, MatrixRing (MatrixType v)) => LinearVectorSpace v where+ type MatrixType v :: *+ (#>) :: MatrixType v -> v -> v+ (<#) :: v -> MatrixType v -> v++++++-- ** LinearVectorSpace + Normed++type V v = (LinearVectorSpace v, Normed v) +++++-- ** Linear systems+ +class LinearVectorSpace v => LinearSystem v where+ (<\>) :: (MonadIO m, MonadThrow m) => MatrixType v -> v -> m v++++++++++ -- * FiniteDim : finite-dimensional objects -class Additive f => FiniteDim f where+class Functor f => FiniteDim f where type FDSize f :: * dim :: f a -> FDSize f +class FiniteDim' f where+ type FDSize' f :: *+ dim' :: f -> FDSize' f --- | unary dimension-checking bracket-withDim :: (FiniteDim f, Show e) =>- f a- -> (FDSize f -> f a -> Bool)- -> (f a -> c)- -> String- -> (f a -> e)- -> c-withDim x p f e ef | p (dim x) x = f x- | otherwise = error e' where e' = e ++ show (ef x) --- | binary dimension-checking bracket-withDim2 :: (FiniteDim f, FiniteDim g, Show e) =>- f a- -> g b- -> (FDSize f -> FDSize g -> f a -> g b -> Bool)- -> (f a -> g b -> c)- -> String- -> (f a -> g b -> e)- -> c-withDim2 x y p f e ef | p (dim x) (dim y) x y = f x y- | otherwise = error e' where e' = e ++ show (ef x y)+-- -- | unary dimension-checking bracket+-- withDim :: (FiniteDim f, Show s) =>+-- f e+-- -> (FDSize f -> f e -> Bool)+-- -> (f e -> c)+-- -> String+-- -> (f e -> s)+-- -> c+-- withDim x p f e ef | p (dim x) x = f x+-- | otherwise = error e' where e' = e ++ show (ef x) +-- -- | binary dimension-checking bracket+-- withDim2 :: (FiniteDim f, FiniteDim g, Show s) =>+-- f e+-- -> g e+-- -> (FDSize f -> FDSize g -> f e -> g e -> Bool)+-- -> (f e -> g e -> c)+-- -> String+-- -> (f e -> g e -> s)+-- -> c+-- withDim2 x y p f e ef | p (dim x) (dim y) x y = f x y+-- | otherwise = error e' where e' = e ++ show (ef x y) + -- * HasData : accessing inner data (do not export) -class Additive f => HasData f a where- type HDData f a :: * +class HasData f a where+ type HDData f a :: *+ nnz :: f a -> Int dat :: f a -> HDData f a +class HasData' f where+ type HDD f :: *+ nnz' :: f -> Int+ dat' :: f -> HDD f + -- * Sparse : sparse datastructures class (FiniteDim f, HasData f a) => Sparse f a where spy :: Fractional b => f a -> b +class (FiniteDim' f, HasData' f) => Sparse' f where+ spy' :: Fractional b => f -> b -- * Set : types that behave as sets class Functor f => Set f where- -- |union binary lift : apply function on _union_ of two Sets+ -- |union binary lift : apply function on _union_ of two "sets" liftU2 :: (a -> a -> a) -> f a -> f a -> f a - -- |intersection binary lift : apply function on _intersection_ of two Sets- liftI2 :: (a -> b -> c) -> f a -> f b -> f c -+ -- |intersection binary lift : apply function on _intersection_ of two "sets"+ liftI2 :: (a -> a -> b) -> f a -> f a -> f b +-- * SpContainer : sparse container datastructures. Insertion, lookup, toList, lookup with 0 default class Sparse c a => SpContainer c a where type ScIx c :: * scInsert :: ScIx c -> a -> c a -> c a scLookup :: c a -> ScIx c -> Maybe a+ scToList :: c a -> [(ScIx c, a)] -- -- | Lookup with default, infix form ("safe" : should throw an exception if lookup is outside matrix bounds) (@@) :: c a -> ScIx c -> a +class SpContainer' c where+ type ScIx' c :: *+ scInsert' :: ScIx' c -> a -> c -> c+ scLookup' :: c -> ScIx' c -> Maybe a+ scToList' :: c -> [a]+ -- (@@') --- * IxContainer : indexed container types --- class IxContainer (c :: * -> *) a where--- type Ix c :: *--- type IxSz c :: *--- ixcLookup :: c a -> Ix c -> Maybe a--- ixcIfilter :: (Ix c -> a -> Bool) -> c a -> c a--- ixcInsert :: Ix c -> a -> c a -> c a--- ixcFromList :: Foldable t => IxSz c -> t (Ix c, a) -> c a--- ixcToList :: c a -> [(Ix c, a)] --- newtype IM_ a = IM (IM.IntMap a) --- instance IxContainer IM_ a where--- type Ix IM_ = Int--- -- -- ixcLookupDefault = lookupDefaultSV--- -- -- ixcFilter = filterSV +-- * SparseVector --- newtype IM2 a = IM2 { unIM2 :: IM.IntMap (IM.IntMap a)}+class SpContainer v e => SparseVector v e where+ type SpvIx v :: *+ svFromList :: Int -> [(SpvIx v, e)] -> v e+ svFromListDense :: Int -> [e] -> v e+ svConcat :: Foldable t => t (v e) -> v e+ -- svZipWith :: (e -> e -> e) -> v e -> v e -> v e --- instance IxContainer IM2 a where--- type Ix IM2 = (Int, Int)--- ixcIfilter f im2 = IM2 $ ifilterIM2 (curry f) (unIM2 im2) --- class Rank2 (c :: * -> *) a where--- type R2IxRow c :: *--- type R2IxCol c :: *--- type R2V c :: * -> *--- rank2Act :: c a -> R2V c a -> R2V c a--- extractRow :: c a -> R2IxRow c -> Maybe (R2V c a)--- extractCol :: c a -> R2IxCol c -> Maybe (R2V c a)--- -- extractRows :: c a -> [R2IxRow c] -> [R2V c]--- -- extractCols :: c a -> [R2IxCol c] -> [R2V c]- --- * SMatrix : sparse matrix types+-- * SparseMatrix --- class (IxContainer c a, Sparse c a, Additive c) => SMatrix c a where+class SpContainer m e => SparseMatrix m e where+ smFromVector :: LexOrd -> (Int, Int) -> V.Vector (IxRow, IxCol, e) -> m e+ -- smFromFoldableDense :: Foldable t => t e -> m e + smTranspose :: m e -> m e+ -- smExtractSubmatrix ::+ -- m e -> (IxRow, IxRow) -> (IxCol, IxCol) -> m e+ encodeIx :: m e -> LexOrd -> (IxRow, IxCol) -> LexIx+ decodeIx :: m e -> LexOrd -> LexIx -> (IxRow, IxCol)+++-- data RowsFirst = RowsFirst+-- data ColsFirst = ColsFirst++-- class SpContainer m e => SparseMatrix m o e where+-- smFromVector :: o -> (Int, Int) -> V.Vector (IxRow, IxCol, e) -> m e+-- -- smFromFoldableDense :: Foldable t => t e -> m e +-- smTranspose :: o -> m e -> m e+-- -- smExtractSubmatrix ::+-- -- m e -> (IxRow, IxRow) -> (IxCol, IxCol) -> m e+-- encodeIx :: m e -> o -> (IxRow, IxCol) -> LexIx+-- decodeIx :: m e -> o -> LexIx -> (IxRow, IxCol)+++++++-- * SparseMatVec++-- | Combining functions for relating (structurally) matrices and vectors, e.g. extracting/inserting rows/columns/submatrices++-- class (SparseMatrix m o e, SparseVector v e) => SparseMatVec m o v e where+-- smvInsertRow :: m e -> v e -> IxRow -> m e+-- smvInsertCol :: m e -> v e -> IxCol -> m e+-- smvExtractRow :: m e -> IxRow -> v e+-- smvExtractCol :: m e -> IxCol -> v e +++++++++-- * Utilities+++-- | Lift a real number onto the complex plane+toC :: Num a => a -> Complex a+toC r = r :+ 0++++++++++++++++-- -- | Instances for AdditiveGroup+-- instance Integral a => AdditiveGroup (Ratio a) where+-- {zero=0; (^+^) = (+); negated = negate}++-- instance (RealFloat v, AdditiveGroup v) => AdditiveGroup (Complex v) where+-- zero = zero :+ zero+-- (^+^) = (+)+-- negated = negate++-- -- | Standard instance for an applicative functor applied to a vector space.+-- instance AdditiveGroup v => AdditiveGroup (a -> v) where+-- zero = pure zero+-- (^+^) = liftA2 (^+^)+-- negated = fmap negated+++-- -- | Instances for VectorSpace+-- instance (RealFloat v, VectorSpace v) => VectorSpace (Complex v) where+-- type Scalar (Complex v) = Scalar v+-- s .* (u :+ v) = s .* u :+ s .* v++++-- #define ScalarType(t) \+-- instance AdditiveGroup (t) where {zero = 0; (^+^) = (+); negated = negate};\+-- instance VectorSpace (t) where {type Scalar (t) = (t); (.*) = (*) };\+-- instance Hilbert (t) where dot = (*)++-- ScalarType(Int)+-- ScalarType(Integer)+-- ScalarType(Float)+-- ScalarType(Double)+-- ScalarType(CSChar)+-- ScalarType(CInt)+-- ScalarType(CShort)+-- ScalarType(CLong)+-- ScalarType(CLLong)+-- ScalarType(CIntMax)+-- ScalarType(CFloat)+-- ScalarType(CDouble)
+ src/Numeric/LinearAlgebra/EigenSolvers/Experimental.hs view
@@ -0,0 +1,51 @@+{-# language TypeFamilies, FlexibleContexts #-}+module Numeric.LinearAlgebra.EigenSolvers.Experimental where++import Control.Exception.Common+import Control.Iterative+import Data.Sparse.Common++import Control.Monad.Catch+import Control.Monad.State.Strict++import Data.VectorSpace++-- | Golub-Kahan-Lanczos bidiagonalization (see "Restarted Lanczos Bidiagonalization for the SVD", SLEPc STR-8, http://slepc.upv.es/documentation/manual.htm )+gklBidiag aa q1nn | dim q1nn == n = return (pp, bb, qq)+ | otherwise = throwM (MatVecSizeMismatchException "hhBidiag" (dim aa) (dim q1nn))+ where+ (m,n) = (nrows aa, ncols aa)+ aat = transpose aa+ bb = fromListSM (n,n) bl+ (ql, _, pl, _, pp, bl, qq) = execState (modifyUntil tf bidiagStep) bidiagInit+ tf (_, _, _, i, _, _, _) = i == n + bidiagInit = (q2n, beta1, p1n, 1 :: Int, pp, bb', qq)+ where+ q1 = normalize2' q1nn+ + p1 = aa #> q1+ alpha1 = norm2' p1+ p1n = p1 ./ alpha1+ + q2 = (aat #> p1) ^-^ (alpha1 .* q1)+ beta1 = norm2' q2+ q2n = q2 ./ beta1+ + pp = insertCol (zeroSM m n) p1n 0+ qq = insertCol (zeroSM n n) q2n 0+ bb' = [(0, 0, alpha1)]+ bidiagStep (qj , betajm, pjm , j , pp, bb, qq ) =+ (qjp, betaj , pj, succ j, pp', bb', qq') where++ u = (aa #> qj) ^-^ (betajm .* pjm)+ alphaj = norm2' u+ pj = u ./ alphaj+ + v = (aat #> pj) ^-^ (alphaj .* qj)+ betaj = norm2' v+ qjp = v ./ betaj+ + pp' = insertCol pp pj j+ bb' = [(j-1, j, betaj),+ (j ,j, alphaj)] ++ bb+ qq' = insertCol qq qjp j
+ src/Numeric/LinearAlgebra/LinearSolvers/Experimental.hs view
@@ -0,0 +1,59 @@+{-# language TypeFamilies, FlexibleContexts #-}+module Numeric.LinearAlgebra.LinearSolvers.Experimental where++import Control.Exception.Common+import Control.Iterative+import Data.Sparse.Common++import Control.Monad.Catch++import Data.VectorSpace+++-- ** TFQMR++tfqmrInit aa b x0 = TFQMR x0 w0 u0 v0 d0 m tau0 theta0 eta0 rho0 alpha0 where+ n = dim b+ r0 = b ^-^ (aa #> x0) -- residual of initial guess solution+ w0 = r0+ u0 = r0+ v0 = aa #> u0+ d0 = zeroSV n+ r0hat = r0+ rho0 = r0hat `dot` r0+ alpha0 = rho0 / (v0 `dot` r0hat)+ m = 0+ tau0 = norm2' r0+ theta0 = 0+ eta0 = 0++tfqmrStep aa r0hat (TFQMR x w u v d m tau theta eta rho alpha) =+ TFQMR x1 w1 u1 v1 d1 (m+1) tau1 theta1 eta1 rho1 alpha1+ where+ w1 = w ^-^ (alpha .* (aa #> u))+ d1 = u ^+^ ((theta**2/alpha*eta) .* d)+ theta1 = norm2' w1 / tau+ c = recip $ sqrt (1 + theta1**2)+ tau1 = tau * theta1 * c+ eta1 = c**2 * alpha+ x1 = x^+^ (eta1 .* d1)+ (alpha1, u1, rho1, v1)+ | even m = let+ alpha' = rho / (v `dot` r0hat)+ u' = u ^-^ (alpha' .* v)+ in+ (alpha', u', rho, v)+ | otherwise = let+ rho' = w1 `dot` r0hat+ beta = rho'/rho+ u' = w1 ^+^ (beta .* u)+ v' = (aa #> u') ^+^ (beta .* (aa #> u ^+^ (beta .* v)) )+ in (alpha, u', rho', v')++data TFQMR a =+ TFQMR { _xTfq, _wTfq, _uTfq, _vTfq, _dTfq :: SpVector a,+ _mTfq :: Int,+ _tauTfq, _thetaTfq, _etaTfq, _rhoTfq, _alphaTfq :: a}+ deriving Eq+instance Show a => Show (TFQMR a) where+ show (TFQMR x _ _ _ _ _ _ _ _ _ _) = "x = " ++ show x ++ "\n"
src/Numeric/LinearAlgebra/Sparse.hs view
@@ -1,1051 +1,1019 @@-{-# LANGUAGE FlexibleContexts, TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-}--- {-# OPTIONS_GHC -O2 -rtsopts -with-rtsopts=-K32m -prof#-}-module Numeric.LinearAlgebra.Sparse- (- -- * Matrix factorizations- qr,- lu,- chol,- -- * Condition number- conditionNumberSM,- -- * Householder reflection- hhMat, hhRefl,- -- * Givens' rotation- givens,- -- * Arnoldi iteration- arnoldi, - -- * Eigensolvers- eigsQR, eigRayleigh,- -- * Linear solvers- -- ** Iterative methods- linSolve, LinSolveMethod(..), (<\>),- pinv,- -- ** Direct methods- luSolve,- -- * Preconditioners- ilu0, mSsor,- -- * Matrix partitioning- diagPartitions,- -- * Random arrays- randArray,- -- * Random matrices and vectors- randMat, randVec, - -- ** Sparse "- randSpMat, randSpVec,- -- * Sparsify data- sparsifySV,- -- * Iteration combinators- modifyInspectN, runAppendN', untilConverged,- diffSqL- )- where---import Data.Sparse.Common---import Control.Monad.Primitive-import Control.Monad (mapM_, forM_, replicateM)-import Control.Monad.State.Strict--- import Control.Monad.Writer--- import Control.Monad.Trans--- import Control.Monad.Trans.State (runStateT)--- import Control.Monad.Trans.Writer (runWriterT)--import qualified Data.IntMap.Strict as IM--- import Data.Utils.StrictFold (foldlStrict) -- hidden in `containers`--import qualified System.Random.MWC as MWC-import qualified System.Random.MWC.Distributions as MWC--import Data.Monoid-import qualified Data.Foldable as F-import qualified Data.Traversable as T---- import qualified Data.List as L-import Data.Maybe--import qualified Data.Vector as V---- --- * Sparsify : remove almost-0 elements (|x| < eps)--- | Sparsify an SpVector-sparsifySV :: Epsilon a => SpVector a -> SpVector a-sparsifySV = filterSV isNz-------- * Matrix condition number---- |uses the R matrix from the QR factorization-conditionNumberSM :: (Epsilon a, RealFloat a) => SpMatrix a -> a-conditionNumberSM m | nearZero lmin = error "Infinite condition number : rank-deficient system"- | otherwise = kappa where- kappa = lmax / lmin- (_, r) = qr m- u = extractDiagDense r -- FIXME : need to extract with default element 0 - lmax = abs (maximum u)- lmin = abs (minimum u)---------- * Householder transformation--hhMat :: Num a => a -> SpVector a -> SpMatrix a-hhMat beta x = eye n ^-^ scale beta (x >< x) where- n = dim x---{-| a vector `x` uniquely defines an orthogonal plane; the Householder operator reflects any point `v` with respect to this plane:- v' = (I - 2 x >< x) v--}-hhRefl :: Num a => SpVector a -> SpMatrix a-hhRefl = hhMat (fromInteger 2)-------------- * Givens rotation matrix---hypot :: Floating a => a -> a -> a-hypot x y = abs x * (sqrt (1 + y/x)**2)--sign :: (Ord a, Num a) => a -> a-sign x- | x > 0 = 1- | x == 0 = 0- | otherwise = -1 ---- | Givens coefficients (using stable algorithm shown in Anderson, Edward (4 December 2000). "Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem". LAPACK Working Note)-givensCoef :: (Ord a, Floating a) => a -> a -> (a, a, a)-givensCoef a b -- returns (c, s, r) where r = norm (a, b)- | b==0 = (sign a, 0, abs a)- | a==0 = (0, sign b, abs b)- | abs a > abs b = let t = b/a- u = sign a * abs ( sqrt (1 + t**2))- in (1/u, - t/u, a*u)- | otherwise = let t = a/b- u = sign b * abs ( sqrt (1 + t**2))- in (t/u, - 1/u, b*u)---{- |-Givens method, row version: choose other row index i' s.t. i' is :-* below the diagonal-* corresponding element is nonzero--QR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is-non-zero but A has zeros in row k for all columns less than j.--NB: the current version is quite inefficient in that:-1. the Givens' matrix `G_i` is different from Identity only in 4 entries-2. at each iteration `i` we multiply `G_i` by the previous partial result `M`. Since this corresponds to a rotation, and the `givensCoef` function already computes the value of the resulting non-zero component (output `r`), `G_i ## M` can be simplified by just changing two entries of `M` (i.e. zeroing one out and changing the other into `r`).--}-{-# inline givens #-}-givens :: (Floating a, Epsilon a, Ord a) => SpMatrix a -> IxRow -> IxCol -> SpMatrix a-givens mm i j - | isValidIxSM mm (i,j) && nrows mm >= ncols mm =- fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))- | otherwise = error "givens : indices out of bounds" - where- (c, s, _) = givensCoef a b- i' = head $ fromMaybe (error $ "givens: no compatible rows for entry " ++ show (i,j)) (candidateRows (immSM mm) i j)- a = mm @@ (i', j)- b = mm @@ (i, j) -- element to zero out---- |Returns a set of rows {k} that satisfy QR.C1-candidateRows :: IM.IntMap (IM.IntMap a) -> IxRow -> IxCol -> Maybe [IM.Key]-candidateRows mm i j | IM.null u = Nothing- | otherwise = Just (IM.keys u) where- u = IM.filterWithKey (\irow row -> irow /= i &&- firstNonZeroColumn row j) mm---- |Is the `k`th the first nonzero column in the row?-{-# inline firstNonZeroColumn #-}-firstNonZeroColumn :: IM.IntMap a -> IxRow -> Bool-firstNonZeroColumn mm k = isJust (IM.lookup k mm) &&- isNothing (IM.lookupLT k mm)---------- * QR decomposition----- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.-qr :: (Epsilon a, Ord a, Floating a) => SpMatrix a -> (SpMatrix a, SpMatrix a)-qr mm = (transposeSM qt, r) where- (qt, r, _) = execState (modifyUntil qf stepf) gminit- qf (_, _, iis) = null iis- stepf (qmatt, m, iis) = (qmatt', m', tail iis) where- (i, j) = head iis- g = givens m i j- qmatt' = g #~# qmatt -- update Q'- m' = g #~# m -- update R- gminit = (eye (nrows mm), mm, subdiagIndicesSM mm)-- - ----------- * Eigenvalue algorithms---- ** QR algorithm---- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm`, and returns a SpVector containing all eigenvalues-eigsQR :: (Epsilon a, Real a, Floating a) => Int -> SpMatrix a -> SpVector a-eigsQR nitermax m = extractDiagDense $ execState (convergtest eigsStep) m where- eigsStep m = r #~# q where (q, r) = qr m- convergtest g = modifyInspectN nitermax f g where- f [m1, m2] = let dm1 = extractDiagDense m1- dm2 = extractDiagDense m2- in nearZero $ norm2 (dm1 ^-^ dm2)--------- ** Rayleigh iteration---- | `eigsRayleigh n mm` performs `n` iterations of the Rayleigh algorithm on matrix `mm` and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair--- eigRayleigh :: Int -- max # iterations--- -> SpMatrix Double -- matrix--- -> (SpVector Double, Double) -- initial guess of (eigenvector, eigenvalue)--- -> (SpVector Double, Double) -- final estimate of (eigenvector, eigenvalue)-eigRayleigh nitermax m = execState (convergtest (rayleighStep m)) where- convergtest g = modifyInspectN nitermax f g where- f [(b1, _), (b2, _)] = nearZero $ norm2 (b2 ^-^ b1)- rayleighStep aa (b, mu) = (b', mu') where- ii = eye (nrows aa)- nom = (aa ^-^ (mu `matScale` ii)) <\> b- b' = normalize 2 nom- mu' = b' `dot` (aa #> b') / (b' `dot` b')------- * Householder vector ---- (Golub & Van Loan, Alg. 5.1.1, function `house`)-hhV :: (Epsilon a, Real a, Floating a) => SpVector a -> (SpVector a, a)-hhV x = (v, beta) where- n = dim x- tx = tailSV x- sigma = tx `dot` tx- vtemp = singletonSV 1 `concatSV` tx- (v, beta) | nearZero sigma = (vtemp, 0)- | otherwise = let mu = sqrt (headSV x**2 + sigma)- xh = headSV x- vh | xh <= 1 = xh - mu- | otherwise = - sigma / (xh + mu)- vnew = (1 / vh) .* insertSpVector 0 vh vtemp - in (vnew, 2 * xh**2 / (sigma + vh**2))-- ------ * Householder bidiagonalization--{- G & VL Alg. 5.4.2 -}--------- * SVD--{- Golub & Van Loan, sec 8.6.2 (p 452 segg.)--SVD of A, Golub-Kahan method--* reduce A to upper bidiagonal form B (Alg. 5.4.2, Householder bidiagonalization)-* compute SVD of B (implicit-shift QR step applied to B^T B, Alg. 8.3.2)---}----------- * Cholesky factorization---- ** Cholesky–Banachiewicz algorithm---- | Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A-chol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> SpMatrix a-chol aa = lfin where- (_, lfin) = execState (modifyUntil q cholUpd) cholInit- q (i, _) = i == nrows aa -- stopping criterion- cholInit = cholUpd (0, zeroSM n n) -- initialization- n = nrows aa- cholUpd (i, ll) = (i + 1, ll') where- ll' = cholDiagUpd (cholSDRowUpd ll) -- first upd subdiagonal entries in the row- cholSDRowUpd ll_ = insertRow ll_ lrs i where- lrs = fromListSV (i + 1) $ onRangeSparse (cholSubDiag ll i) [0 .. i-1]- cholDiagUpd ll_ = insertSpMatrix i i (cholDiag ll_ i) ll_ - cholSubDiag ll i j = 1/ljj*(aij - inn) where- ljj = ll@@(j, j)- aij = aa@@(i, j)- inn = contractSub ll ll i j (j - 1)- cholDiag ll i | i == 0 = sqrt aai- | i > 0 = sqrt $ aai - sum (fmap (**2) lrow)- | otherwise = error "cholDiag : index must be nonnegative" where- aai = aa@@(i,i)- lrow = ifilterSV (\j _ -> j < i) (extractRow ll i) -- sub-diagonal elems of L-------------------- * LU factorization--- ** Doolittle algorithm-{- Doolittle algorithm for factoring A' = P A, where P is a permutation matrix such that A' has a nonzero as its (0, 0) entry -}---- | Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = A-lu :: (Epsilon a, Fractional a, Real a) => SpMatrix a -> (SpMatrix a, SpMatrix a)-lu aa = (lf, ufin) where- (ixf, lf, uf) = execState (modifyUntil q luUpd) luInit- ufin = uUpdSparse (ixf, lf, uf) -- final U update- q (i, _, _) = i == (nrows aa - 1)- n = nrows aa- luInit = (1, l0, u0) where- l0 = insertCol (eye n) ((1/u00) .* extractSubCol aa 0 (1,n - 1)) 0 -- initial L- u0 = insertRow (zeroSM n n) (extractRow aa 0) 0 -- initial U- u00 = u0 @@ (0,0) -- make sure this is non-zero by applying permutation- luUpd (i, l, u) = (i + 1, l', u') where- u' = uUpdSparse (i, l, u) -- update U- l' = lUpdSparse (i, l, u') -- update L- uUpdSparse (ix, lmat, umat) = insertRow umat (fromListSV n us) ix where- us = onRangeSparse (solveForUij ix) [ix .. n - 1]- solveForUij i j = a - p where- a = aa @@! (i, j)- p = contractSub lmat umat i j (i - 1)- lUpdSparse (ix, lmat, umat) = insertCol lmat (fromListSV n ls) ix where- ls = onRangeSparse (`solveForLij` ix) [ix + 1 .. n - 1]- solveForLij i j- | isNz ujj = (a - p)/ujj- | otherwise =- error $ unwords ["solveForLij : U",- show (j ,j ),- "is close to 0. Permute rows in order to have a nonzero diagonal of U"]- where- a = aa @@! (i, j)- ujj = umat @@! (j , j) -- NB this must be /= 0- p = contractSub lmat umat i j (i - 1)------ | Apply a function over a range of integer indices, zip the result with it and filter out the almost-zero entries-onRangeSparse :: Epsilon b => (Int -> b) -> [Int] -> [(Int, b)]-onRangeSparse f ixs = filter (isNz . snd) $ zip ixs $ map f ixs----------------------- -- Produces the permutation matrix necessary to have a nonzero in position (iref, jref). This is used in the LU factorization--- permutAA :: Num b => IxRow -> IxCol -> SpMatrix a -> Maybe (SpMatrix b)--- permutAA iref jref (SM (nro,_) mm) --- | isJust (lookupIM2 iref jref mm) = Nothing -- eye nro--- | otherwise = Just $ permutationSM nro [head u] where--- u = IM.keys (ifilterIM2 ff mm)--- ff i j _ = i /= iref &&--- j == jref------------------- * Arnoldi iteration---- | Given a matrix A, a vector b and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q.--- At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.--arnoldi ::- (Epsilon a, Floating a, Eq a) => SpMatrix a -> SpVector a -> Int -> (SpMatrix a, SpMatrix a)-arnoldi aa b kn = (fromCols qvfin, fromListSM (nmax + 1, nmax) hhfin)- where- (qvfin, hhfin, nmax, _) = execState (modifyUntil tf arnoldiStep) arnInit - tf (_, _, ii, fbreak) = ii == kn || fbreak -- termination criterion- (m, n) = dim aa- arnInit = (qv1, hh1, 1, False) where- q0 = normalize 2 b -- starting basis vector- aq0 = aa #> q0 -- A q0- h11 = q0 `dot` aq0 - q1nn = (aq0 ^-^ (h11 .* q0))- hh1 = V.fromList [(0, 0, h11), (1, 0, h21)] where - h21 = norm 2 q1nn- q1 = normalize 2 q1nn -- q1 `dot` q0 ~ 0- qv1 = V.fromList [q0, q1]- arnoldiStep (qv, hh, i, _) = (qv', hh', i + 1, fb') where- qi = V.last qv- aqi = aa #> qi- hhcoli = fmap (`dot` aqi) qv -- H_{1, i}, H_{2, i}, .. , H_{m + 1, i}- zv = zeroSV m- qipnn =- aqi ^-^ (V.foldl' (^+^) zv (V.zipWith (.*) hhcoli qv)) -- unnormalized q_{i+1}- qipnorm = norm 2 qipnn -- normalization factor H_{i+1, i}- qip = normalize 2 qipnn -- q_{i + 1}- hh' = (V.++) hh (indexed2 $ V.snoc hhcoli qipnorm) where -- update H- indexed2 v = V.zip3 ii jj v- ii = V.fromList [0 .. n] -- nth col of upper Hessenberg has `n+1` nz- jj = V.replicate (n + 1) i -- `n+1` replicas of `i`- qv' = V.snoc qv qip -- append q_{i+1} to Krylov basis Q_i- fb' | nearZero qipnorm = True -- breakdown condition- | otherwise = False---------- * Preconditioning---- | Partition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts-diagPartitions :: SpMatrix a -> (SpMatrix a, SpMatrix a, SpMatrix a)-diagPartitions aa = (e,d,f) where- e = extractSubDiag aa- d = extractDiag aa- f = extractSuperDiag aa----- -- ** Jacobi preconditioner---- -- | Returns the reciprocal of the diagonal --- jacobiPreconditioner :: SpMatrix Double -> SpMatrix Double--- jacobiPreconditioner = reciprocal . extractDiag----- ** Incomplete LU---- | Used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2`-ilu0 :: (Epsilon a, Real a, Fractional a) => SpMatrix a -> (SpMatrix a, SpMatrix a)-ilu0 aa = (lh, uh) where- (l, u) = lu aa- lh = sparsifyLU l aa- uh = sparsifyLU u aa- sparsifyLU m m2 = ifilterSM f m where- f i j _ = isJust (lookupSM m2 i j)----- ** SSOR---- | `mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditioner-mSsor :: Fractional a => SpMatrix a -> a -> (SpMatrix a, SpMatrix a)-mSsor aa omega = (l, r) where- (e, d, f) = diagPartitions aa- n = nrows e- l = (eye n ^-^ scale omega e) ## reciprocal d- r = d ^-^ scale omega f ----------- * Linear solver, LU-based---- | Direct solver based on a triangular factorization of the system matrix.-luSolve ::- (Fractional a, Eq a, Epsilon a) => SpMatrix a -> SpMatrix a -> SpVector a -> SpVector a-luSolve ll uu b- | isLowerTriSM ll && isUpperTriSM uu = triUpperSolve uu (triLowerSolve ll b)- | otherwise = error "luSolve : factors must be triangular matrices" --triLowerSolve :: (Epsilon a, Fractional a) => SpMatrix a -> SpVector a -> SpVector a-triLowerSolve ll b = sparsifySV v where- (v, _) = execState (modifyUntil q lStep) lInit where- q (_, i) = i == dim b- lStep (ww, i) = (ww', i + 1) where- lii = ll @@ (i, i)- bi = b @@ i- wi = (bi - r)/lii where- r = extractSubRow ll i (0, i-1) `dot` takeSV i ww- ww' = insertSpVector i wi ww- lInit = (ww0, 1) where- l00 = ll @@ (0, 0)- b0 = b @@ 0- w0 = b0 / l00- ww0 = insertSpVector 0 w0 $ zeroSV (dim b) ---- | NB in the computation of `xi` we must rebalance the subrow indices because `dropSV` does that too, in order to take the inner product with consistent index pairs-triUpperSolve :: (Epsilon a, Fractional a) => SpMatrix a -> SpVector a -> SpVector a-triUpperSolve uu w = sparsifySV x where- (x, _) = execState (modifyUntil q uStep) uInit- q (_, i) = i == (- 1)- uStep (xx, i) = (xx', i - 1) where- uii = uu @@ (i, i)- wi = w @@ i- xi = (wi - r) / uii where- r = extractSubRow_RK uu i (i + 1, dim w - 1) `dot` dropSV (i + 1) xx- xx' = insertSpVector i xi xx- uInit = (xx0, i - 1) where- i = dim w - 1- u00 = uu @@ (i, i)- w0 = w @@ i- x0 = w0 / u00- xx0 = insertSpVector i x0 $ zeroSV (dim w)--------- * Iterative linear solvers----- ** GMRES---- | Given a linear system `A x = b` where `A` is an (m x m) real-valued matrix, the GMRES method finds an approximate solution `xhat` such that the Euclidean norm of the residual `A xhat - b` is minimized. `xhat` is spanned by the order-`n` Krylov subspace of (A, b).--- In this implementation:--- 1) the Arnoldi iteration is carried out until numerical breakdown (therefore yielding _at_most_ `m` Krylov basis vectors)--- 2) the resulting Hessenberg matrix is factorized in QR form--- 3) the Krylov-subspace solution `yhat` is found by backsubstitution--- 4) the approximate solution in the original space `xhat` is computed using the Krylov basis, `xhat = Q_n yhat`------ Many optimizations are possible, for example interleaving the QR factorization with the Arnoldi process (and employing an updating QR factorization which only requires one Givens' rotation at every update). --gmres :: (Epsilon a, Ord a, Floating a) => SpMatrix a -> SpVector a -> SpVector a-gmres aa b = qa' #> yhat where- m = ncols aa- (qa, ha) = arnoldi aa b m -- at most m steps of Arnoldi (aa, b)- -- b' = (transposeSe qa) #> b- b' = norm2 b .* (ei mp1 1) -- b rotated back to canonical basis by Q^T- where mp1 = nrows ha -- 1 + # Arnoldi iterations- (qh, rh) = qr ha -- QR fact.of H- yhat = triUpperSolve rh' rhs' where- rhs' = takeSV (dim b' - 1) (transposeSM qh #> b')- rh' = takeRows (nrows rh - 1) rh- qa' = takeCols (ncols qa - 1) qa---------- *** Left-preconditioning---------- ** CGNE--cgneStep :: SpMatrix Double -> CGNE -> CGNE-cgneStep aa (CGNE x r p) = CGNE x1 r1 p1 where- alphai = r `dot` r / (p `dot` p)- x1 = x ^+^ (alphai .* p)- r1 = r ^-^ (alphai .* (aa #> p))- beta = r1 `dot` r1 / (r `dot` r)- p1 = transposeSM aa #> r ^+^ (beta .* p)--data CGNE =- CGNE {_xCgne , _rCgne, _pCgne :: SpVector Double} deriving Eq-instance Show CGNE where- show (CGNE x r p) = "x = " ++ show x ++ "\n" ++- "r = " ++ show r ++ "\n" ++- "p = " ++ show p ++ "\n"--cgne :: SpMatrix Double -> SpVector Double -> SpVector Double -> CGNE-cgne aa b x0 = execState (untilConverged _xCgne (cgneStep aa)) cgneInit where- r0 = b ^-^ (aa #> x0) -- residual of initial guess solution- p0 = transposeSM aa #> r0- cgneInit = CGNE x0 r0 p0----- ** TFQMR---- | one step of TFQMR-tfqmrStep :: SpMatrix Double -> SpVector Double -> TFQMR -> TFQMR-tfqmrStep aa r0hat (TFQMR x w u v d m tau theta eta rho alpha) =- TFQMR x1 w1 u1 v1 d1 (m+1) tau1 theta1 eta1 rho1 alpha1- where- -- alpham = alpha- w1 = w ^-^ (alpha .* (aa #> u))- d1 = u ^+^ ((theta**2/alpha*eta) .* d)- theta1 = norm2 w1 / tau- c = recip $ sqrt (1 + theta1**2)- tau1 = tau * theta1 * c- eta1 = c**2 * alpha- x1 = x^+^ (eta1 .* d1)- (alpha1, u1, rho1, v1)- | even m = let- alpha' = rho / (v `dot` r0hat)- u' = u ^-^ (alpha' .* v)- in- (alpha', u', rho, v)- | otherwise = let- rho' = w1 `dot` r0hat- beta = rho'/rho- u' = w1 ^+^ (beta .* u)- v' = (aa #> u') ^+^ (beta .* (aa #> u ^+^ (beta .* v)) )- in (alpha, u', rho', v')--tfqmr :: SpMatrix Double -> SpVector Double -> SpVector Double -> TFQMR -tfqmr aa b x0 = execState (untilConverged _xTfq (tfqmrStep aa r0)) tfqmrInit where- n = dim b- r0 = b ^-^ (aa #> x0) -- residual of initial guess solution- w0 = r0- u0 = r0- v0 = aa #> u0- d0 = zeroSV n- r0hat = r0- rho0 = r0hat `dot` r0- alpha0 = rho0 / (v0 `dot` r0hat)- m = 0- tau0 = norm2 r0- theta0 = 0- eta0 = 0- tfqmrInit = TFQMR x0 w0 u0 v0 d0 m tau0 theta0 eta0 rho0 alpha0--data TFQMR =- TFQMR { _xTfq, _wTfq, _uTfq, _vTfq, _dTfq :: SpVector Double,- _mTfq :: Int,- _tauTfq, _thetaTfq, _etaTfq, _rhoTfq, _alphaTfq :: Double}- deriving Eq-instance Show TFQMR where- show (TFQMR x _ _ _ _ _ _ _ _ _ _) = "x = " ++ show x ++ "\n"----- ** BCG---- | one step of BCG-bcgStep :: SpMatrix Double -> BCG -> BCG-bcgStep aa (BCG x r rhat p phat) = BCG x1 r1 rhat1 p1 phat1 where- aap = aa #> p- alpha = (r `dot` rhat) / (aap `dot` phat)- x1 = x ^+^ (alpha .* p)- r1 = r ^-^ (alpha .* aap)- rhat1 = rhat ^-^ (alpha .* (transposeSM aa #> phat))- beta = (r1 `dot` rhat1) / (r `dot` rhat)- p1 = r1 ^+^ (beta .* p)- phat1 = rhat1 ^+^ (beta .* phat)--data BCG =- BCG { _xBcg, _rBcg, _rHatBcg, _pBcg, _pHatBcg :: SpVector Double } deriving Eq--bcg :: SpMatrix Double -> SpVector Double -> SpVector Double -> BCG-bcg aa b x0 = execState (untilConverged _xBcg (bcgStep aa)) bcgInit where- r0 = b ^-^ (aa #> x0) -- residual of initial guess solution- r0hat = r0- p0 = r0- p0hat = r0- bcgInit = BCG x0 r0 r0hat p0 p0hat--instance Show BCG where- show (BCG x r rhat p phat) = "x = " ++ show x ++ "\n" ++- "r = " ++ show r ++ "\n" ++- "r_hat = " ++ show rhat ++ "\n" ++- "p = " ++ show p ++ "\n" ++- "p_hat = " ++ show phat ++ "\n"----- ** CGS---- | one step of CGS-cgsStep :: SpMatrix Double -> SpVector Double -> CGS -> CGS-cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1- where- aap = aa #> p- alphaj = (r `dot` rhat) / (aap `dot` rhat)- q = u ^-^ (alphaj .* aap)- xj1 = x ^+^ (alphaj .* (u ^+^ q)) -- updated solution- rj1 = r ^-^ (alphaj .* (aa #> (u ^+^ q))) -- updated residual- betaj = (rj1 `dot` rhat) / (r `dot` rhat)- uj1 = rj1 ^+^ (betaj .* q)- pj1 = uj1 ^+^ (betaj .* (q ^+^ (betaj .* p)))--data CGS = CGS { _x, _r, _p, _u :: SpVector Double} deriving Eq---- | iterate solver until convergence or until max # of iterations is reached-cgs ::- SpMatrix Double ->- SpVector Double ->- SpVector Double ->- SpVector Double ->- CGS-cgs aa b x0 rhat =- execState (untilConverged _x (cgsStep aa rhat)) cgsInit where- r0 = b ^-^ (aa #> x0) -- residual of initial guess solution- p0 = r0- u0 = r0- cgsInit = CGS x0 r0 p0 u0---instance Show CGS where- show (CGS x r p u) = "x = " ++ show x ++ "\n" ++- "r = " ++ show r ++ "\n" ++- "p = " ++ show p ++ "\n" ++- "u = " ++ show u ++ "\n"-------- ** BiCGSTAB---- _aa :: SpMatrix Double, -- matrix--- _b :: SpVector Double, -- rhs--- _r0 :: SpVector Double, -- initial residual--- _r0hat :: SpVector Double, -- candidate solution: r0hat `dot` r0 >= 0---- | one step of BiCGSTAB-bicgstabStep :: SpMatrix Double -> SpVector Double -> BICGSTAB -> BICGSTAB-bicgstabStep aa r0hat (BICGSTAB x r p) = BICGSTAB xj1 rj1 pj1 where- aap = aa #> p- alphaj = (r `dot` r0hat) / (aap `dot` r0hat)- sj = r ^-^ (alphaj .* aap)- aasj = aa #> sj- omegaj = (aasj `dot` sj) / (aasj `dot` aasj)- xj1 = x ^+^ (alphaj .* p) ^+^ (omegaj .* sj) -- updated solution- rj1 = sj ^-^ (omegaj .* aasj)- betaj = (rj1 `dot` r0hat)/(r `dot` r0hat) * alphaj / omegaj- pj1 = rj1 ^+^ (betaj .* (p ^-^ (omegaj .* aap)))--data BICGSTAB =- BICGSTAB { _xBicgstab, _rBicgstab, _pBicgstab :: SpVector Double} deriving Eq---- | iterate solver until convergence or until max # of iterations is reached-bicgstab- :: SpMatrix Double- -> SpVector Double- -> SpVector Double- -> SpVector Double- -> BICGSTAB-bicgstab aa b x0 r0hat =- execState (untilConverged _xBicgstab (bicgstabStep aa r0hat)) bicgsInit where- r0 = b ^-^ (aa #> x0) -- residual of initial guess solution- p0 = r0- bicgsInit = BICGSTAB x0 r0 p0--instance Show BICGSTAB where- show (BICGSTAB x r p) = "x = " ++ show x ++ "\n" ++- "r = " ++ show r ++ "\n" ++- "p = " ++ show p ++ "\n"-------- * Moore-Penrose pseudoinverse--- | Least-squares approximation of a rectangular system of equaitons. Uses <\\> for the linear solve-pinv :: SpMatrix Double -> SpVector Double -> SpVector Double-pinv aa b = aa #^# aa <\> atb where- atb = transposeSM aa #> b---------- * Linear solver interface--data LinSolveMethod = GMRES_ | CGNE_ | TFQMR_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) ---- -- | Linear solve with _random_ starting vector--- linSolveM ::--- PrimMonad m =>--- LinSolveMethod -> SpMatrix Double -> SpVector Double -> m (SpVector Double)--- linSolveM method aa b = do--- let (m, n) = dim aa--- nb = dim b--- if n /= nb then error "linSolve : operand dimensions mismatch" else do--- x0 <- randVec nb--- case method of CGS_ -> return $ _xBicgstab (bicgstab aa b x0 x0)--- BICGSTAB_ -> return $ _x (cgs aa b x0 x0)---- | Linear solve with _deterministic_ starting vector (every component at 0.1) -linSolve ::- LinSolveMethod -> SpMatrix Double -> SpVector Double -> SpVector Double-linSolve method aa b- | n /= nb = error "linSolve : operand dimensions mismatch"- | otherwise = solve aa b where- solve aa' b' | isDiagonalSM aa' = reciprocal aa' #> b' -- diagonal solve is easy- | otherwise = solveWith aa' b' - solveWith aa' b' = case method of- GMRES_ -> gmres aa' b'- CGNE_ -> _xCgne (cgne aa' b' x0)- TFQMR_ -> _xTfq (tfqmr aa' b' x0)- BCG_ -> _xBcg (bcg aa' b' x0)- BICGSTAB_ -> _xBicgstab (bicgstab aa' b' x0 x0)- CGS_ -> _x (cgs aa' b' x0 x0)- x0 = mkSpVectorD n $ replicate n 0.1 - (m, n) = dim aa- nb = dim b---- | <\\> : linSolve using the BiCGSTAB method as default-(<\>) :: SpMatrix Double -> SpVector Double -> SpVector Double -(<\>) = linSolve GMRES_ - ------------ | TODO : if system is poorly conditioned, is it better to warn the user or just switch solvers (e.g. via the pseudoinverse) ?---- linSolveQR aa b init f1 stepf--- | isInfinite k = do--- tell "linSolveQR : rank-deficient system"--- | otherwise = do--- solv aa b init--- where--- (q, r) = qr aa--- k = conditionNumberSM r--- solv aa b init = execState (untilConverged f1 stepf) init----------------------- * Control primitives for bounded iteration with convergence check---- | transform state until a condition is met-modifyUntil :: MonadState s m => (s -> Bool) -> (s -> s) -> m s-modifyUntil q f = do- x <- get- let y = f x- put y- if q y then return y- else modifyUntil q f ---- | Keep a moving window buffer (length 2) of state `x` to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached -loopUntilAcc :: Int -> ([t] -> Bool) -> (t -> t) -> t -> t-loopUntilAcc nitermax q f x = go 0 [] x where- go i ll xx | length ll < 2 = go (i + 1) (y : ll) y - | otherwise = if q ll || i == nitermax- then xx- else go (i + 1) (take 2 $ y:ll) y- where y = f xx---- | Keep a moving window buffer (length 2) of state `x` to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached (i.e. same thing as `loopUntilAcc` but this one runs in the State monad)-modifyInspectN ::- MonadState s m =>- Int -> -- iteration budget- ([s] -> Bool) -> -- convergence criterion- (s -> s) -> -- state stepping function- m s-modifyInspectN nitermax q f - | nitermax > 0 = go 0 []- | otherwise = error "modifyInspectN : n must be > 0" where- go i ll = do- x <- get- let y = f x- if length ll < 2- then do put y- go (i + 1) (y : ll)- else if q ll || i == nitermax- then do put y- return y- else do put y- go (i + 1) (take 2 $ y : ll)----- helper functions for estimating convergence-meanl :: (Foldable t, Fractional a) => t a -> a-meanl xx = 1/fromIntegral (length xx) * sum xx--norm2l :: (Foldable t, Functor t, Floating a) => t a -> a-norm2l xx = sqrt $ sum (fmap (**2) xx)--diffSqL :: Floating a => [a] -> a-diffSqL xx = (x1 - x2)**2 where [x1, x2] = [head xx, xx!!1]---------- | iterate until convergence is verified or we run out of a fixed iteration budget-untilConverged :: MonadState a m => (a -> SpVector Double) -> (a -> a) -> m a-untilConverged fproj = modifyInspectN 200 (normDiffConverged fproj)---- | convergence check (FIXME)-normDiffConverged :: (Foldable t, Functor t) =>- (a -> SpVector Double) -> t a -> Bool-normDiffConverged fp xx = nearZero $ normSq (foldrMap fp (^-^) (zeroSV 0) xx)--- ------ | run `niter` iterations and append the state `x` to a list `xs`, stop when either the `xs` satisfies a predicate `q` or when the counter reaches 0-runAppendN :: ([t] -> Bool) -> (t -> t) -> Int -> t -> [t]-runAppendN qq ff niter x0 | niter<0 = error "runAppendN : niter must be > 0"- | otherwise = go qq ff niter x0 [] where- go q f n z xs = - let x = f z in- if n <= 0 || q xs then xs- else go q f (n-1) x (x : xs)---- | ", NO convergence check -runAppendN' :: (t -> t) -> Int -> t -> [t]-runAppendN' ff niter x0 | niter<0 = error "runAppendN : niter must be > 0"- | otherwise = go ff niter x0 [] where- go f n z xs = - let x = f z in- if n <= 0 then xs- else go f (n-1) x (x : xs)-- --------------- * Random arrays--randArray :: PrimMonad m => Int -> Double -> Double -> m [Double]-randArray n mu sig = do- g <- MWC.create- replicateM n (MWC.normal mu sig g)- ------ * Random matrices and vectors---- |Dense SpMatrix-randMat :: PrimMonad m => Int -> m (SpMatrix Double)-randMat n = do- g <- MWC.create- aav <- replicateM (n^2) (MWC.normal 0 1 g)- let ii_ = [0 .. n-1]- (ix_,iy_) = unzip $ concatMap (zip ii_ . replicate n) ii_- return $ fromListSM (n,n) $ zip3 ix_ iy_ aav---- | Dense SpVector -randVec :: PrimMonad m => Int -> m (SpVector Double)-randVec n = do- g <- MWC.create- bv <- replicateM n (MWC.normal 0 1 g)- let ii_ = [0..n-1]- return $ fromListSV n $ zip ii_ bv------ | Sparse SpMatrix-randSpMat :: Int -> Int -> IO (SpMatrix Double)-randSpMat n nsp | nsp > n^2 = error "randSpMat : nsp must be < n^2 "- | otherwise = do- g <- MWC.create- aav <- replicateM nsp (MWC.normal 0 1 g)- ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)- jj <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)- return $ fromListSM (n,n) $ zip3 ii jj aav---- | Sparse SpVector-randSpVec :: Int -> Int -> IO (SpVector Double)-randSpVec n nsp | nsp > n = error "randSpVec : nsp must be < n"- | otherwise = do- g <- MWC.create- aav <- replicateM nsp (MWC.normal 0 1 g)- ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)- return $ fromListSV n $ zip ii aav-+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE FlexibleContexts, TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-}+{-# language ApplicativeDo #-}+-- {-# OPTIONS_GHC -O2 -rtsopts -with-rtsopts=-K32m -prof#-}+{-|++This module exposes the high-level functionality of the library.++-}++module Numeric.LinearAlgebra.Sparse+ (+ -- * Linear solvers+ -- ** Iterative methods+ linSolve0, LinSolveMethod(..), (<\>),+ -- ** Moore-Penrose pseudoinverse+ pinv,+ -- ** Preconditioners+ jacobiPre, ilu0, mSsor, + -- ** Direct methods+ luSolve,+ -- *** Forward substitution+ triLowerSolve,+ -- *** Backward substitution+ triUpperSolve,+ -- * Eigensolvers+ eigsQR,+ eigRayleigh,+ -- * Matrix factorization algorithms+ -- ** QR+ qr,+ -- ** LU+ lu,+ -- ** Cholesky+ chol, + -- ** Arnoldi iteration+ arnoldi, + -- * Matrix partitioning+ diagPartitions,+ -- * Utilities+ -- ** Givens' rotation+ givens,+ -- ** Condition number+ conditionNumberSM,+ -- ** Householder reflection+ hhMat, hhRefl,+ -- -- * Householder bidiagonalization + -- -- * Random arrays+ -- randArray,+ -- -- * Random matrices and vectors+ -- randMat, randVec, + -- -- ** Sparse "+ -- randSpMat, randSpVec,+ -- * From/to SpVector+ fromListSV, toListSV,+ -- * From/to SpMatrix+ fromListSM, toListSM,+ -- * Iteration combinators+ untilConvergedG0, untilConvergedG, untilConvergedGM,+ modifyInspectGuarded, modifyInspectGuardedM, IterationConfig (..),+ modifyUntil, modifyUntilM+ )+ where+++import Control.Exception.Common+import Control.Iterative+import Data.Sparse.Common++import Control.Monad.Catch+import Data.Typeable++-- import Control.Applicative ((<|>))++-- import Control.Monad (replicateM)+import Control.Monad.State.Strict+-- import Control.Monad.Writer+-- import Control.Monad.Trans.Class+import qualified Control.Monad.Trans.State as MTS -- (runStateT)+-- import Control.Monad.Trans.Writer (runWriterT)+import Data.Complex++import Data.VectorSpace hiding (magnitude)++import qualified Data.Sparse.Internal.IntM as I+-- import Data.Utils.StrictFold (foldlStrict) -- hidden in `containers`++-- import qualified System.Random.MWC as MWC+-- import qualified System.Random.MWC.Distributions as MWC++-- import Data.Monoid+-- import qualified Data.Foldable as F+-- import qualified Data.Traversable as T++-- import qualified Data.List as L++import Data.Maybe++import qualified Data.Vector as V++++-- | A lumped constraint for numerical types+type Num' x = (Epsilon x, Elt x, Show x, Ord x)++++++-- * Matrix condition number++-- |uses the R matrix from the QR factorization+conditionNumberSM :: (MonadThrow m, MatrixRing (SpMatrix a), Num' a, Typeable a) =>+ SpMatrix a -> m a+conditionNumberSM m = do+ (_, r) <- qr m+ let+ u = extractDiagDense r + lmax = abs (maximum u)+ lmin = abs (minimum u)+ kappa = lmax / lmin + if nearZero lmin+ then throwM (HugeConditionNumber "conditionNumberSM" kappa)+ else return kappa + ++++++-- * Householder transformation++hhMat :: Num a => a -> SpVector a -> SpMatrix a+hhMat beta x = eye n ^-^ beta `scale` (x >< x) where+ n = dim x+++-- | Householder reflection: a vector `x` uniquely defines an orthogonal plane; the Householder operator reflects any point `v` with respect to this plane:+-- v' = (I - 2 x >< x) v+hhRefl :: Num a => SpVector a -> SpMatrix a+hhRefl = hhMat (fromInteger 2)++++++++++++-- * Givens rotation matrix++-- -- -- | Givens coefficients (using stable algorithm shown in Anderson, Edward (4 December 2000). "Discontinuous Plane Rotations and the Symmetric Eigenvalue Problem". LAPACK Working Note)+-- -- -- givensCoef0 :: (Ord b, Floating a, Eq a) => (a -> b) -> a -> a -> (a, a, a)+-- -- givensCoef0+-- -- :: (Signum t, Ord a, Floating t, Eq t) =>+-- -- (t -> a) -> t -> t -> (t, t, t)+-- givensCoef0 ff a b -- returns (c, s, r) where r = norm (a, b)+-- | b==0 = (signum' a, 0, abs a)+-- | a==0 = (0, signum' b, abs b)+-- | ff a > ff b = let t = b/a+-- u = signum' a * abs ( sqrt (1 + t**2))+-- in (1/u, - t/u, a*u)+-- | otherwise = let t = a/b+-- u = signum' b * abs ( sqrt (1 + t**2))+-- in (t/u, - 1/u, b*u)++-- -- | Givens coefficients, real-valued+-- givensCoef :: (Floating a, Ord a) => a -> a -> (a, a, a)+-- givensCoef = givensCoef0 abs++-- -- | Givens coefficients, complex-valued+-- givensCoefC :: RealFloat a => Complex a -> Complex a -> (Complex a, Complex a, Complex a)+-- givensCoefC = givensCoef0 magnitude++++ +{- |+Givens method, row version: choose other row index i' s.t. i' is :+* below the diagonal+* corresponding element is nonzero++QR.C1 ) To zero out entry A(i, j) we must find row k such that A(k, j) is+non-zero but A has zeros in row k for all columns less than j.++NB: the current version is quite inefficient in that:+1. the Givens' matrix `G_i` is different from Identity only in 4 entries+2. at each iteration `i` we multiply `G_i` by the previous partial result `M`. Since this corresponds to a rotation, and the `givensCoef` function already computes the value of the resulting non-zero component (output `r`), `G_i ## M` can be simplified by just changing two entries of `M` (i.e. zeroing one out and changing the other into `r`).+-}+{-# inline givens #-}+givens :: (Elt a, MonadThrow m) => SpMatrix a -> Int -> Int -> m (SpMatrix a)+givens aa i j + | isValidIxSM aa (i,j) && nrows aa >= ncols aa = do+ i' <- candidateRows' (immSM aa) i j+ return $ givensMat aa i i' j+ | otherwise = throwM (OOBIxsError "Givens" [i, j])+ where+ givensMat mm i i' j =+ fromListSM'+ [(i,i, c), (j,j, conj c), (j,i, - conj s), (i,j, s)]+ (eye (nrows mm))+ where+ (c, s, _) = givensCoef a b+ a = mm @@ (i', j)+ b = mm @@ (i, j) -- element to zero out+ candidateRows' mm i j | null u = throwM (OOBNoCompatRows "Givens" (i,j))+ | otherwise = return $ head (I.keys u) where+ u = I.filterWithKey (\irow row -> irow /= i &&+ firstNZColumn row j) mm+ firstNZColumn m k = isJust (I.lookup k m) &&+ isNothing (I.lookupLT k m)++rotMat :: Elt e => e -> e -> IxRow -> IxRow -> Int -> SpMatrix e+rotMat a b i j n =+ fromListSM' [(i,i, c), (j,j, conj c), (j,i, - conj s), (i,j, s)] (eye n)+ where+ (c, s, _) = givensCoef a b+ ++-- | Givens coefficients and norm of associated vector+givensCoef :: Elt t => t -> t -> (t, t, t)+givensCoef a b = (c0/r, s0/r, r) where+ c0 = conj a+ s0 = - conj b+ r = hypot c0 s0 + hypot x y = abs x * sqrt (1 + (y/x)**2)++{-++ ( c s )+G =( )+ ( -s* c*)++-}+++ +++++-- * QR decomposition++-- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, where Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements.+{-# inline qr #-}+qr :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+ SpMatrix a -> m (SpMatrix a, SpMatrix a)+qr mm = do + (qt, r, _) <- MTS.execStateT (modifyUntilM haltf qrstepf) gminit+ return (transpose qt, r) + where+ gminit = (eye (nrows mm), mm, subdiagIndicesSM mm)+ haltf (_, _, iis) = null iis+ qrstepf (qmatt, m, iis) = do+ let (i, j) = head iis+ g <- givens m i j+ let+ qmatt' = g #~# qmatt -- update Q'+ m' = g #~# m -- update R+ return (qmatt', m', tail iis)++ +++++++++-- * Eigenvalue algorithms++-- ** QR algorithm++-- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm`, and returns a SpVector containing all eigenvalues+eigsQR :: (MonadThrow m, MonadIO m, Elt a, Normed (SpVector a), MatrixRing (SpMatrix a), Epsilon a, Typeable (Magnitude (SpVector a)), Typeable a, Show a) =>+ Int+ -> Bool -- ^ Print debug information + -> SpMatrix a -- ^ Operand matrix+ -> m (SpVector a) -- ^ Eigenvalues+eigsQR nitermax debq m = pf <$> untilConvergedGM "eigsQR" c (const True) stepf m+ where+ pf = extractDiagDense+ c = IterConf nitermax debq pf prd+ stepf mm = do+ (q, _) <- qr mm+ return $ q #~^# (m ## q) -- r #~# q++++ +++++-- ** Rayleigh iteration++-- | `eigsRayleigh n mm` performs `n` iterations of the Rayleigh algorithm on matrix `mm` and returns the eigenpair closest to the initialization. It displays cubic-order convergence, but it also requires an educated guess on the initial eigenpair.++eigRayleigh nitermax debq prntf m = untilConvergedGM "eigRayleigh" config (const True) (rayStep m)+ where+ ii = eye (nrows m)+ config = IterConf nitermax debq fst prntf+ rayStep aa (b, mu) = do+ nom <- (m ^-^ (mu `matScale` ii)) <\> b+ let b' = normalize2' nom+ mu' = (b' <.> (aa #> b')) / (b' <.> b')+ return (b', mu')+++-- | `eigArnoldi n aa b` computes at most n iterations of the Arnoldi algorithm to find a Krylov subspace of (A, b), along with a Hessenberg matrix of coefficients H. After that, it computes the QR decomposition of H, and the eigenvalues of A are listed on the diagonal of the R factor.+eigArnoldi :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,+ Elt t, Normed (SpVector t), MatrixRing (SpMatrix t),+ LinearVectorSpace (SpVector t), Epsilon t, MonadThrow m) =>+ Int+ -> SpMatrix t+ -> SpVector t+ -> m (SpMatrix t, SpMatrix t, SpVector t) -- ^ Q, O, R+eigArnoldi nitermax aa b = do+ (q, h) <- arnoldi aa b nitermax+ (o, r) <- qr h+ return (q, o, extractDiagDense r)+ ++++-- * Householder vector ++-- (Golub & Van Loan, Alg. 5.1.1, function `house`)+hhV :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t), Epsilon t) =>+ SpVector t -> (SpVector t, t)+hhV x = (v, beta) where+ tx = tailSV x+ sigma = tx <.> tx + vtemp = singletonSV 1 `concatSV` tx+ (v, beta) | nearZero sigma = (vtemp, 0)+ | otherwise = let mu = sqrt (headSV x**2 + sigma)+ xh = headSV x+ vh | mag xh <= 1 = xh - mu+ | otherwise = - sigma / (xh + mu)+ vnew = (1 / vh) `scale` insertSpVector 0 vh vtemp + in (vnew, 2 * xh**2 / (sigma + vh**2))++ ++++-- -- * Bidiagonalization++++-- -- * SVD++{- Golub & Van Loan, sec 8.6.2 (p 452 segg.)++SVD of A, Golub-Kahan method++* reduce A to upper bidiagonal form B (Alg. 5.4.2, Householder bidiagonalization)+* compute SVD of B (implicit-shift QR step applied to B^T B, Alg. 8.3.2)++-}+++++++++-- * Cholesky factorization++-- | Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A . This is an implementation of the Cholesky–Banachiewicz algorithm, i.e. proceeding row by row from the upper-left corner.+chol :: (Elt a, Epsilon a, MonadThrow m) =>+ SpMatrix a+ -> m (SpMatrix a) -- ^ L+chol aa = do+ let n = nrows aa+ q (i, _) = i == n+ l0 <- cholUpd aa (0, zeroSM n n)+ (_, lfin) <- MTS.execStateT (modifyUntilM q (cholUpd aa)) l0+ return lfin+ where+ oops i = throwM (NeedsPivoting "chol" (unwords ["L", show (i,i)]) :: MatrixException Double)+ cholUpd aa (i, ll) = do + sd <- cholSDRowUpd aa ll i -- update subdiagonal entries+ ll' <- cholDiagUpd aa sd i -- update diagonal entries+ return (i + 1, ll')+ cholSDRowUpd aa ll i = do + lrs <- fromListSV (i + 1) <$> onRangeSparseM cholSubDiag [0 .. i-1]+ return $ insertRow ll lrs i where+ cholSubDiag j | isNz ljj = return $ 1/ljj*(aij - inn)+ | otherwise = oops j+ where+ ljj = ll @@! (j, j)+ aij = aa @@! (i, j)+ inn = contractSub ll ll i j (j - 1)+ cholDiagUpd aa ll i = do+ cd <- cholDiag + return $ insertSpMatrix i i cd ll where+ cholDiag | i == 0 = sqrt <$> aai+ | otherwise = do+ a <- aai+ let l = sum (fmap (**2) lrow)+ return $ sqrt (a - l)+ where+ lrow = ifilterSV (\j _ -> j < i) (extractRow ll i) -- sub-diagonal elems of L+ aai | isNz aaii = return aaii+ | otherwise = oops i+ where+ aaii = aa @@! (i,i)+ +++++++++++++++++-- * LU factorization+-- ** Doolittle algorithm+{- Doolittle algorithm for factoring A' = P A, where P is a permutation matrix such that A' has a nonzero as its (0, 0) entry -}++-- | Given a matrix A, returns a pair of matrices (L, U) where L is lower triangular and U is upper triangular such that L U = A+lu :: (Scalar (SpVector t) ~ t, Elt t, VectorSpace (SpVector t), Epsilon t,+ MonadThrow m) =>+ SpMatrix t+ -> m (SpMatrix t, SpMatrix t) -- ^ L, U+lu aa = do+ let oops j = throwM (NeedsPivoting "solveForLij" ("U" ++ show (j, j)) :: MatrixException Double)+ n = nrows aa+ q (i, _, _) = i == n - 1+ luInit | isNz u00 = return (1, l0, u0)+ | otherwise = oops (0 :: Int)+ where+ l0 = insertCol (eye n) ((extractSubCol aa 0 (1, n-1)) ./ u00 ) 0+ u0 = insertRow (zeroSM n n) (extractRow aa 0) 0 -- initial U+ u00 = u0 @@! (0,0) -- make sure this is non-zero by applying permutation+ luUpd (i, l, u) = do -- (i + 1, l', u') + u' <- uUpd aa n (i, l, u) -- update U+ l' <- lUpd (i, l, u') -- update L+ return (i + 1, l', u')+ uUpd aa n (ix, lmat, umat) = do+ let us = onRangeSparse (solveForUij ix) [ix .. n - 1]+ solveForUij i j = a - p where+ a = aa @@! (i, j)+ p = contractSub lmat umat i j (i - 1)+ return $ insertRow umat (fromListSV n us) ix+ lUpd (ix, lmat, umat) = do -- insertCol lmat (fromListSV n ls) ix+ ls <- lsm+ return $ insertCol lmat (fromListSV n ls) ix+ where+ lsm = onRangeSparseM (`solveForLij` ix) [ix + 1 .. n - 1]+ solveForLij i j+ | isNz ujj = return $ (a - p)/ujj+ | otherwise = oops j+ where+ a = aa @@! (i, j)+ ujj = umat @@! (j , j) -- NB this must be /= 0+ p = contractSub lmat umat i j (i - 1) + s0 <- luInit+ (ixf, lf, uf) <- MTS.execStateT (modifyUntilM q luUpd) s0+ ufin <- uUpd aa n (ixf, lf, uf) -- final U update+ return (lf, ufin)+ +++++++ ++++++++-- -- Produces the permutation matrix necessary to have a nonzero in position (iref, jref). This is used in the LU factorization+-- permutAA :: Num b => IxRow -> IxCol -> SpMatrix a -> Maybe (SpMatrix b)+-- permutAA iref jref (SM (nro,_) mm) +-- | isJust (lookupIM2 iref jref mm) = Nothing -- eye nro+-- | otherwise = Just $ permutationSM nro [head u] where+-- u = IM.keys (ifilterIM2 ff mm)+-- ff i j _ = i /= iref &&+-- j == jref+++++++++++++++++-- * Arnoldi iteration++-- | Given a matrix A, a vector b and a positive integer `n`, this procedure finds the basis of an order `n` Krylov subspace (as the columns of matrix Q), along with an upper Hessenberg matrix H, such that A = Q^T H Q.+-- At the i`th iteration, it finds (i + 1) coefficients (the i`th column of the Hessenberg matrix H) and the (i + 1)`th Krylov vector.+arnoldi :: (MatrixType (SpVector a) ~ SpMatrix a, V (SpVector a) ,+ Scalar (SpVector a) ~ a, Epsilon a, MonadThrow m) =>+ SpMatrix a -- ^ System matrix+ -> SpVector a -- ^ r.h.s.+ -> Int -- ^ Max. # of iterations+ -> m (SpMatrix a, SpMatrix a) -- ^ Q, H+arnoldi aa b kn | n == nb = return (fromCols qvfin, fromListSM (nmax + 1, nmax) hhfin)+ | otherwise = throwM (MatVecSizeMismatchException "arnoldi" (m,n) nb)+ where+ (qvfin, hhfin, nmax, _) = execState (modifyUntil tf arnoldiStep) arnInit + tf (_, _, ii, fbreak) = ii == kn || fbreak -- termination criterion+ (m, n) = (nrows aa, ncols aa)+ nb = dim b+ arnInit = (qv1, hh1, 1, False) where+ q0 = normalize2 b -- starting basis vector+ aq0 = aa #> q0 -- A q0+ h11 = q0 `dot` aq0 + q1nn = aq0 ^-^ (h11 .* q0)+ hh1 = V.fromList [(0, 0, h11), (1, 0, h21)] where + h21 = norm2' q1nn+ q1 = normalize2 q1nn -- q1 `dot` q0 ~ 0+ qv1 = V.fromList [q0, q1]+ arnoldiStep (qv, hh, i, _) = (qv', hh', i + 1, breakf) where+ qi = V.last qv+ aqi = aa #> qi+ hhcoli = fmap (`dot` aqi) qv -- H_{1, i}, H_{2, i}, .. , H_{m + 1, i}+ zv = zeroSV m+ qipnn =+ aqi ^-^ V.foldl' (^+^) zv (V.zipWith (.*) hhcoli qv) -- unnormalized q_{i+1}+ qipnorm = norm2' qipnn -- normalization factor H_{i+1, i}+ qip = normalize2 qipnn -- q_{i + 1}+ hh' = (V.++) hh (indexed2 $ V.snoc hhcoli qipnorm) where -- update H+ indexed2 v = V.zip3 ii jj v+ ii = V.fromList [0 .. n] -- nth col of upper Hessenberg has `n+1` nz+ jj = V.replicate (n + 1) i -- `n+1` replicas of `i`+ qv' = V.snoc qv qip -- append q_{i+1} to Krylov basis Q_i+ breakf | nearZero qipnorm = True -- breakdown condition+ | otherwise = False++++++++-- * Preconditioning++-- | Partition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts+diagPartitions :: SpMatrix a+ -> (SpMatrix a, SpMatrix a, SpMatrix a) -- ^ Subdiagonal, diagonal, superdiagonal partitions+diagPartitions aa = (e,d,f) where+ e = extractSubDiag aa+ d = extractDiag aa+ f = extractSuperDiag aa+++-- ** Jacobi preconditioner++-- | The Jacobi preconditioner is just the reciprocal of the diagonal +jacobiPre :: Fractional a => SpMatrix a -> SpMatrix a+jacobiPre x = recip <$> extractDiag x+++-- ** Incomplete LU++-- | Used for Incomplete LU : remove entries in `m` corresponding to zero entries in `m2` (this is called ILU(0) in the preconditioner literature)+ilu0 :: (Scalar (SpVector t) ~ t, Elt t, VectorSpace (SpVector t),+ Epsilon t, MonadThrow m) =>+ SpMatrix t+ -> m (SpMatrix t, SpMatrix t) -- ^ L, U (with holes)+ilu0 aa = do+ (l, u) <- lu aa+ let lh = sparsifyLU l aa+ uh = sparsifyLU u aa+ sparsifyLU m m2 = ifilterSM f m where+ f i j _ = isJust (lookupSM m2 i j)+ return (lh, uh)+++-- ** SSOR++-- | Symmetric Successive Over-Relaxation. `mSsor aa omega` : if `omega = 1` it returns the symmetric Gauss-Seidel preconditioner. When ω = 1, the SOR reduces to Gauss-Seidel; when ω > 1 and ω < 1, it corresponds to over-relaxation and under-relaxation, respectively.+mSsor :: (MatrixRing (SpMatrix b), Fractional b) =>+ SpMatrix b+ -> b -- ^ relaxation factor+ -> (SpMatrix b, SpMatrix b) -- ^ Left, right factors+mSsor aa omega = (l, r) where+ (e, d, f) = diagPartitions aa+ n = nrows e+ l = (eye n ^-^ scale omega e) ## reciprocal d+ r = d ^-^ scale omega f +++++++++-- * Linear solver, LU-based++-- | Direct solver based on a triangular factorization of the system matrix.+luSolve :: (Scalar (SpVector t) ~ t, MonadThrow m, Elt t, InnerSpace (SpVector t), Epsilon t) =>+ SpMatrix t -- ^ Lower triangular+ -> SpMatrix t -- ^ Upper triangular+ -> SpVector t -- ^ r.h.s.+ -> m (SpVector t)+luSolve ll uu b+ | isLowerTriSM ll && isUpperTriSM uu = do+ w <- triLowerSolve ll b+ triUpperSolve uu w+ | otherwise = throwM (NonTriangularException "luSolve")++-- | Forward substitution solver+triLowerSolve :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),+ Epsilon t, MonadThrow m) =>+ SpMatrix t -- ^ Lower triangular+ -> SpVector t+ -> m (SpVector t)+triLowerSolve ll b = do+ let q (_, i) = i == nb+ nb = svDim b+ oops i = throwM (NeedsPivoting "triLowerSolve" (unwords ["L", show (i, i)]) :: MatrixException Double)+ lStep (ww, i) = do -- (ww', i + 1) where+ let+ lii = ll @@ (i, i)+ bi = b @@ i+ wi = (bi - r)/lii where+ r = extractSubRow ll i (0, i-1) `dot` takeSV i ww+ if isNz lii+ then return (insertSpVector i wi ww, i + 1)+ else oops i+ lInit = do -- (ww0, 1) where+ let+ l00 = ll @@ (0, 0)+ b0 = b @@ 0+ w0 = b0 / l00+ if isNz l00+ then return (insertSpVector 0 w0 $ zeroSV (dim b), 1)+ else oops (0 :: Int)+ l0 <- lInit + (v, _) <- MTS.execStateT (modifyUntilM q lStep) l0+ return $ sparsifySV v++++-- NB in the computation of `xi` we must rebalance the subrow indices (extractSubRow_RK) because `dropSV` does that too, in order to take the inner product with consistent index pairs+-- | Backward substitution solver+triUpperSolve :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),+ Epsilon t, MonadThrow m) =>+ SpMatrix t -- ^ Upper triangular+ -> SpVector t+ -> m (SpVector t)+triUpperSolve uu w = do + let q (_, i) = i == (- 1)+ nw = svDim w+ oops i = throwM (NeedsPivoting "triUpperSolve" (unwords ["U", show (i, i)]) :: MatrixException Double)+ uStep (xx, i) = do+ let uii = uu @@ (i, i) + wi = w @@ i+ r = extractSubRow_RK uu i (i + 1, nw - 1) `dot` dropSV (i + 1) xx + xi = (wi - r) / uii+ if isNz uii+ then return (insertSpVector i xi xx, i - 1)+ else oops i + uInit = do + let i = nw - 1+ u00 = uu @@! (i, i)+ w0 = w @@ i+ x0 = w0 / u00+ if isNz u00+ then return (insertSpVector i x0 (zeroSV nw), i - 1)+ else oops (0 :: Int)+ u0 <- uInit + (x, _) <- MTS.execStateT (modifyUntilM q uStep) u0+ return $ sparsifySV x ++++++ ++++++-- * Iterative linear solvers+++-- ** GMRES++-- | Given a linear system `A x = b` where `A` is an (m x m) real-valued matrix, the GMRES method finds an approximate solution `xhat` such that the Euclidean norm of the residual `A xhat - b` is minimized. `xhat` is spanned by the order-`n` Krylov subspace of (A, b).+-- In this implementation:+-- 1) the Arnoldi iteration is carried out until numerical breakdown (therefore yielding _at_most_ `m+1` Krylov basis vectors)+-- 2) the resulting Hessenberg matrix H is factorized in QR form (H = Q R)+-- 3) the Krylov-subspace solution `yhat` is found by backsubstitution (since R is upper-triangular)+-- 4) the approximate solution in the original space `xhat` is computed using the Krylov basis, `xhat = Q_n yhat`+--+-- Many optimizations are possible, for example interleaving the QR factorization (and the subsequent triangular solve) with the Arnoldi process (and employing an updating QR factorization which only requires one Givens' rotation at every update). ++gmres :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,+ Elt t, Normed (SpVector t), LinearVectorSpace (SpVector t), Epsilon t,+ MonadThrow m) =>+ SpMatrix t -> SpVector t -> m (SpVector t)+gmres aa b = do+ let m = ncols aa+ (qa, ha) <- arnoldi aa b m -- at most m steps of Arnoldi (aa, b)+ -- b' = (transposeSe qa) #> b+ let b' = norm2' b .* ei mp1 1 -- b rotated back to canonical basis by Q^T+ where mp1 = nrows ha -- = 1 + (# Arnoldi iterations)+ (qh, rh) <- qr ha+ let rhs' = takeSV (dim b' - 1) (transpose qh #> b')+ rh' = takeRows (nrows rh - 1) rh -- last row of `rh` is 0+ yhat <- triUpperSolve rh' rhs'+ let qa' = takeCols (ncols qa - 1) qa -- we don't use last column of Krylov basis + return $ qa' #> yhat++++ +-- ** CGNE++data CGNE a =+ CGNE {_xCgne , _rCgne, _pCgne :: SpVector a} deriving Eq+instance Show a => Show (CGNE a) where+ show (CGNE x r p) = "x = " ++ show x ++ "\n" +++ "r = " ++ show r ++ "\n" +++ "p = " ++ show p ++ "\n"++cgneInit :: (MatrixType (SpVector a) ~ SpMatrix a,+ LinearVectorSpace (SpVector a)) =>+ SpMatrix a -> SpVector a -> SpVector a -> CGNE a+cgneInit aa b x0 = CGNE x0 r0 p0 where+ r0 = b ^-^ (aa #> x0) -- residual of initial guess solution+ p0 = transposeSM aa #> r0++cgneStep :: (MatrixType (SpVector a) ~ SpMatrix a,+ LinearVectorSpace (SpVector a), InnerSpace (SpVector a),+ Fractional (Scalar (SpVector a))) =>+ SpMatrix a -> CGNE a -> CGNE a+cgneStep aa (CGNE x r p) = CGNE x1 r1 p1 where+ alphai = (r `dot` r) / (p `dot` p)+ x1 = x ^+^ (alphai .* p)+ r1 = r ^-^ (alphai .* (aa #> p))+ beta = (r1 `dot` r1) / (r `dot` r)+ p1 = transpose aa #> r ^+^ (beta .* p)++++++-- ** BCG++data BCG a =+ BCG { _xBcg, _rBcg, _rHatBcg, _pBcg, _pHatBcg :: SpVector a } deriving Eq++bcgInit :: LinearVectorSpace (SpVector a) =>+ MatrixType (SpVector a) -> SpVector a -> SpVector a -> BCG a+bcgInit aa b x0 = BCG x0 r0 r0hat p0 p0hat where+ r0 = b ^-^ (aa #> x0) + r0hat = r0+ p0 = r0+ p0hat = r0++bcgStep :: (MatrixType (SpVector a) ~ SpMatrix a,+ LinearVectorSpace (SpVector a), InnerSpace (SpVector a),+ Fractional (Scalar (SpVector a))) =>+ SpMatrix a -> BCG a -> BCG a +bcgStep aa (BCG x r rhat p phat) = BCG x1 r1 rhat1 p1 phat1 where+ aap = aa #> p+ alpha = (r `dot` rhat) / (aap `dot` phat)+ x1 = x ^+^ (alpha .* p)+ r1 = r ^-^ (alpha .* aap)+ rhat1 = rhat ^-^ (alpha .* (transpose aa #> phat))+ beta = (r1 `dot` rhat1) / (r `dot` rhat)+ p1 = r1 ^+^ (beta .* p)+ phat1 = rhat1 ^+^ (beta .* phat)++instance Show a => Show (BCG a) where+ show (BCG x r rhat p phat) = "x = " ++ show x ++ "\n" +++ "r = " ++ show r ++ "\n" +++ "r_hat = " ++ show rhat ++ "\n" +++ "p = " ++ show p ++ "\n" +++ "p_hat = " ++ show phat ++ "\n"+++-- ** CGS++data CGS a = CGS { _x, _r, _p, _u :: SpVector a} deriving Eq++cgsInit :: LinearVectorSpace (SpVector a) =>+ MatrixType (SpVector a) -> SpVector a -> SpVector a -> CGS a+cgsInit aa b x0 = CGS x0 r0 r0 r0 where+ r0 = b ^-^ (aa #> x0) -- residual of initial guess solution++cgsStep :: (V (SpVector a), Fractional (Scalar (SpVector a))) =>+ MatrixType (SpVector a) -> SpVector a -> CGS a -> CGS a+cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1+ where+ aap = aa #> p+ alphaj = (r `dot` rhat) / (aap `dot` rhat)+ q = u ^-^ (alphaj .* aap)+ xj1 = x ^+^ (alphaj .* (u ^+^ q)) -- updated solution+ rj1 = r ^-^ (alphaj .* (aa #> (u ^+^ q))) -- updated residual+ betaj = (rj1 `dot` rhat) / (r `dot` rhat)+ uj1 = rj1 ^+^ (betaj .* q)+ pj1 = uj1 ^+^ (betaj .* (q ^+^ (betaj .* p)))+++instance (Show a) => Show (CGS a) where+ show (CGS x r p u) = "x = " ++ show x ++ "\n" +++ "r = " ++ show r ++ "\n" +++ "p = " ++ show p ++ "\n" +++ "u = " ++ show u ++ "\n"++++++-- ** BiCGSTAB++data BICGSTAB a =+ BICGSTAB { _xBicgstab, _rBicgstab, _pBicgstab :: SpVector a} deriving Eq++bicgsInit :: LinearVectorSpace (SpVector a) =>+ MatrixType (SpVector a) -> SpVector a -> SpVector a -> BICGSTAB a+bicgsInit aa b x0 = BICGSTAB x0 r0 r0 where+ r0 = b ^-^ (aa #> x0) -- residual of initial guess solution++bicgstabStep :: (V (SpVector a), Fractional (Scalar (SpVector a))) =>+ MatrixType (SpVector a) -> SpVector a -> BICGSTAB a -> BICGSTAB a+bicgstabStep aa r0hat (BICGSTAB x r p) = BICGSTAB xj1 rj1 pj1 where+ aap = aa #> p+ alphaj = (r <.> r0hat) / (aap <.> r0hat)+ sj = r ^-^ (alphaj .* aap)+ aasj = aa #> sj+ omegaj = (aasj <.> sj) / (aasj <.> aasj)+ xj1 = x ^+^ (alphaj .* p) ^+^ (omegaj .* sj) -- updated solution+ rj1 = sj ^-^ (omegaj .* aasj)+ betaj = (rj1 <.> r0hat)/(r <.> r0hat) * alphaj / omegaj+ pj1 = rj1 ^+^ (betaj .* (p ^-^ (omegaj .* aap)))++instance Show a => Show (BICGSTAB a) where+ show (BICGSTAB x r p) = "x = " ++ show x ++ "\n" +++ "r = " ++ show r ++ "\n" +++ "p = " ++ show p ++ "\n"+++++++-- * Moore-Penrose pseudoinverse+-- | Least-squares approximation of a rectangular system of equaitons. Uses <\\> for the linear solve+pinv :: (MatrixType v ~ SpMatrix a, LinearSystem v, Epsilon a,+ MonadThrow m, MonadIO m) =>+ SpMatrix a -> v -> m v+pinv aa b = aa #~^# aa <\> atb where+ atb = transpose aa #> b+++++-- | Interface method to individual linear solvers+linSolve0 method aa b x0+ | m /= nb = throwM (MatVecSizeMismatchException "linSolve0" dm nb)+ | otherwise = solve aa b where+ solve aa' b' | isDiagonalSM aa' = return $ reciprocal aa' #> b' -- diagonal solve+ | otherwise = xHat+ xHat = case method of+ BICGSTAB_ -> solver "BICGSTAB" nits _xBicgstab (bicgstabStep aa r0hat) (bicgsInit aa b x0)+ BCG_ -> solver "BCG" nits _xBcg (bcgStep aa) (bcgInit aa b x0)+ CGS_ -> solver "CGS" nits _x (cgsStep aa r0hat) (cgsInit aa b x0)+ GMRES_ -> gmres aa b+ CGNE_ -> solver "CGNE" nits _xCgne (cgneStep aa) (cgneInit aa b x0)+ r0hat = b ^-^ (aa #> x0)+ nits = 200+ dm@(m,n) = dim aa+ nb = dim b+ solver fname nitermax fproj stepf initf = do+ xf <- untilConvergedG fname config (const True) stepf initf+ return $ fproj xf+ where+ config = IterConf nitermax True fproj prd+ ++-- * Linear solver interface++-- -- TFQMR is in LinearSolvers.Experimental for now+data LinSolveMethod = GMRES_ | CGNE_ | BCG_ | CGS_ | BICGSTAB_ deriving (Eq, Show) +++++-- -- -- | linSolve using the GMRES method as default+instance LinearSystem (SpVector Double) where+ aa <\> b = linSolve0 GMRES_ aa b (mkSpVR n $ replicate n 0.1)+ where n = ncols aa++instance LinearSystem (SpVector (Complex Double)) where+ aa <\> b = linSolve0 GMRES_ aa b (mkSpVC n $ replicate n 0.1)+ where n = ncols aa++++++++++++-- | TODO : if system is poorly conditioned, is it better to warn the user or just switch solvers (e.g. via the pseudoinverse) ?++-- linSolveQR aa b init f1 stepf+-- | isInfinite k = do+-- tell "linSolveQR : rank-deficient system"+-- | otherwise = do+-- solv aa b init+-- where+-- (q, r) = qr aa+-- k = conditionNumberSM r+-- solv aa b init = execState (untilConverged f1 stepf) init++++ ++++++-- test data++-- aa4 : eigenvalues 1 (mult.=2) and -1+aa4 :: SpMatrix Double+aa4 = fromListDenseSM 3 [3,2,-2,2,2,-1,6,5,-4] ++aa4c :: SpMatrix (Complex Double)+aa4c = toC <$> aa4+++++++++-- -- * Random arrays++-- randArray :: PrimMonad m => Int -> Double -> Double -> m [Double]+-- randArray n mu sig = do+-- g <- MWC.create+-- replicateM n (MWC.normal mu sig g)+ ++++-- -- * Random matrices and vectors++-- -- |Dense SpMatrix+-- randMat :: PrimMonad m => Int -> m (SpMatrix Double)+-- randMat n = do+-- g <- MWC.create+-- aav <- replicateM (n^2) (MWC.normal 0 1 g)+-- let ii_ = [0 .. n-1]+-- (ix_,iy_) = unzip $ concatMap (zip ii_ . replicate n) ii_+-- return $ fromListSM (n,n) $ zip3 ix_ iy_ aav++-- -- | Dense SpVector +-- randVec :: PrimMonad m => Int -> m (SpVector Double)+-- randVec n = do+-- g <- MWC.create+-- bv <- replicateM n (MWC.normal 0 1 g)+-- let ii_ = [0..n-1]+-- return $ fromListSV n $ zip ii_ bv++++-- -- | Sparse SpMatrix+-- randSpMat :: Int -> Int -> IO (SpMatrix Double)+-- randSpMat n nsp | nsp > n^2 = error "randSpMat : nsp must be < n^2 "+-- | otherwise = do+-- g <- MWC.create+-- aav <- replicateM nsp (MWC.normal 0 1 g)+-- ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)+-- jj <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)+-- return $ fromListSM (n,n) $ zip3 ii jj aav++-- -- | Sparse SpVector+-- randSpVec :: Int -> Int -> IO (SpVector Double)+-- randSpVec n nsp | nsp > n = error "randSpVec : nsp must be < n"+-- | otherwise = do+-- g <- MWC.create+-- aav <- replicateM nsp (MWC.normal 0 1 g)+-- ii <- replicateM nsp (MWC.uniformR (0, n-1) g :: IO Int)+-- return $ fromListSV n $ zip ii aav
test/LibSpec.hs view
@@ -1,4 +1,5 @@-{-# language ScopedTypeVariables #-}+{-# LANGUAGE FlexibleContexts, TypeFamilies #-}+{-# language ScopedTypeVariables, FlexibleInstances #-} ----------------------------------------------------------------------------- -- | -- Copyright : (C) 2016 Marco Zocca@@ -10,25 +11,32 @@ ----------------------------------------------------------------------------- module LibSpec where +import Control.Exception.Common+import Data.Sparse.Common import Numeric.LinearAlgebra.Sparse -- import Numeric.LinearAlgebra.Class -import Control.Applicative (liftA2)--- import Control.Monad (liftM, liftM2, replicateM)-import Control.Monad.Primitive-import Data.Foldable (foldrM)+-- import Control.Applicative (liftA2)+-- import Control.Monad (replicateM)+import Control.Monad.Catch -- (MonadThrow (..))+import Control.Monad.IO.Class+import Control.Monad.State+-- import Data.Foldable (foldrM) -import Data.Sparse.Common+import Data.Complex+-- import Data.Either (either)+import Data.Typeable -import Control.Monad.State.Strict (execState)+import Data.VectorSpace hiding (magnitude) -import qualified System.Random.MWC as MWC-import qualified System.Random.MWC.Distributions as MWC+-- 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--+import Test.Hspec.QuickCheck+import Test.QuickCheck main :: IO ()@@ -38,19 +46,26 @@ 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+ describe "Numeric.LinearAlgebra.Sparse : Library" $ do+ prop "Subtraction is cancellative" $ \(x :: SpVector Double) ->+ norm2Sq (x ^-^ x) `shouldBe` zeroV+ it "<.> : inner product (Real)" $+ tv0 <.> tv0 `shouldBe` 61+ it "<.> : inner product (Complex)" $+ tvc2 <.> tvc3 `shouldBe` 2 :+ (-2) + it "transpose : sparse matrix transpose" $+ transpose m1 `shouldBe` m1t+ it "(#>) : matrix-vector product (Real)" $+ nearZero ( norm2Sq ((aa0 #> x0true) ^-^ b0 )) `shouldBe` True+ it "(<#) : vector-matrix product (Real)" $+ nearZero ( norm2Sq ((x0true <# aa0) ^-^ aa0tx0 ))`shouldBe` True + it "(##) : matrix-matrix product (Real, square)" $ + (m1 ## m2) `shouldBe` m1m2+ it "(##) : matrix-matrix product (Real, rectangular)" $ do+ (m1' ## m2') `shouldBe` m1m2'+ (m2' ## m1') `shouldBe` m2m1'+ it "(##) : matrix-matrix product (Complex)" $+ (aa3c ## aa3c) `shouldBe` aa3cx it "eye : identity matrix" $ infoSM (eye 10) `shouldBe` SMInfo 10 0.1 it "insertCol : insert a column in a SpMatrix" $@@ -65,109 +80,553 @@ 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+ 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 "untilConvergedG0 : early termination by iteration count and termination by convergence" $ + let+ n1 = 4+ nexp1 = fromIntegral n1 / fromIntegral (2^n1) -- 0.25+ f x = x/2+ mm1 = untilConvergedG0 "blah"+ (IterConf n1 False id print) (1/(2^n1)) f (fromIntegral n1 :: Double)+ n2 = 2^16+ mm2 = untilConvergedG0 "blah"+ (IterConf n2 False id print) (1/(2^n2)) f (fromIntegral n1 :: Double)+ eh (NotConvergedE _ _ x) = return x+ in+ do x1 <- mm1 `catch` eh+ x1 `shouldBe` nexp1+ x2 <- mm2 `catch` eh+ nearZero x2 `shouldBe` True+ + ++ describe "QuickCheck properties:" $ do+ -- prop "prop_matSPD_vec : (m #^# m) is symmetric positive definite" $+ -- \(PropMatSPDVec (m :: SpMatrix Double) v) -> prop_spd m v+ prop "prop_dot : (v <.> v) ~= 1 if ||v|| == 1" $+ \(v :: SpVector Double) -> prop_dot v+ prop "prop_matMat1 : (A ## B)^T == (B^T ## A^T)" $+ \p@(PropMatMat (_ :: SpMatrix Double) _) -> prop_matMat1 p+ prop "prop_matMat2 : M^T ##^ M == M #^# M^T" $+ \p@(PropMat (_ :: SpMatrix Double)) -> prop_matMat2 p+ -- prop "prop_QR : Q R = A, Q is orthogonal, R is upper triangular" $+ -- \p@(PropMatI (_ :: SpMatrix Double)) -> prop_QR p+ -- prop "prop_Cholesky" $ \p@(PropMat_SPD (_ :: SpMatrix Double)) -> prop_Cholesky p+ -- prop "prop_linSolve GMRES" $ prop_linSolve GMRES_+ -- prop "aa2 is positive semidefinite" $ \(v :: SpVector Double) ->+ -- prop_psd aa2 v+ + describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers (Real)" $ do -- it "TFQMR (2 x 2 dense)" $- -- normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True it "GMRES (2 x 2 dense)" $- nearZero (normSq (linSolve GMRES_ aa0 b0 ^-^ x0true)) `shouldBe` True- it "GMRES (3 x 3 sparse, s.p.d.)" $- nearZero (normSq (linSolve GMRES_ aa2 b2 ^-^ x2)) `shouldBe` True+ checkLinSolveR GMRES_ aa0 b0 x0true >>= (`shouldBe` True)+ it "GMRES (3 x 3 sparse, symmetric pos.def.)" $+ checkLinSolveR GMRES_ aa2 b2 x2 >>= (`shouldBe` True) it "GMRES (4 x 4 sparse)" $- nearZero (normSq (linSolve GMRES_ aa1 b1 ^-^ x1)) `shouldBe` True + checkLinSolveR GMRES_ aa1 b1 x1 >>= (`shouldBe` True)+ it "GMRES (5 x 5 sparse)" $+ checkLinSolveR GMRES_ tm7 tvb7 tvx7 >>= (`shouldBe` True) it "BCG (2 x 2 dense)" $- nearZero (normSq (linSolve BCG_ aa0 b0 ^-^ x0true)) `shouldBe` True- it "BCG (3 x 3 sparse, s.p.d.)" $- nearZero (normSq (linSolve BCG_ aa2 b2 ^-^ x2)) `shouldBe` True + checkLinSolveR BCG_ aa0 b0 x0true >>= (`shouldBe` True)+ it "BCG (3 x 3 sparse, symmetric pos.def.)" $+ checkLinSolveR BCG_ aa2 b2 x2 >>= (`shouldBe` True) -- it "BiCGSTAB (2 x 2 dense)" $ -- nearZero (normSq (linSolve BICGSTAB_ aa0 b0 ^-^ x0true)) `shouldBe` True- it "BiCGSTAB (3 x 3 sparse, s.p.d.)" $ - nearZero (normSq (linSolve BICGSTAB_ aa2 b2 ^-^ x2)) `shouldBe` True + it "BiCGSTAB (3 x 3 sparse, symmetric pos.def.)" $ + checkLinSolveR BICGSTAB_ aa2 b2 x2 >>= (`shouldBe` True) it "CGS (2 x 2 dense)" $ - nearZero (normSq (linSolve CGS_ aa0 b0 ^-^ x0true)) `shouldBe` True- it "CGS (3 x 3 sparse, s.p.d.)" $ - nearZero (normSq (linSolve CGS_ aa2 b2 ^-^ x2)) `shouldBe` True - describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ + checkLinSolveR CGS_ aa0 b0 x0true >>= (`shouldBe` True)+ it "CGS (3 x 3 sparse, SPD)" $ + checkLinSolveR CGS_ aa2 b2 x2 >>= (`shouldBe` True)+ + describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers (Real)" $ it "luSolve (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+ checkLuSolve aa1 b1 >>= (`shouldBe` True)+ describe "Numeric.LinearAlgebra.Sparse : QR factorization (Real)" $ do it "qr (3 x 3 dense)" $ - checkQr tm2 `shouldBe` True- describe "Numeric.LinearAlgebra.Sparse : LU decomposition" $ do+ checkQr tm2 >>= (`shouldBe` True)+ it "qr (4 x 4 sparse)" $+ checkQr tm4 >>= (`shouldBe` True)+ it "qr (5 x 5 sparse)" $+ checkQr tm7 >>= (`shouldBe` True)+ describe "Numeric.LinearAlgebra.Sparse : QR factorization (Complex)" $ do+ it "qr (2 x 2 dense)" $+ checkQr aa3cx >>= (`shouldBe` True)+ it "qr (3 x 3 dense)" $+ checkQr aa4c >>= (`shouldBe` True) + describe "Numeric.LinearAlgebra.Sparse : LU factorization (Real)" $ do+ it "lu (3 x 3 dense)" $+ checkLu tm2 >>= (`shouldBe` True) 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)" $ + checkLu tm6 >>= (`shouldBe` True)+ it "lu (5 x 5 sparse)" $+ checkLu tm7 >>= (`shouldBe` True)+ describe "Numeric.LinearAlgebra.Sparse : Cholesky factorization (Real, symmetric pos.def.)" $ it "chol (5 x 5 sparse)" $- checkChol tm7 `shouldBe` True- describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration, early breakdown detection" $ do + checkChol tm7 >>= (`shouldBe` True)+ describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration (Real)" $ do it "arnoldi (4 x 4 dense)" $- checkArnoldi tm6 4 `shouldBe` True+ checkArnoldi tm6 4 >>= (`shouldBe` True) it "arnoldi (5 x 5 sparse)" $- checkArnoldi tm7 5 `shouldBe` True + checkArnoldi tm7 5 >>= (`shouldBe` True) +{- linear systems -}++-- checkLinSolve method aa b x x0r =+-- either+-- (error . show)+-- (\xhat -> nearZero (norm2Sq (x ^-^ xhat)))+-- (linSolve0 method aa b x0r)++checkLinSolve' method aa b x x0r =+ nearZero . norm2 <$> linSolve0 method aa b x0r -- `catch` eh+ -- where+ -- eh (NotConvergedE _ i xhat) = return $ xhat ^-^ x++checkLinSolve method aa b x x0r = do+ xhat <- linSolve0 method aa b x0r+ return $ nearZero $ norm2 (x ^-^ xhat)++checkLinSolveR+ :: (MonadIO m, MonadCatch m) =>+ LinSolveMethod + -> SpMatrix Double -- ^ operator+ -> SpVector Double -- ^ r.h.s+ -> SpVector Double -- ^ candidate solution+ -> m Bool+checkLinSolveR method aa b x = checkLinSolve method aa b x x0r where+ x0r = mkSpVR n $ replicate n 0.1+ n = ncols aa++checkLinSolveC+ :: (MonadIO m, MonadCatch m) =>+ LinSolveMethod+ -> SpMatrix (Complex Double)+ -> SpVector (Complex Double)+ -> SpVector (Complex Double)+ -> m Bool+checkLinSolveC method aa b x = checkLinSolve method aa b x x0r where+ x0r = mkSpVC n $ replicate n (0.1 :+ 0.1)+ n = ncols aa+++++++-- {- Givens rotation-}+checkGivens1 :: (MonadThrow m, Elt a, MatrixRing (SpMatrix a), Epsilon a) =>+ SpMatrix a -> IxRow -> IxCol -> m (a, Bool)+checkGivens1 tm i j = do -- (rij, nearZero rij) where+ g <- givens tm i j+ let r = g ## tm+ rij = r @@ (i, j)+ return (rij, nearZero rij)++ {- QR-} +checkQr :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+ SpMatrix a -> m Bool+checkQr a = do+ (q, r) <- qr a+ let c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)+ c2 = isOrthogonalSM q+ c3 = isUpperTriSM r+ return $ c1 && c2 && c3 -checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkQr a = c1 && c2 && c3 where- (q, r) = qr a- c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)- c2 = isOrthogonalSM q- c3 = isUpperTriSM r +stepQR a = do+ (q, r) <- qr a+ return $ r #~# q +stepQRix (i, a) = do+ a' <- stepQR a+ return (i + 1, a')+ + {- LU -} -checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkLu a = c1 && c2 where- (l, u) = lu a- c1 = nearZero (normFrobenius ((l #~# u) ^-^ a))- c2 = isUpperTriSM u && isLowerTriSM l+checkLu :: (Scalar (SpVector t) ~ t, Elt t, MatrixRing (SpMatrix t),+ VectorSpace (SpVector t), Epsilon t, MonadThrow m) =>+ SpMatrix t -> m Bool+checkLu a = do + (l, u) <- lu a+ let c1 = nearZero (normFrobenius ((l #~# u) ^-^ a))+ c2 = isUpperTriSM u && isLowerTriSM l+ return (c1 && c2) {- Cholesky -} -checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkChol a = c1 && c2 where- l = chol a- c1 = nearZero $ normFrobenius ((l ##^ l) ^-^ a)- c2 = isLowerTriSM l+checkChol :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+ SpMatrix a -> m Bool+checkChol a = do -- c1 && c2 where+ l <- chol a+ let c1 = nearZero $ normFrobenius ((l ##^ l) ^-^ a)+ c2 = isLowerTriSM l+ return $ c1 && c2 {- 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- +checkLuSolve :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,+ Elt t, Normed (SpVector t), LinearVectorSpace (SpVector t),+ Epsilon t, MonadThrow m) =>+ SpMatrix t -> SpVector t -> m Bool+checkLuSolve amat rhs = do+ (lmat, umat) <- lu amat+ xlu <- luSolve lmat umat rhs+ return $ nearZero (norm2Sq ( (lmat #> (umat #> xlu)) ^-^ rhs ))+ {- Arnoldi iteration -}-checkArnoldi :: (Epsilon a, Floating a, Eq a) => SpMatrix a -> Int -> Bool-checkArnoldi aa kn = nearZero (normFrobenius $ lhs ^-^ rhs) where- b = onesSV (nrows aa)- (q, h) = arnoldi aa b kn- (m, n) = dim q- q' = extractSubmatrix q (0, m - 1) (0, n - 2) -- q' = all but one column of q- rhs = q #~# h- lhs = aa #~# q' +checkArnoldi :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,+ Normed (SpVector t), MatrixRing (SpMatrix t),+ LinearVectorSpace (SpVector t), Epsilon t, MonadThrow m) =>+ SpMatrix t -> Int -> m Bool+checkArnoldi aa kn = do -- nearZero (normFrobenius $ lhs ^-^ rhs) where+ let b = onesSV (nrows aa)+ (q, h) <- arnoldi aa b kn+ let (m, n) = dim q+ q' = extractSubmatrix q (0, m - 1) (0, n - 2) -- q' = all but one column of q+ rhs = q #~# h+ lhs = aa #~# q'+ return $ nearZero (normFrobenius $ lhs ^-^ rhs) ++++++++++++-- * Arbitrary newtypes and instances for QuickCheck++-- | helpers+sized2 :: (Int -> Int -> Gen a) -> Gen a+sized2 f = sized $ \i -> sized $ \j -> f i j++sized3 :: (Int -> Int -> Int -> Gen a) -> Gen a+sized3 f = sized $ \i -> sized $ \j -> sized $ \k -> f i j k+++whenFail1 :: Testable prop => (t -> IO ()) -> (t -> prop) -> t -> Property+whenFail1 io p x = whenFail (io x) (property $ p x)+++sampleWith :: (a -> IO b) -> Gen a -> IO ()+sampleWith pf g = do+ cases <- sample' g+ mapM_ pf cases++sampleSp :: PrintDense a => Gen a -> IO ()+sampleSp = sampleWith prd++++-- | (m * n) random sparse matrix having d elements+genSpM0 :: Arbitrary a => Int -> Int -> Int -> Gen (SpMatrix a)+genSpM0 m n d = do+ -- let d = floor (sqrt $ fromIntegral (m * n)) :: Int+ i_ <- vectorOf d $ choose (0, m-1)+ j_ <- vectorOf d $ choose (0, n-1) + x_ <- vector d+ return $ fromListSM (m,n) $ zip3 i_ j_ x_++-- | Random (m * n) sparse matrix having sqrt(m * n) elements+genSpM :: Arbitrary a => Int -> Int -> Gen (SpMatrix a) +genSpM m n = genSpM0 m n $ floor (sqrt $ fromIntegral (m * n))+++++-- | (m * n) random DENSE matrix+genSpMDense :: (Arbitrary a, Num a) => Int -> Int -> Gen (SpMatrix a)+genSpMDense m n = do+ xs <- vector (m*n)+ let ii = concatMap (replicate n) [0..m-1]+ jj = concat $ replicate m [0..n-1]+ return $ fromListSM (m,n) $ zip3 ii jj xs++++-- | Order n diagonal SpMatrix with constant elements+genSpMConstDiagonal ::+ (Arbitrary a, Ord a, Num a) => (a -> Bool) -> Int -> Gen (SpMatrix a)+genSpMConstDiagonal f n = do+ x <- arbitrary `suchThat` f+ return $ mkDiagonal n (replicate n x)+++-- | Order-n diagonal SpMatrix+genSpMDiagonal :: Arbitrary a => ([a] -> Bool) -> Int -> Gen (SpMatrix a)+genSpMDiagonal f n = do+ xs <- vector n `suchThat` f+ return $ mkDiagonal n xs++++++++++++ +++++++++++++++++++++-- | Random sparse vector+genSpV0 :: Arbitrary a => Int -> Int -> Gen (SpVector a)+genSpV0 d n = do+ i_ <- vectorOf d $ choose (0, n -1)+ v_ <- vector d+ return $ fromListSV n (zip i_ v_)++genSpV :: Arbitrary a => Int -> Gen (SpVector a)+genSpV n = genSpV0 (floor (sqrt $ fromIntegral n) :: Int) n+++-- | Random dense vector+genSpVDense :: (Epsilon a, Arbitrary a) => Int -> Gen (SpVector a)+genSpVDense n = do+ v <- vector n `suchThat` any isNz+ return $ fromListDenseSV n v++++++-- | An Arbitrary SpVector such that at least one entry is nonzero+instance Arbitrary (SpVector Double) where+ arbitrary = sized genSpV `suchThat` any isNz +++-- | An arbitrary square SpMatrix+newtype PropMat0 a = PropMat0 (SpMatrix a) deriving (Eq, Show)+instance Arbitrary (PropMat0 Double) where+ arbitrary = sized (\n -> PropMat0 <$> genSpM n n) ++ +-- | An arbitrary SpMatrix+newtype PropMat a = PropMat { unPropMat :: SpMatrix a} deriving (Eq, Show)+instance Arbitrary (PropMat Double) where+ arbitrary = sized2 (\m n -> PropMat <$> genSpM m n) `suchThat` ((> 2) . nrows . unPropMat)++-- nzDim :: SpMatrix a -> Bool+-- nzDim mm = let (m, n) = dim mm in m > 2 && n > 2++-- sizedCon :: (a -> Bool) -> (Int -> Gen a) -> Gen a+-- sizedCon f genf = sized genf `suchThat` f+++++++-- | An arbitrary DENSE SpMatrix+newtype PropMatDense a = PropMatDense {unPropMatDense :: SpMatrix a} deriving (Eq, Show)+instance Arbitrary (PropMatDense Double) where+ arbitrary = sized2 (\m n -> PropMatDense <$> genSpM m n) `suchThat` ((> 2) . nrows . unPropMatDense)++++-- | An arbitrary SpMatrix with identity diagonal +newtype PropMatI a = PropMatI {unPropMatI :: SpMatrix a} deriving (Eq)+instance Show a => Show (PropMatI a) where show = show . unPropMatI+instance Arbitrary (PropMatI Double) where+ arbitrary = sized (\m -> PropMatI <$> genSpMI m) `suchThat` ((> 2) . nrows . unPropMatI)++genSpMI :: (Num a, Arbitrary a) => Int -> Gen (SpMatrix a)+genSpMI m = do+ mm <- genSpM m m+ return $ mm ^+^ eye m++ +++++-- | An arbitrary Householder reflection matrix +newtype PropHhReflMat a = PropHhReflMat { unHh :: SpMatrix a} deriving Show+instance Arbitrary (PropHhReflMat Double) where+ arbitrary = sized (\n -> PropHhReflMat <$> genReflMat n) `suchThat` ((> 2) . nrows . unHh)++genReflMat :: Int -> Gen (SpMatrix Double)+genReflMat n = do+ v <- normalize2 <$> (genSpVDense n :: Gen (SpVector Double))+ return $ hhRefl v+++++-- a product of a "large" number of random Householder matrices+newtype PropMatUnitary a = PropMatUnitary {unUnitary :: SpMatrix a} deriving Show+instance Arbitrary (PropMatUnitary Double) where+ arbitrary = sized (\n -> PropMatUnitary <$> genReflMat n) `suchThat` ((> 2) . nrows . unUnitary)+ +genUnitary :: Int -> Gen (SpMatrix Double)+genUnitary n = do+ q <- vectorOf (10 * n) $ genReflMat n -- random Householder matrix+ return $ foldr (##) (eye n) q ++prop_unitary :: (MatrixRing (SpMatrix a), Epsilon a, Eq a) =>+ PropMatUnitary a -> Bool+prop_unitary (PropMatUnitary m) = isOrthogonalSM m+++++-- | A symmetric, positive-definite matrix+newtype PropMatSPD a = PropMatSPD {unPropMatSPD :: SpMatrix a} deriving (Show)+instance Arbitrary (PropMatSPD Double) where+ arbitrary = sized genSpM_SPD `suchThat` ((> 2) . nrows . unPropMatSPD)++genSpM_SPD :: Int -> Gen (PropMatSPD Double)+genSpM_SPD n = do+ q <- genUnitary n + d <- genSpMDiagonal (all (> 0)) n -- positive diagonal+ return $ PropMatSPD ( q ## (d ##^ q) )+++++++-- | A pair of arbitrary SpMatrix, having compliant dimensions+data PropMatMat a = PropMatMat (SpMatrix a) (SpMatrix a) deriving (Eq, Show)+instance Arbitrary (PropMatMat Double) where+ arbitrary = sized3 genf where+ genf m n o = do+ mat1 <- genSpM m n+ mat2 <- genSpM n o+ return $ PropMatMat mat1 mat2++++-- | A square matrix and vector of compatible size+data PropMatVec a = PropMatVec (SpMatrix a) (SpVector a) deriving (Eq, Show)+instance Arbitrary (PropMatVec Double) where+ arbitrary = sized genf `suchThat` \(PropMatVec _ v) -> dim v > 2 where+ genf n = do+ mm <- genSpM n n+ v <- genSpV n+ return $ PropMatVec mm v++++-- -- | A symmetric positive definite matrix and vector of compatible size+-- data PropMatSPDVec a = PropMatSPDVec (SpMatrix a) (SpVector a) deriving (Eq, Show)+-- instance Arbitrary (PropMatSPDVec Double) where+-- arbitrary = do+-- PropMatVec m v <- arbitrary -- :: Gen (PropMatVec Double)+-- return $ PropMatSPDVec (m #^# m) v+++ ++++++++++++-- | QuickCheck properties++-- | Dot product of a normalized vector with itself is ~= 1+prop_dot :: (Normed v, Epsilon (Scalar v)) => v -> Bool+prop_dot v = let v' = normalize2 v in nearOne (v' <.> v')++-- | Positive semidefinite matrix. +prop_spd :: (LinearVectorSpace v, InnerSpace v, Ord (Scalar v), Num (Scalar v)) =>+ MatrixType v -> v -> Bool+prop_spd mm v = (v <.> (mm #> v)) >= 0++-- prop_spd' :: PropMatSPDVec Double -> Bool+-- prop_spd' (PropMatSPDVec m v) = prop_spd m v++++-- | (A B)^T == (B^T A^T)+prop_matMat1 :: (MatrixRing (SpMatrix t), Eq t) => PropMatMat t -> Bool+prop_matMat1 (PropMatMat a b) =+ transpose (a ## b) == (transpose b ##^ a)++-- | Implementation of transpose, (##), (##^) and (#^#) is consistent+prop_matMat2 :: (MatrixRing (SpMatrix t), Eq t) => PropMat t -> Bool+prop_matMat2 (PropMat m) = transpose m ##^ m == m #^# transpose m++++-- -- | The composition of a large number of random Householder reflection operator is an orthogonal matrix+-- prop_refl_orthog :: (MatrixRing (SpMatrix a), Epsilon a, Eq a) => PropHhReflMat a -> Bool+-- prop_refl_orthog (PropMatSPD m) = isOrthogonalSM m++++++-- | Cholesky factorization of a random SPD matrix +prop_Cholesky :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+ PropMatSPD a -> m Bool+prop_Cholesky (PropMatSPD m) = checkChol m+++-- | QR decomposition+prop_QR :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+ PropMatI a -> m Bool+prop_QR (PropMatI m) = checkQr m+++-- -- | check a random linear system+-- prop_linSolve :: LinSolveMethod -> PropMatVec Double -> Bool+-- prop_linSolve method (PropMatVec aa x) = do+-- let+-- aai = aa ^+^ eye (nrows aa) -- for invertibility+-- b = aai #> x+-- checkLinSolveR method aai b x++-- -- test data++++ {- example 0 : 2x2 linear system@@ -175,7 +634,6 @@ [1 2] [2] = [8] [3 4] [3] [18] - [1 3] [2] = [11] [2 4] [3] [16] @@ -186,19 +644,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]+b0, x0, x0true, aa0tx0 :: SpVector Double+b0 = mkSpVR 2 [8,18]+x0 = mkSpVR 2 [0.3,1.4] -- x0true : true solution-x0true = mkSpVectorD 2 [2,3]--+x0true = mkSpVR 2 [2,3] -aa0tx0 = mkSpVectorD 2 [11,16]+aa0tx0 = mkSpVR 2 [11,16] @@ -212,9 +667,9 @@ 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]+x1 = mkSpVR 4 [1,2,3,4] -b1 = mkSpVectorD 4 [30,56,60,101]+b1 = mkSpVR 4 [30,56,60,101] @@ -222,74 +677,43 @@ 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]+x2 = mkSpVR 3 [3,2,3] -b2 = mkSpVectorD 3 [4,-2,4]+b2 = mkSpVR 3 [4,-2,4] aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double --- -- -{--example 1 : random linear system --} +{- 2x2 Complex system -} +aa0c :: SpMatrix (Complex Double)+aa0c = fromListDenseSM 2 [ 3 :+ 1, (-3) :+ 2, (-2) :+ (-1), 1 :+ (-2)] --- -- 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)+b0c = mkSpVC 2 [3 :+ (-4), (-1) :+ 0.5] --- -- 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+x1c = mkSpVC 2 [2 :+ 2, 2 :+ 3]+b1c = mkSpVC 2 [4 :+ (-2), (-10) :+ 1] +aa2c :: SpMatrix (Complex Double)+aa2c = fromListDenseSM 2 [3, -3, -2, 1] --- 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 ::--- (Additive f1, Applicative f, Num a) => f (f1 a) -> f (f1 a) -> f (f1 a) --- plusM = liftA2 (^+^) -+-- matlab : aa = [1, 2-j; 2+j, 1-j]+aa3c, aa3cx :: SpMatrix (Complex Double)+aa3c = fromListDenseSM 2 [1, 2 :+ 1, 2 :+ (-1), 1 :+ (-1)] +-- matlab : aaxaa = aa * aa+aa3cx = fromListDenseSM 2 [6, 5, 3 :+ (-4), 5:+ (-2)] @@ -297,8 +721,6 @@ ---- {- matMat @@ -306,12 +728,18 @@ [3, 4] [7, 8] [43, 50] -} +m1, m2, m1m2, m1', m2', m1m2', m2m1' :: SpMatrix Double m1 = fromListDenseSM 2 [1,3,2,4] m2 = fromListDenseSM 2 [5, 7, 6, 8] m1m2 = fromListDenseSM 2 [19, 43, 22, 50] --- transposeSM+m1' = fromListSM (2,3) [(0,0,2), (1,0,3), (1,2,4), (1,2,1)]+m2' = fromListSM (3,2) [(0,0,5), (0,1,3), (2,1,4)]+m1m2' = fromListDenseSM 2 [10,15,6,13] +m2m1' = fromListSM (3,3) [(0,0,19),(2,0,12),(0,2,3),(2,2,4)] +-- transposeSM+m1t :: SpMatrix Double m1t = fromListDenseSM 2 [1,2,3,4] @@ -320,52 +748,48 @@ {- countSubdiagonalNZ -}-+m3 :: SpMatrix Double 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+aa3 :: SpMatrix Double+aa3 = fromListDenseSM 3 [1,1,3,2,2,2,3,1,1] -b3 = mkSpVectorD 3 [1,1,1] :: SpVector Double+b3 = mkSpVR 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+aa4 :: SpMatrix Double+aa4 = fromListDenseSM 3 [3,2,-2,2,2,-1,6,5,-4] +aa4c :: SpMatrix (Complex Double)+aa4c = toC <$> aa4+ b4 = fromListDenseSV 3 [-3,-3,-3] :: SpVector Double +aa5c :: SpMatrix (Complex Double)+aa5c = fromListDenseSM 4 cv where+ cv = zipWith (:+) [1..16] [16,15..1] --- test data -tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double++tm0, tm1, tm2, tm3, tm4, tm5, tm6 :: 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]+tv0 = mkSpVR 2 [5, 6] tv1 = fromListSV 2 [(0,1)] @@ -373,27 +797,12 @@ 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 --@@ -401,10 +810,10 @@ 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+tm5 = fromListDenseSM 3 [2, -4, -4, -1, 6, -2, -2, 3, 8] -tm6 = fromListDenseSM 4 [1,3,4,2,2,5,2,10,3,6,8,11,4,7,9,12] :: SpMatrix Double+tm6 = fromListDenseSM 4 [1,3,4,2,2,5,2,10,3,6,8,11,4,7,9,12] tm7 :: SpMatrix Double tm7 = a ^+^ b ^+^ c where@@ -413,11 +822,90 @@ b = mkSubDiagonal n 0 $ replicate n 2 c = mkSubDiagonal n (-1) $ replicate n (-1) +tvx7 = mkSpVR 5 [3,8,-12,4,9] +tvb7 = tm7 #> tvx7 +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)]++++tm9 :: SpMatrix Double+tm9 = fromListSM (4, 3) [(0,0,pi), (1,1, 3), (2,2,4), (3,2, 1), (3,1, 5)]++++++-- tvc0 <.> tvc1 = 5 +tvc0, tvc1, tvc2, tvc3 :: SpVector (Complex Double)+tvc0 = fromListSV 2 [(0,0), (1,2 :+ 1)]+tvc1 = fromListSV 2 [(0,0), (1, 2 :+ (-1))] +++-- dot([1+i, 2-i], [3-2i, 1+i]) = 2 - 2i+tvc2 = fromListDenseSV 2 [1 :+ 1, 2 :+ (-1)]+tvc3 = fromListDenseSV 2 [3 :+ (-2), 1 :+ 1 ]++++++--++-- | Example 5.4.2 from G & VL+-- aa1 :: SpMatrix Double+-- aa1 = transpose $ fromListDenseSM 3 [1..12]++-- aa1 :: SpMatrix Double+-- aa1 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]+++++++-- l0 = [1,2,4,5,8]+-- l1 = [2,3,6]+-- l2 = [7]++-- v0,v1 :: V.Vector Int+-- v0 = V.fromList [0,1,2,5,6]+-- v1 = V.fromList [0,3,4,6]++-- -- e1, e2 :: V.Vector (Int, Double)+-- -- e1 = V.indexed $ V.fromList [1,0,0]+-- -- e2 = V.indexed $ V.fromList [0,1,0]++-- e1, e2:: CsrVector Double+-- e1 = fromListCV 4 [(0, 1)] +-- e2 = fromListCV 4 [(1, 1)]+-- e3 = fromListCV 4 [(0, 1 :+ 2)] :: CsrVector (Complex Double)++-- e1c = V.indexed $ V.fromList [1,0,0] :: V.Vector (Int, Complex Double)++-- m0,m1,m2,m3 :: CsrMatrix Double+-- m0 = toCSR 2 2 $ V.fromList [(0,0, pi), (1,0,3), (1,1,2)]+-- m1 = toCSR 4 4 $ V.fromList [(0,0,1), (0,2,5), (1,0,2), (1,1,3), (2,0,4), (2,3,1), (3,2,2)]+-- m2 = toCSR 4 4 $ V.fromList [(0,0,1), (0,2,5), (2,0,4), (2,3,1), (3,2,2)]+-- m3 = toCSR 4 4 $ V.fromList [(1,0,5), (1,1,8), (2,2,3), (3,1,6)]++++++++-- --++ -- -- run N iterations -- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB@@ -450,13 +938,50 @@ -- 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)]+{-+random linear system +-} -tm9 :: SpMatrix Double-tm9 = fromListSM (4, 3) [(0,0,pi), (1,1, 3), (2,2,4), (3,2, 1), (3,1, 5)]++-- -- 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