sparse-linear-algebra 0.1.0.0 → 0.1.0.1

raw patch · 4 files changed

+518/−297 lines, 4 filesdep −monad-loopsdep −transformersdep ~basePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies removed: monad-loops, transformers

Dependency ranges changed: base

API changes (from Hackage documentation)

Files

• app/Main.hs +0/−8
• sparse-linear-algebra.cabal +17/−19
• src/Math/Linear/Sparse.hs +181/−270
• test/LibSpec.hs +320/−0

@@ -1,8 +0,0 @@-module Main where--import Lib (ourAdd)--import Text.Printf (printf)--main :: IO ()-main = printf "2 + 3 = %d\n" (ourAdd 2 3)

@@ -1,14 +1,14 @@ name:                sparse-linear-algebra-version:             0.1.0.0-synopsis:            Sparse linear algebra datastructures and algorithms+version:             0.1.0.1+synopsis:            Sparse linear algebra datastructures and algorithms. Currently it provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. description:         Please see README.md homepage:            https://github.com/ocramz/sparse-linear-algebra-license:             BSD3+license:             GPL-3 license-file:        LICENSE author:              Marco Zocca maintainer:          zocca.marco gmail copyright:           2016 Marco Zocca-category:            Math+category:            Numeric build-type:          Simple extra-source-files:  README.md cabal-version:       >=1.10@@ -26,29 +26,28 @@                      , containers                      , hspec                      , primitive >= 0.6.1.0-                     , transformers >= 0.5.2.0-                     -- , lens                      , mtl >= 2.2.1                      , mwc-random-                     , monad-loops -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-                     , mtl >= 2.2.1-                     , mwc-random-                     , primitive >= 0.6.1.0-                     , sparse-linear-algebra-                     , transformers >= 0.5.2.0 +-- 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+--                      , mtl >= 2.2.1+--                      , mwc-random+--                      , primitive >= 0.6.1.0+--                      , sparse-linear-algebra+--                      , transformers >= 0.5.2.0+ test-suite spec   default-language:    Haskell2010   ghc-options:         -Wall   type:                exitcode-stdio-1.0   hs-source-dirs:      test+  other-modules:       LibSpec   main-is:             Spec.hs   build-depends:       base                      , containers@@ -57,7 +56,6 @@                      , mwc-random                      , primitive >= 0.6.1.0                      , sparse-linear-algebra-                     , transformers >= 0.5.2.0                      , criterion                      -- , QuickCheck 

@@ -9,15 +9,13 @@ import Control.Monad.Primitive  import Control.Monad (mapM_, forM_, replicateM)-import Control.Monad.Loops -import Control.Monad.Cont import Control.Monad.State.Strict-import Control.Monad.Writer-import Control.Monad.Trans+-- import Control.Monad.Writer+-- import Control.Monad.Trans -import Control.Monad.Trans.State (runStateT)-import Control.Monad.Trans.Writer (runWriterT)+-- 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`@@ -30,33 +28,35 @@ import qualified Data.Traversable as T  -import qualified Data.List as L+-- import qualified Data.List as L import Data.Maybe   --- | ========= CLASSES and common operations+{-|  CLASSES and common operations -} --- | Additive ring +-- * Additive ring  class Functor f => Additive f where-  -- | zero element+  -- | Ring zero element   zero :: Num a => f a   -  -- | componentwise operations+  -- | Ring +   (^+^) :: Num a => f a -> f a -> f a-  (^-^) :: Num a => f a -> f a -> f a   + -- | negate the values in a functor negated :: (Num a, Functor f) => f a -> f a negated = fmap negate -x `minus` y = x ^+^ negated y+-- | subtract two Additive objects+(^-^) :: (Additive f, Num a) => f a -> f a -> f a+x ^-^ y = x ^+^ negated y   --- | Vector space+-- * Vector space class Additive f => VectorSpace f where   -- | multiplication by a scalar   (.*) :: Num a => a -> f a -> f a@@ -67,44 +67,45 @@ lerp a u v = a .* u ^+^ ((1-a) .* v)  --- | Hilbert space (inner product)+-- * Hilbert space (inner product) class VectorSpace f => Hilbert f where   -- | inner product   dot :: Num a => f a -> f a -> a +-- * Normed vector space class Hilbert f => Normed f where   norm :: (Floating a, Eq a) => a -> f a -> a  --- some norms and related results+-- ** Norms and related results --- squared norm +-- *** squared 2-norm normSq :: (Hilbert f, Num a) => f a -> a normSq v = v `dot` v  --- L1 norm+-- *** L1 norm norm1 :: (Foldable t, Num a, Functor t) => t a -> a norm1 v = sum (fmap abs v) --- Euclidean norm+-- *** Euclidean norm norm2 :: (Hilbert f, Floating a) => f a -> a norm2 v = sqrt (normSq v) --- Lp norm (p > 0)+-- *** 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+-- *** infinity-norm normInfty :: (Foldable t, Ord a) => t a -> a normInfty = maximum   --- normalize+-- *** normalize w.r.t. p-norm (p finite) normalize :: (Normed f, Floating a, Eq a) => a -> f a -> f a-normalize n v = (1 / norm n v) .* v+normalize p v = (1 / norm p v) .* v   @@ -112,19 +113,19 @@   --- -- 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  --- scale+-- *** scale scale :: (Num b, Functor f) => b -> f b -> f b scale n = fmap (* n) @@ -134,24 +135,20 @@   --- | FiniteDim : finite-dimensional objects+-- ** FiniteDim : finite-dimensional objects  class Additive f => FiniteDim f where   type FDSize f :: *   dim :: f a -> FDSize f  -instance FiniteDim SpVector where-  type FDSize SpVector = Int-  dim = svDim  -instance FiniteDim SpMatrix where-  type FDSize SpMatrix = (Rows, Cols)-  dim = smDim  --- unary dimension-checking bracket+++-- | unary dimension-checking bracket withDim :: (FiniteDim f, Show e) =>      f a      -> (FDSize f -> f a -> Bool)@@ -162,7 +159,7 @@ 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+-- | binary dimension-checking bracket withDim2 :: (FiniteDim f, FiniteDim g, Show e) =>      f a      -> g b@@ -179,40 +176,28 @@   --- | HasData : accessing inner data (do not export)+-- ** HasData : accessing inner data (do not export)  class Additive f => HasData f a where   type HDData f a :: *    dat :: f a -> HDData f a -instance HasData SpVector a where-  type HDData SpVector a = IM.IntMap a-  dat = svData -instance HasData SpMatrix a where-  type HDData SpMatrix a = IM.IntMap (IM.IntMap a)-  dat = smData  ------- | Sparse : sparse datastructures+-- ** Sparse : sparse datastructures  class (FiniteDim f, HasData f a) => Sparse f a where   spy :: Fractional b => f a -> b  -instance Sparse SpVector a where-  spy = spySV -instance Sparse SpMatrix a where-  spy = spySM    +-- ** Set : things that behave as sets (e.g. of which we can take the union and the intersection)+ class Functor f => Set f where   -- |union binary lift   liftU2 :: (a -> a -> a) -> f a -> f a -> f a@@ -228,7 +213,7 @@   --- | =======================================================+-- * IntMap implementation  instance Set IM.IntMap where   liftU2 = IM.unionWith@@ -236,15 +221,13 @@   liftI2 = IM.intersectionWith   {-# INLINE liftI2 #-} --- | IntMap implementation instance Additive IM.IntMap where   zero = IM.empty   {-# INLINE zero #-}   (^+^) = liftU2 (+)   {-# INLINE (^+^) #-}-  x ^-^ y = x ^+^ negated y-  {-# INLINE (^-^) #-} + instance VectorSpace IM.IntMap where   n .* im = IM.map (* n) im   @@ -259,27 +242,22 @@   --- | =======================================================+-- * Sparse Vector --- | Sparse Vector data SpVector a = SV { svDim :: Int ,                        svData :: IM.IntMap a} deriving Eq -dimSV :: SpVector a -> Int-dimSV = svDim-+-- | SpVector sparsity spySV :: Fractional b => SpVector a -> b-spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (svDim s)+spySV s = fromIntegral (IM.size (dat s)) / fromIntegral (dim s)  --- internal : projection functions, do not export-imSV :: SpVector a -> IM.IntMap a-imSV = svData    --- | instances for SpVector++-- ** instances for SpVector instance Functor SpVector where   fmap f (SV n x) = SV n (fmap f x) @@ -293,12 +271,25 @@ instance Additive SpVector where   zero = SV 0 IM.empty   (^+^) = liftU2 (+)-  (^-^) = liftU2 (-) +                        instance VectorSpace SpVector where   n .* v = scale n v ++instance FiniteDim SpVector where+  type FDSize SpVector = Int+  dim = svDim  ++instance HasData SpVector a where+  type HDData SpVector a = IM.IntMap a+  dat = svData++instance Sparse SpVector a where+  spy = spySV++ instance Hilbert SpVector where   a `dot` b | dim a == dim b = dot (dat a) (dat b)             | otherwise =@@ -323,6 +314,7 @@ singletonSV :: a -> SpVector a singletonSV x = SV 1 (IM.singleton 0 x) +-- ** Create new sparse vector  -- | create a sparse vector from an association list while discarding all zero entries mkSpVector :: (Num a, Eq a) => Int -> IM.IntMap a -> SpVector a@@ -336,8 +328,9 @@ mkSpVector1 :: Int -> IM.IntMap a -> SpVector a mkSpVector1 d ll = SV d $ IM.filterWithKey (\ k _ -> inBounds0 d k) ll -mkSpVector1D :: Int -> [a] -> SpVector a-mkSpVector1D d ll = mkSpVector1 d (IM.fromList $ denseIxArray (take d 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))   @@ -352,7 +345,7 @@   --- insert+-- |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)@@ -362,11 +355,11 @@ fromListSV :: Int -> [(Int, a)] -> SpVector a fromListSV d iix = SV d (IM.fromList (filter (inBounds0 d . fst) iix )) --- toList+-- |toList toListSV :: SpVector a -> [(IM.Key, a)]-toListSV sv = IM.toList (imSV sv)+toListSV sv = IM.toList (dat sv) --- to dense list (default = 0)+-- |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] @@ -396,13 +389,14 @@  -- | SV manipulation +-- | Tail elements tailSV :: SpVector a -> SpVector a tailSV (SV n sv) = SV (n-1) 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 (imSV sv))+headSV sv = fromMaybe 0 (IM.lookup 0 (dat sv))   @@ -431,7 +425,7 @@       --- | outer vector product+-- *** Outer vector product  outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a outerProdSV v1 v2 = fromListSM (m, n) ixy where@@ -449,17 +443,14 @@   --- | =======================================================---+-- * Sparse Matrix  data SpMatrix a = SM {smDim :: (Rows, Cols),                       smData :: IM.IntMap (IM.IntMap a)} deriving Eq   --- | instances for SpMatrix+-- ** Instances for SpMatrix instance Show a => Show (SpMatrix a) where   show sm@(SM _ x) = "SM " ++ sizeStr sm ++ " "++ show (IM.toList x) @@ -473,9 +464,20 @@ instance Additive SpMatrix where   zero = SM (0,0) IM.empty   (^+^) = liftU2 (+)-  (^-^) = liftU2 (-)  +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)+  dat = smData++instance Sparse SpMatrix a where+  spy = spySM+  + -- | TODO : use semilattice properties instead maxTup, minTup :: Ord t => (t, t) -> (t, t) -> (t, t) maxTup (x1,y1) (x2,y2) = (max x1 x2, max y1 y2)@@ -487,11 +489,11 @@   --- multiply matrix by a scalar+-- *** multiply matrix by a scalar matScale :: Num a => a -> SpMatrix a -> SpMatrix a matScale a = fmap (*a) --- Frobenius norm (sqrt of trace of M^T M)+-- *** Frobenius norm (sqrt of trace of M^T M) normFrobenius :: SpMatrix Double -> Double normFrobenius m = sqrt $ foldlSM (+) 0 m' where   m' | nrows m > ncols m = transposeSM m ## m@@ -502,7 +504,7 @@   --- | ========= MATRIX METADATA+-- ** MATRIX METADATA  -- type synonyms type Rows = Int@@ -511,16 +513,16 @@ type IxRow = Int type IxCol = Int --- -- predicates--- are the supplied indices within matrix bounds?+-- *** predicates+-- |Are the supplied indices within matrix bounds? validIxSM :: SpMatrix a -> (Int, Int) -> Bool validIxSM mm = inBounds02 (dim mm) --- is the matrix square?+-- |Is the matrix square? isSquareSM :: SpMatrix a -> Bool isSquareSM m = nrows m == ncols m --- is the matrix diagonal?+-- |Is the matrix diagonal? isDiagonalSM :: SpMatrix a -> Bool isDiagonalSM m = IM.size d == nrows m where   d = IM.filterWithKey ff (immSM m)@@ -565,7 +567,7 @@ spySM s = fromIntegral (nzSM s) / fromIntegral (nelSM s)  --- # NZ in row i+-- *** Non-zero elements in a row  nzRowU :: SpMatrix a -> IM.Key -> Int nzRowU s i = maybe 0 IM.size (IM.lookup i $ immSM s)@@ -577,7 +579,7 @@   --- | bandwidth bounds (min, max)+-- *** bandwidth bounds (min, max)  bwMinSM :: SpMatrix a -> Int bwMinSM = fst . bwBoundsSM@@ -603,7 +605,7 @@   --- | ========= SPARSE MATRIX BUILDERS+-- ** SPARSE MATRIX BUILDERS  zeroSM :: Int -> Int -> SpMatrix a zeroSM m n = SM (m,n) IM.empty @@ -617,15 +619,17 @@       d = dim s  --- | from list (row, col, value)+-- | 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   ins t (i,j,x) = insertSpMatrix i j x t +-- | Create new SpMatrix using data from list (row, col, value) fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a fromListSM (m,n) iix = fromListSM' iix (zeroSM m n)  +-- | 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   n = length ll `div` m@@ -634,7 +638,7 @@  -- | to List --- toDenseListSM : populate missing entries with 0+-- |Populate missing entries with 0 toDenseListSM :: Num t => SpMatrix t -> [(IxRow, IxCol, t)] toDenseListSM m =   [(i, j, m @@ (i, j)) | i <- [0 .. nrows m - 1], j <- [0 .. ncols m- 1]]@@ -643,11 +647,11 @@   --- -- create diagonal and identity matrix+-- ** Diagonal matrix mkDiagonal :: Int -> [a] -> SpMatrix a mkDiagonal n = mkSubDiagonal n 0 -+-- *** Identity matrix eye :: Num a => Int -> SpMatrix a eye n = mkDiagonal n (ones n) @@ -658,7 +662,7 @@     --- super- and sub- diagonal+-- *** Create Super- or sub- diagonal matrix  mkSubDiagonal :: Int -> Int -> [a] -> SpMatrix a mkSubDiagonal n o xx | abs o < n = if o >= 0@@ -693,7 +697,7 @@   --- | ========= SUB-MATRICES+-- ** SUB-MATRICES   extractSubmatrixSM :: SpMatrix a -> (Int, Int) -> (Int, Int) -> SpMatrix a@@ -714,16 +718,7 @@       inBounds0 c j2 &&             i2 >= i1 --- extract row / column-extractRowSM :: SpMatrix a -> Int -> SpMatrix a-extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)--extractColSM :: SpMatrix a -> Int -> SpMatrix a-extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)------ demote (n x 1) or (1 x n) SpMatrix to SpVector+-- |Demote (n x 1) or (1 x n) SpMatrix to SpVector toSV :: SpMatrix a -> SpVector a toSV (SM (m,n) im) = SV d $ snd . head $ IM.toList im where   d | m==1 && n==1 = 1@@ -732,12 +727,22 @@     | otherwise = error $ "toSV : incompatible dimensions " ++ show (m,n)  --- extract row or column and place into SpVector+-- *** Extract jth column+extractColSM :: SpMatrix a -> Int -> SpMatrix a+extractColSM sm j = extractSubmatrixSM sm (0, nrows sm - 1) (j, j)++-- |", and place into SpVector extractCol :: SpMatrix a -> Int -> SpVector a-extractCol m i = toSV $ extractColSM m i+extractCol m j = toSV $ extractColSM m j ++-- *** Extract ith row+extractRowSM :: SpMatrix a -> Int -> SpMatrix a+extractRowSM sm i = extractSubmatrixSM sm (i, i) (0, ncols sm - 1)++-- |", and place into SpVector extractRow :: SpMatrix a -> Int -> SpVector a-extractRow m j = toSV $ extractRowSM m j+extractRow m i = toSV $ extractRowSM m i   @@ -745,8 +750,9 @@   --- | ========= MATRIX STACKING+-- ** MATRIX STACKING +-- | Vertical stacking vertStackSM, (-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a vertStackSM mm1 mm2 = SM (m, n) $ IM.union u1 u2 where   nro1 = nrows mm1@@ -757,7 +763,7 @@  (-=-) = vertStackSM -+-- | Horizontal stacking horizStackSM, (-||-) :: SpMatrix a -> SpMatrix a -> SpMatrix a horizStackSM mm1 mm2 = t (t mm1 -=- t mm2) where   t = transposeSM@@ -772,12 +778,12 @@   --- | ========= LOOKUP+-- ** MATRIX ELEMENT LOOKUP  lookupSM :: SpMatrix a -> IM.Key -> IM.Key -> Maybe a lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j --- | Looks up an element in the matrix (if not found, zero is returned)+-- | Looks up an element in the matrix with a default (if the element is not found, zero is returned)  lookupWD_SM, (@@) :: Num a => SpMatrix a -> (IM.Key, IM.Key) -> a lookupWD_SM sm (i,j) =@@ -786,6 +792,7 @@ lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IM.Key, IM.Key) -> a lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j) +-- | Zero-default lookup, infix form (@@) = lookupWD_SM  @@ -799,7 +806,7 @@   --- | ========= MISC SpMatrix OPERATIONS+-- *** MISC SpMatrix OPERATIONS  foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b foldlSM f n (SM _ m)= foldlIM2 f n m@@ -830,9 +837,9 @@ -- extractDiagonalSM :: (Num a, Eq a) => SpMatrix a -> SpVector a -- extractDiagonalSM (SM (m,n) im) = mkSpVectorD m $ extractDiagonalIM2 im --- extract with default 0+-- | extract with default 0 extractDiagonalDSM :: Num a => SpMatrix a -> SpVector a-extractDiagonalDSM mm = mkSpVector1D n $ foldr ins [] ll  where+extractDiagonalDSM mm = fromListDenseSV n $ foldr ins [] ll  where   ll = [0 .. n - 1]   n = nrows mm   ins i acc = mm@@(i,i) : acc@@ -852,7 +859,7 @@   --- | sparsify : remove 0s (!!!)+-- ** sparsify : remove 0s (!!!)  sparsifyIM2 :: IM.IntMap (IM.IntMap Double) -> IM.IntMap (IM.IntMap Double) sparsifyIM2 = ifilterIM2 (\_ _ x -> abs x >= eps)@@ -862,7 +869,7 @@   --- | ROUNDING operations (!!!)+-- ** ROUNDING operations (!!!)                                roundZeroOneSM :: SpMatrix Double -> SpMatrix Double roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im@@ -877,10 +884,10 @@   --- | ========= ALGEBRAIC PRIMITIVE OPERATIONS+-- * ALGEBRAIC PRIMITIVE OPERATIONS  --- | transpose+-- ** Matrix transpose   transposeSM, (#^) :: SpMatrix a -> SpMatrix a@@ -904,7 +911,7 @@   --- | matrix action on a vector+-- ** matrix action on a vector  {-  FIXME : matVec is more general than SpVector's :@@ -915,7 +922,7 @@   --- matrix on vector+-- |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@@ -923,7 +930,7 @@  (#>) = matVec --- vector on matrix (FIXME : transposes matrix: more costly than `matVec`)+-- |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)   | n == nr = SV nc $ fmap (`dot` sv) (transposeIM2 mdata)@@ -940,17 +947,12 @@   --- matVec' mm vv =---   withDim2 mm vv (\(nro, nco) nv _ _ -> nco == nv) matVecU "matVec : mismatching dimensions"---    (\ m v -> unwords [show (dim m), show (dim v)])  --- asdfm ll = unwords (map (show . dim) ll)   ---- | matrix-matrix product+-- ** Matrix-matrix product  -- unsafe matMat matMatU :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a@@ -959,12 +961,6 @@     im = fmap (\vm1 -> (`dot` vm1) <$> transposeIM2 (immSM m2)) (immSM m1)  --- matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a--- matMat (SM (nr1,nc1) m1) (SM (nr2,nc2) m2)---   | nc1 == nr2 = SM (nr1, nc2) $---       fmap (\vm1 -> fmap (`dot` vm1) (transposeIM2 m2)) m1---   | otherwise = error "matMat : incompatible matrix sizes"- matMat, (##) :: Num a => SpMatrix a -> SpMatrix a -> SpMatrix a matMat m1 m2   | c1 == r2 = matMatU m1 m2@@ -986,7 +982,7 @@   --- | sparsified matrix-matrix product (prunes all elements `x` for which `abs x <= eps`)+-- ** Matrix-matrix product, sparsified (prunes all elements `x` for which `abs x <= eps`) matMatSparsified, (#~#)  :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2 @@ -997,9 +993,9 @@   --- | ========= predicates+-- * Predicates --- is the matrix orthogonal? i.e. Q^t ## Q == I+-- |is the matrix orthogonal? i.e. Q^t ## Q == I isOrthogonalSM :: SpMatrix Double -> Bool isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where   rsm = roundZeroOneSM $ transposeSM sm ## sm@@ -1012,9 +1008,9 @@   --- | ========= condition number+-- ** Condition number --- uses the R matrix from the QR factorization+-- |uses the R matrix from the QR factorization conditionNumberSM :: SpMatrix Double -> Double conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system"                     | otherwise = kappa where@@ -1030,15 +1026,16 @@   --- | ========= Householder transformation+-- ** 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 +{-| 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 :: SpVector Double -> SpMatrix Double hhRefl = hhMat 2.0 @@ -1052,7 +1049,7 @@   --- | ========= Givens rotation matrix+-- ** Givens rotation matrix   hypot :: Floating a => a -> a -> a@@ -1064,6 +1061,7 @@   | 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)@@ -1076,7 +1074,7 @@                 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@@ -1096,12 +1094,12 @@     a = mm @@ (i', j)     b = mm @@ (i, j)   -- element to zero out --- is the `k`th the first nonzero column in the row?+-- |Is the `k`th the first nonzero column in the row? firstNonZeroColumn :: IM.IntMap a -> IM.Key -> Bool firstNonZeroColumn mm k = isJust (IM.lookup k mm) &&                           isNothing (IM.lookupLT k mm) --- returns a set of rows {k} that satisfy QR.C1+-- |Returns a set of rows {k} that satisfy QR.C1 candidateRows :: IM.IntMap (IM.IntMap a) -> IM.Key -> IM.Key -> Maybe [IM.Key] candidateRows mm i j | IM.null u = Nothing                      | otherwise = Just (IM.keys u) where@@ -1112,20 +1110,17 @@   --- | ========= QR algorithm--{--applies Givens rotation iteratively to zero out sub-diagonal elements--}+-- ** QR decomposition  +-- | Applies Givens rotation iteratively to zero out sub-diagonal elements qr :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double) qr mm = (transposeSM qmatt, rmat)  where   qmatt = F.foldl' (#~#) ee $ gmats mm -- Q^T = (G_n * G_n-1 ... * G_1)   rmat = qmatt #~# mm                  -- R = Q^T A   ee = eye (nrows mm)       --- Givens matrices in order [G1, G2, .. , G_N ]+-- | Givens matrices in order [G1, G2, .. , G_N ] gmats :: SpMatrix Double -> [SpMatrix Double] gmats mm = gm mm (subdiagIndicesSM mm) where  gm m ((i,j):is) = let g = givens m i j@@ -1148,9 +1143,12 @@   --- | ========= Eigenvalues, using QR +-- ** Eigenvalue algorithms +-- *** All eigenvalues using QR algorithm++ eigsQR :: Int -> SpMatrix Double -> SpVector Double eigsQR nitermax m = extractDiagonalDSM $ execState (convergtest eigsStep) m where   eigsStep m = r #~# q where (q, r) = qr m@@ -1164,11 +1162,12 @@   --- | ========= Eigenvalues, using Rayleigh iteration+-- *** One eigenvalue and corresponding eigenvector, using Rayleigh iteration  +-- | Cubic-order convergence, but it requires a mildly educated guess on the initial eigenpair rayleighStep ::   SpMatrix Double ->-  (SpVector Double, Double) ->+  (SpVector Double, Double) ->    (SpVector Double, Double)    -- updated estimate of (eigenvector, eigenvalue) rayleighStep aa (b, mu) = (b', mu') where   ii = eye (nrows aa)@@ -1188,9 +1187,8 @@   --- | ========= Householder vector (G & VL Alg. 5.1.1, function `house`)+-- ** Householder vector (G & VL Alg. 5.1.1, function `house`) --- hhV :: (Ord a, Floating a) => SpVector a -> (SpVector a, a) hhV :: SpVector Double -> (SpVector Double, Double) hhV x = (v, beta) where   n = dim x@@ -1212,7 +1210,7 @@   --- | ========= SVD+-- * SVD  {- Golub & Van Loan, sec 8.6.2 (p 452 segg.) @@ -1240,10 +1238,10 @@   --- | ======================================================= --- | LINEAR SOLVERS : solve A x = b +-- * LINEAR SOLVERS : solve A x = b+ -- | numerical tolerance for e.g. solution convergence eps :: Double eps = 1e-8@@ -1257,7 +1255,7 @@   --- | CGS+-- ** CGS  -- | one step of CGS cgsStep :: SpMatrix Double -> SpVector Double -> CGS -> CGS@@ -1277,15 +1275,14 @@                  _p :: SpVector Double,                  _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 ->-  -- Int ->   CGS cgs aa b x0 rhat =-  -- execState (replicateM n (modify (cgsStep aa rhat))) cgsInit where   execState (untilConverged _x (cgsStep aa rhat)) cgsInit where   r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution   p0 = r0@@ -1303,7 +1300,7 @@     --- | BiCSSTAB+-- ** BiCSSTAB  -- _aa :: SpMatrix Double,    -- matrix -- _b :: SpVector Double,     -- rhs@@ -1327,21 +1324,18 @@                            _rBicgstab :: SpVector Double,                            _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-     -- -> Int      -> BICGSTAB bicgstab aa b x0 r0hat =-  -- execState (replicateM n (modify (bicgstabStep aa r0hat))) bicgsInit where   execState (untilConverged _xBicgstab (bicgstabStep aa r0hat)) bicgsInit where    r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution    p0 = r0    bicgsInit = BICGSTAB x0 r0 p0-   -- q (BICGSTAB xi _ _) = nor  instance Show BICGSTAB where   show (BICGSTAB x r p) = "x = " ++ show x ++ "\n" ++@@ -1358,11 +1352,11 @@   --- | ========= LINEAR SOLVERS INTERFACE+-- * LINEAR SOLVERS INTERFACE  data LinSolveMethod = CGS_ | BICGSTAB_ deriving (Eq, Show)  --- random starting vector+-- | linear solve with _random_ starting vector linSolveM ::   PrimMonad m =>     LinSolveMethod -> SpMatrix Double -> SpVector Double -> m (SpVector Double)@@ -1374,7 +1368,7 @@     case method of CGS_ -> return $ _xBicgstab (bicgstab aa b x0 x0)                    BICGSTAB_ -> return $ _x (cgs aa b x0 x0) --- deterministic starting vector (every component at 0.1) +-- | linear solve with _deterministic_ starting vector (every component at 0.1)  linSolve ::   LinSolveMethod -> SpMatrix Double -> SpVector Double -> SpVector Double linSolve method aa b@@ -1422,9 +1416,8 @@   --- | ========= PRETTY PRINTING+-- * PRETTY PRINTING --- | Show details and contents of sparse matrix  sizeStr :: SpMatrix a -> String sizeStr sm =@@ -1484,6 +1477,7 @@       rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]            | otherwise = rr_ +-- ** Pretty printer typeclass class PrintDense a where   prd :: a -> IO () @@ -1550,9 +1544,13 @@                            else go (i + 1) (take 2 $ y:ll) y                 where y = f xx --- modify state and append, until max # of iterations is reached+-- | keep state `x` in a moving window of length 2 to assess convergence, stop when either a condition on that list is satisfied or when max # of iterations is reached modifyInspectN ::-  MonadState s m => Int -> ([s] -> Bool) -> (s -> s) -> m s+  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@@ -1569,13 +1567,14 @@                        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]  @@ -1594,9 +1593,12 @@ normDiffConverged :: (Foldable t, Functor t) =>      (a -> SpVector Double) -> t a -> Bool normDiffConverged fp xx = normSq (foldrMap fp (^-^) (zeroSV 0) xx) <= eps-                +  +++ -- 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]@@ -1616,8 +1618,6 @@     if n <= 0 then xs               else go f (n-1) x (x : xs) --- runN :: Int -> (a -> a) -> a -> a--- runN n stepf x0 = runAppendN' stepf n x0     @@ -1739,104 +1739,15 @@  -- -tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double-tm0 = fromListSM (2,2) [(0,0,pi), (1,0,sqrt 2), (0,1, exp 1), (1,1,sqrt 5)] -tv0, tv1 :: SpVector Double-tv0 = mkSpVectorD 2 [5, 6]  -tv1 = SV 2 $ IM.singleton 0 1 --- wikipedia test matrix for Givens rotation -tm1 = sparsifySM $ fromListDenseSM 3 [6,5,0,5,1,4,0,4,3] -tm1g1 = givens tm1 1 0-tm1a2 = tm1g1 ## tm1 -tm1g2 = givens tm1a2 2 1-tm1a3 = tm1g2 ## tm1a2 -tm1q = transposeSM (tm1g2 ## tm1g1) ---- wp test matrix for QR decomposition via Givens rotation--tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41]-----tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9]--tm3g1 = fromListDenseSM 3 [1, 0,0, 0,c,-s, 0, s, c]-  where c= 0.4961-        s = 0.8682-------tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]----- playground---- | terminate after n iterations or when q becomes true, whichever comes first-untilC :: (a -> Bool) -> Int ->  (a -> a) -> a -> a-untilC p n f = go n-  where-    go m x | p x || m <= 0 = x-           | otherwise     = True `seq` go (m-1) (f x)-------- testing State----- data T0 = T0 {unT :: Int} deriving Eq--- instance Show T0 where---   show (T0 x) = show x---- -- modifyT :: MonadState T0 m => (Int -> Int) -> m String--- modifyT f = state (\(T0 i) -> (i, T0 (f i)))-  ---- t00 = T0 0---- testT n = execState $ replicateM n (modifyT (+1)) ----- testT2 = execState $ when -  ---- replicateSwitch p m f = loop m where---       loop n | n <= 0 || p = pure (#)---              | otherwise = f *> loop (n-1)------ testing Writer-               --- asdfw n = runWriter $ do---   tell $ "potato " ++ show n---   tell "jam"---   return (n+1)----- ------- testing State and Writer------ runMyApp runA k maxDepth =---     let config = maxDepth---         state =  0---     in runStateT (runWriterT (runA k) config) state     

@@ -0,0 +1,320 @@+{-# language ScopedTypeVariables #-}+module LibSpec where++import qualified Data.IntMap as IM++import Control.Monad (replicateM)+import Control.Monad.State.Strict (execState)++import qualified System.Random.MWC as MWC+import qualified System.Random.MWC.Distributions as MWC+       +import Test.Hspec+-- import Test.Hspec.QuickCheck++import Lib+import Math.Linear.Sparse++++main :: IO ()+main = hspec spec++-- niter = 5++spec :: Spec+spec = do+  describe "Math.Linear.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" $+      normSq ((aa0 #> x0true) ^-^ b0 ) <= eps `shouldBe` True+    it "vecMat : vector-matrix product" $+      normSq ((x0true <# aa0) ^-^ aa0tx0 ) <= eps `shouldBe` True  +    it "matMat : matrix-matrix product" $+      (m1 `matMat` m2) `shouldBe` m1m2+    it "eye : identity matrix" $+      infoSM (eye 10) `shouldBe` SMInfo 10 0.1+    it "countSubdiagonalNZ : # of nonzero elements below the diagonal" $+      countSubdiagonalNZSM m3 `shouldBe` 1+    it "modifyInspectN : early termination by iteration count" $+      execState (modifyInspectN 2 ((< eps) . diffSqL) (/2)) 1 `shouldBe` 1/8+    it "modifyInspectN : termination by value convergence" $+      execState (modifyInspectN (2^16) ((< eps) . head) (/2)) 1 < eps `shouldBe` True +  describe "Math.Linear.Sparse : Linear solvers" $ do    +    it "BiCGSTAB (2 x 2 dense)" $ +      -- normSq (_xBicgstab (bicgstab aa0 b0 x0 x0) ^-^ x0true) <= eps `shouldBe` True+      normSq (aa0 <\> b0 ^-^ x0true) <= eps `shouldBe` True+    it "CGS (2 x 2 dense)" $ +      normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true) <= eps `shouldBe` True+  describe "Math.Linear.Sparse : QR decomposition" $ do    +    it "QR (4 x 4 sparse)" $+      checkQr tm4 `shouldBe` True+    it "QR (3 x 3 dense)" $ +      checkQr tm2 `shouldBe` True+    +  -- let n = 10+  --     nsp = 3+  -- describe ("random sparse linear system of size " ++ show n ++ " and sparsity " ++ show (fromIntegral nsp/fromIntegral n)) $ it "<\\>" $ do+  --   aa <- randSpMat n nsp+  --   xtrue <- randSpVec n nsp+  --   b <- randSpVec n nsp    +  --   let b = aa #> xtrue+  --   printDenseSM aa+  --   normSq (aa <\> b ^-^ xtrue) <= eps `shouldBe` True+  -- --     normSq (_xBicgstab (bicgstab aa b x0 x0) ^-^ x) <= eps `shouldBe` True++++-- -- run N iterations ++-- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB+runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where+   x0 = mkSpVectorD nd $ replicate nd 0.9+   nd = dim r0+   r0 = b ^-^ (aa #> x0)    +   p0 = r0+   bicgsInit = BICGSTAB x0 r0 p0++-- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS+runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where+  x0 = mkSpVectorD nd $ replicate nd 0.1+  nd = dim r0+  r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution+  p0 = r0+  u0 = r0+  cgsInit = CGS x0 r0 p0 u0  +++{-++example 0 : 2x2 linear system++[1 2] [2] = [8]+[3 4] [3]   [18]+++[1 3] [2] = [11]+[2 4] [3]   [16]+++-}++aa0 :: SpMatrix Double+aa0 = SM (2,2) im where+  im = IM.fromList [(0, aa0r0), (1, aa0r1)]++aa0r0, aa0r1 :: IM.IntMap Double+aa0r0 = IM.fromList [(0,1),(1,2)]+aa0r1 = IM.fromList [(0,3),(1,4)]+++-- b0, x0 : r.h.s and initial solution resp.+b0, x0, x0true :: SpVector Double+b0 = mkSpVectorD 2 [8,18]+x0 = mkSpVectorD 2 [0.3,1.4]+++-- x0true : true solution+x0true = mkSpVectorD 2 [2,3]++++aa0tx0 = mkSpVectorD 2 [11,16]++++++++{- 4x4 system -}++aa1 :: SpMatrix Double+aa1 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]++x1, b1 :: SpVector Double+x1 = mkSpVectorD 4 [1,2,3,4]++b1 = mkSpVectorD 4 [30,56,60,101]++++{- 3x3 system -}+aa2 :: SpMatrix Double+aa2 = sparsifySM $ fromListDenseSM 3 [2, -1, 0, -1, 2, -1, 0, -1, 2]+x2, b2 :: SpVector Double+x2 = mkSpVectorD 3 [3,2,3]++b2 = mkSpVectorD 3 [4,-2,4]++++-- --++{-+example 1 : random linear system++-}++++-- dense+solveRandom n = do+  aa0 <- randMat n+  let aa = aa0 ^+^ eye n+  xtrue <- randVec n+  -- x0 <- randVec n+  let b = aa #> xtrue+      dx = aa <\> b ^-^ xtrue+  return $ normSq dx+  -- let xhatB = _xBicgstab (bicgstab aa b x0 x0)+  --     xhatC = _x (cgs aa b x0 x0)+  -- return (aa, x, x0, b, xhatB, xhatC)++-- sparse+solveSpRandom :: Int -> Int -> IO Double+solveSpRandom n nsp = do+  aa0 <- randSpMat n nsp+  let aa = aa0 ^+^ eye n+  xtrue <- randSpVec n nsp+  let b = (aa ^+^ eye n) #> xtrue+      dx = aa <\> b ^-^ xtrue+  return $ normSq dx++++-- `ndim` iterations++solveRandomN ndim nsp niter = do+  aa0 <- randSpMat ndim (nsp ^ 2)+  let aa = aa0 ^+^ eye ndim+  xtrue <- randSpVec ndim nsp+  let b = aa #> xtrue+      xhatB = head $ runNBiC niter aa b+      xhatC = head $ runNCGS niter aa b+  printDenseSM aa    +  return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))++--++{-+matMat++[1, 2] [5, 6] = [19, 22]+[3, 4] [7, 8]   [43, 50]+-}++m1 = fromListDenseSM 2 [1,3,2,4]+m2 = fromListDenseSM 2 [5, 7, 6, 8]     +m1m2 = fromListDenseSM 2 [19, 43, 22, 50]++-- transposeSM++m1t = fromListDenseSM 2 [1,2,3,4]+++--++{-+countSubdiagonalNZ+-}++m3 = fromListSM (3,3) [(0,2,3),(2,0,4),(1,1,3)] +++++{- mkSubDiagonal -}++testLaplacian1 :: Int -> SpMatrix Double+testLaplacian1 n = m where+  m :: SpMatrix Double+  m = mksd (-1) l1 ^+^+       mksd 0 l2 ^+^+       mksd 1 l3+    where+    mksd = mkSubDiagonal n+    l1 = replicate n (-1)+    l2 = replicate n 2+    l3 = l1+  -- x :: SpVector Double+  -- x = mkSpVectorD n (replicate n 2)+  -- b = m #> x++-- t3 n = normSq $ (aa <\> b) ^-^ xhat where+--   aa = testLaplacian1 n :: SpMatrix Double+--   xhat = mkSpVectorD n (concat $ replicate 20 [1,2,3,4,5]) :: SpVector Double+--   b = aa #> xhat++{- QR-}++checkQr :: SpMatrix Double -> Bool+checkQr a = c1 && c2 where+  (q, r) = qr a+  c1 = normFrobenius ((q #~# r) ^-^ a) <= eps+  c2 = isOrthogonalSM q+++aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double+++++{- eigenvalues -}+++aa3 = fromListDenseSM 3 [1,1,3,2,2,2,3,1,1] :: SpMatrix Double++b3 = mkSpVectorD 3 [1,1,1] :: SpVector Double+++++++-- test data++tm0, tm1, tm2, tm3, tm4 :: SpMatrix Double+tm0 = fromListSM (2,2) [(0,0,pi), (1,0,sqrt 2), (0,1, exp 1), (1,1,sqrt 5)]++tv0, tv1 :: SpVector Double+tv0 = mkSpVectorD 2 [5, 6]+++tv1 = SV 2 $ IM.singleton 0 1++-- wikipedia test matrix for Givens rotation++tm1 = sparsifySM $ fromListDenseSM 3 [6,5,0,5,1,4,0,4,3]++tm1g1 = givens tm1 1 0+tm1a2 = tm1g1 ## tm1++tm1g2 = givens tm1a2 2 1+tm1a3 = tm1g2 ## tm1a2++tm1q = transposeSM (tm1g2 ## tm1g1)+++-- wp test matrix for QR decomposition via Givens rotation++tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41]+++++tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9]++tm3g1 = fromListDenseSM 3 [1, 0,0, 0,c,-s, 0, s, c]+  where c= 0.4961+        s = 0.8682+++--++tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]