sparse-linear-algebra 0.2.0.9 → 0.2.1.0

raw patch · 10 files changed

+340/−173 lines, 10 filesPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

API changes (from Hackage documentation)

Files

• README.md +5/−3
• sparse-linear-algebra.cabal +1/−1
• src/Data/Sparse/Common.hs +21/−0
• src/Data/Sparse/IntMap2/IntMap2.hs +2/−3
• src/Data/Sparse/SpMatrix.hs +49/−20
• src/Data/Sparse/SpVector.hs +32/−4
• src/Numeric/Eps.hs +56/−17
• src/Numeric/LinearAlgebra/Class.hs +32/−9
• src/Numeric/LinearAlgebra/Sparse.hs +104/−102
• test/LibSpec.hs +38/−14

@@ -18,11 +18,13 @@      * Transpose-Free Quasi-Minimal Residual (TFQMR) -* Matrix decompositions+* Matrix factorization algorithms -    * QR factorization+    * QR -    * LU factorization+    * LU++    * Cholesky  * Eigenvalue algorithms 

@@ -1,5 +1,5 @@ name:                sparse-linear-algebra-version:             0.2.0.9+version:             0.2.1.0 synopsis:            Numerical computation in native Haskell  description:         Currently the library provides iterative linear solvers, matrix decompositions, eigenvalue computations and related utilities. Please see README.md for details homepage:            https://github.com/ocramz/sparse-linear-algebra

@@ -1,6 +1,16 @@+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (C) 2016 Marco Zocca+-- License     :  GPL-3 (see LICENSE)+-- Maintainer  :  zocca.marco gmail+-- Stability   :  provisional+-- Portability :  portable+--+----------------------------------------------------------------------------- module Data.Sparse.Common        ( module X,          insertRowWith, insertRow, insertColWith, insertCol,+         diagonalSM,          outerProdSV, (><), toSV, svToSM,           extractCol, extractRow,          extractVectorDenseWith, extractRowDense, extractColDense,@@ -58,6 +68,7 @@  -- * Outer vector product +-- | Outer product (all-with-all matrix) outerProdSV, (><) :: Num a => SpVector a -> SpVector a -> SpMatrix a outerProdSV v1 v2 = fromListSM (m, n) ixy where   m = dim v1@@ -65,6 +76,16 @@   ixy = [(i,j, x * y) | (i,x) <- toListSV v1 , (j, y) <- toListSV v2]  (><) = outerProdSV++++-- * Diagonal matrix++-- | Fill the diagonal of a SpMatrix with the components of a SpVector+diagonalSM :: SpVector a -> SpMatrix a+diagonalSM sv = ifoldSV iins (zeroSM n n) sv where+  n = dim sv+  iins i = insertSpMatrix i i   

@@ -35,10 +35,9 @@   --- * set-like brackets+-- set-like brackets  --- unionWithKeyIM2 f im1 im2 = undefined where   @@ -110,7 +109,7 @@   --- * filtering+-- * Filtering  -- |Map over outer IM and filter all inner IM's ifilterIM2 ::

@@ -1,4 +1,13 @@ {-# language FlexibleInstances, MultiParamTypeClasses, TypeFamilies #-}+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (C) 2016 Marco Zocca+-- License     :  GPL-3 (see LICENSE)+-- Maintainer  :  zocca.marco gmail+-- Stability   :  provisional+-- Portability :  portable+--+----------------------------------------------------------------------------- module Data.Sparse.SpMatrix where  import Data.Sparse.Utils@@ -15,6 +24,20 @@   +-- *++-- instance IxContainer SpMatrix a where+--   type Ix SpMatrix = (Int, Int)+--   type IxSz SpMatrix = (Int, Int)+--   ixcLookup m (i,j) = lookupSM m i j+--   ixcIfilter g = ifilterSM g' where g' i j = g (i, j)+--   ixcInsert (i, j) = insertSpMatrix i j+--   -- ixcFromList = fromListSM++-- instance SMatrix SpMatrix a where+++ -- * Sparse Matrix  data SpMatrix a = SM {smDim :: (Rows, Cols),@@ -88,7 +111,7 @@ permutationSM n iis = permutPairsSM n (zip [0 .. n-1] iis)  -- | Permutation matrix from a (possibly incomplete) list of row pair swaps--- e.g. `permutPairs 5 [(2,4)]` swaps rows (2, 4) :+-- e.g. `permutPairs 5 [(2,4)]` swaps rows 2 and 4 : -- -- [1,0,0,0,0] -- [0,1,0,0,0]@@ -159,7 +182,7 @@  -- ** toList --- |Populate list with SpMatrix contents and populate missing entries with 0+-- | Populate list with SpMatrix contents and 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]]@@ -263,6 +286,7 @@   -- *** Extract i'th row+-- | Extract whole row extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a extractRowSM sm i = extractSubmatrix sm (i, i) (0, ncols sm - 1) @@ -279,7 +303,7 @@   -- *** Extract j'th column--- | Extract all column+-- | Extract whole column extractColSM :: SpMatrix a -> IxCol -> SpMatrix a extractColSM sm j = extractSubmatrix sm (0, nrows sm - 1) (j, j) @@ -323,8 +347,8 @@   ff irow row = IM.size row == 1 &&                 IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1 --- |is the matrix orthogonal? i.e. Q^t ## Q == I-isOrthogonalSM :: SpMatrix Double -> 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 @@ -497,6 +521,11 @@  -- ** Misc. SpMatrix operations +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 foldlSM :: (a -> b -> b) -> b -> SpMatrix a -> b foldlSM f n (SM _ m)= foldlIM2 f n m@@ -536,11 +565,12 @@   -- ** Sparsify : remove almost-0 elements (|x| < eps)-sparsifyIM2 :: IM.IntMap (IM.IntMap Double) -> IM.IntMap (IM.IntMap Double)-sparsifyIM2 = ifilterIM2 (\_ _ x -> abs x >= eps)+sparsifyIM2 ::+  Epsilon a => IM.IntMap (IM.IntMap a) -> IM.IntMap (IM.IntMap a)+sparsifyIM2 = ifilterIM2 (\_ _ x -> isNz x)  -- | Sparsify an SpMatrix-sparsifySM :: SpMatrix Double -> SpMatrix Double+sparsifySM :: Epsilon a => SpMatrix a -> SpMatrix a sparsifySM (SM d im) = SM d $ sparsifyIM2 im  @@ -548,7 +578,7 @@  -- ** Value rounding -- | Round almost-0 and almost-1 to 0 and 1 respectively-roundZeroOneSM :: SpMatrix Double -> SpMatrix Double+roundZeroOneSM :: Epsilon a => SpMatrix a -> SpMatrix a roundZeroOneSM (SM d im) = sparsifySM $ SM d $ mapIM2 roundZeroOne im    @@ -560,14 +590,14 @@   -- ** Matrix row swap--- | swap two rows of a SpMatrix (bounds not checked)+-- | 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 --- | swap two rows of a SpMatrix (bounds checked)  +-- | Swap two rows of a SpMatrix (bounds checked)   swapRowsSafe :: IxRow -> IxRow -> SpMatrix a -> SpMatrix a swapRowsSafe i1 i2 m   | inBounds02 (nro, nro) (i1, i2) = swapRows i1 i2 m@@ -581,11 +611,10 @@   -- ** Matrix transpose--- | transposeSM, (#^) : Matrix transpose-transposeSM, (#^) :: SpMatrix a -> SpMatrix a+-- | transposeSM : Matrix transpose+transposeSM :: SpMatrix a -> SpMatrix a transposeSM (SM (m, n) im) = SM (n, m) (transposeIM2 im) -(#^) = transposeSM   @@ -603,7 +632,7 @@   -- ** Frobenius norm-normFrobenius :: SpMatrix Double -> Double+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 @@ -655,7 +684,7 @@  -- ** Matrix-matrix product, sparsified -- | Removes all elements `x` for which `| x | <= eps`)-matMatSparsified, (#~#)  :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+matMatSparsified, (#~#) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a matMatSparsified m1 m2 = sparsifySM $ matMat m1 m2  (#~#) = matMatSparsified@@ -666,12 +695,12 @@ -- *** Sparsified matrix products of two matrices  -- | A^T B-(#^#) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+(#^#) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a a #^# b = transposeSM a #~# b   -- | A B^T-(##^) :: SpMatrix Double -> SpMatrix Double -> SpMatrix Double+(##^) :: Epsilon a => SpMatrix a -> SpMatrix a -> SpMatrix a a ##^ b = a #~# transposeSM b  @@ -686,8 +715,8 @@   --- *** Matrix contraction--- | Contract two matrices A and B up to an index `n`, i.e. summing over repeated indices: +-- ** 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 a b i j n

@@ -1,12 +1,23 @@ {-# language TypeFamilies, MultiParamTypeClasses, FlexibleInstances #-}+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (C) 2016 Marco Zocca+-- License     :  GPL-3 (see LICENSE)+-- Maintainer  :  zocca.marco gmail+-- Stability   :  provisional+-- Portability :  portable+--+----------------------------------------------------------------------------- module Data.Sparse.SpVector where  import Data.Sparse.Utils import Data.Sparse.Types+import Data.Sparse.IntMap2.IntMap2 +import Numeric.Eps import Numeric.LinearAlgebra.Class-import Data.Sparse.IntMap2.IntMap2 + import Data.Maybe  import qualified Data.IntMap as IM@@ -94,11 +105,11 @@   -- | create a sparse vector from an association list while discarding all zero entries-mkSpVector :: (Num a, Eq a) => Int -> IM.IntMap a -> SpVector a-mkSpVector d im = SV d $ IM.filterWithKey (\k v -> v /= 0 && inBounds0 d k) im+mkSpVector :: Epsilon a => Int -> IM.IntMap a -> SpVector a+mkSpVector d im = SV d $ IM.filterWithKey (\k v -> isNz v && inBounds0 d k) im  -- | ", from logically dense array (consecutive indices)-mkSpVectorD :: (Num a, Eq a) => Int -> [a] -> SpVector a+mkSpVectorD :: Epsilon a => Int -> [a] -> SpVector a mkSpVectorD d ll = mkSpVector d (IM.fromList $ denseIxArray (take d ll))  -- ", don't filter zero elements@@ -110,6 +121,18 @@ fromListDenseSV d ll = SV d (IM.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+spVectorDenseIx :: Epsilon a => (Int -> a) -> UB -> [Int] -> SpVector a+spVectorDenseIx f n ix =+  fromListSV n $ filter q $ zip ix $ map f ix where+    q (i, v) = inBounds0 n i && isNz v++-- | ", using just the integer bounds of the interval+spVectorDenseLoHi :: Epsilon a => (Int -> a) -> UB -> Int -> Int -> SpVector a+spVectorDenseLoHi f n lo hi = spVectorDenseIx f n [lo .. hi]    +++ -- | 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)@@ -152,6 +175,11 @@ toDenseListSV (SV d 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   

@@ -1,36 +1,75 @@-module Numeric.Eps where+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (C) 2016 Marco Zocca, 2012-2015 Edward Kmett+-- License     :  GPL-3 (see LICENSE)+-- Maintainer  :  zocca.marco gmail+-- Stability   :  provisional+-- Portability :  portable+--+-- Testing for values "near" zero+-----------------------------------------------------------------------------+module Numeric.Eps+  ( Epsilon(..), nearZero, isNz, roundZero, roundOne, roundZeroOne+  ) where+import Foreign.C.Types (CFloat, CDouble) --- * Numerical tolerance for "near-0" tests--- | eps = 1e-8 -eps :: Double-eps = 1e-8+-- | Provides a test to see if a quantity is near zero.+--+-- >>> nearZero (1e-11 :: Double)+-- False+--+-- >>> nearZero (1e-17 :: Double)+-- True+--+-- >>> nearZero (1e-5 :: Float)+-- False+--+-- >>> nearZero (1e-7 :: Float)+-- True+class Num a => Epsilon a where+  -- | Determine if a quantity is near zero.+  nearZero :: a -> Bool +-- | @'abs' a '<=' 1e-6@+instance Epsilon Float where+  nearZero a = abs a <= 1e-6 +-- | @'abs' a '<=' 1e-12@+instance Epsilon Double where+  nearZero a = abs a <= 1e-12 +-- | @'abs' a '<=' 1e-6@+instance Epsilon CFloat where+  nearZero a = abs a <= 1e-6 +-- | @'abs' a '<=' 1e-12@+instance Epsilon CDouble where+  nearZero a = abs a <= 1e-12++++ -- * Rounding operations --- | Rounding rule-almostZero, almostOne :: Double -> Bool-almostZero x = abs x <= eps-almostOne x = x >= (1-eps) && x < (1+eps) -isNz :: Double -> Bool-isNz = not . almostZero+-- | Rounding rule+almostZero, almostOne, isNz :: Epsilon a => a -> Bool+almostZero = nearZero+almostOne x = nearZero (1 - x)+isNz x = not (almostZero x)  withDefault :: (t -> Bool) -> t -> t -> t withDefault q d x | q x = d                   | otherwise = x -roundZero, roundOne :: Double -> Double-roundZero = withDefault almostZero 0-roundOne = withDefault almostOne 1+roundZero, roundOne, roundZeroOne :: Epsilon a => a -> a+roundZero = withDefault almostZero (fromIntegral 0)+roundOne = withDefault almostOne (fromIntegral 1)  with2Defaults :: (t -> Bool) -> (t -> Bool) -> t -> t -> t -> t with2Defaults q1 q2 d1 d2 x | q1 x = d1                             | q2 x = d2                             | otherwise = x --- | Round to respectively 0 or 1 within some predefined numerical precision eps-roundZeroOne :: Double -> Double-roundZeroOne = with2Defaults almostZero almostOne 0 1+-- | Round to respectively 0 or 1+roundZeroOne = with2Defaults almostZero almostOne (fromIntegral 0) (fromIntegral 1)

@@ -10,8 +10,12 @@   (^+^) :: Num a => f a -> f a -> f a  +  one :: Num 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@@ -185,15 +189,14 @@  -- * IxContainer : indexed container types -class IxContainer (c :: * -> *) a where-  type Ix c :: *-  ixcLookup :: Ix c -> c a -> Maybe a-  ixcLookupDefault :: a -> Ix c -> c a -> a-  ixcFilter :: (a -> Bool) -> c a -> c a-  ixcIfilter :: (Ix c -> a -> Bool) -> c a -> c a-  ixcInsert :: Ix c -> a -> c a -> c a-  ixcFromList :: [(Ix c, a)] -> c a-  ixcToList :: c a -> [(Ix c, a)]+-- 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) @@ -208,3 +211,23 @@ -- 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++-- class (IxContainer c a, Sparse c a, Additive c) => SMatrix c a where

@@ -4,6 +4,7 @@        (          -- * Matrix factorizations          qr, lu,+         chol,          -- * Incomplete LU          ilu0,          -- * Condition number@@ -21,6 +22,8 @@          _xCgne, _xTfq, _xBicgstab, _x, _xBcg,          cgsStep, bicgstabStep,          CGNE, TFQMR, BICGSTAB, CGS, BCG,+         -- * Preconditioners+         ilu0, mSsor,          -- * Matrix partitioning          diagPartitions,          -- * Random arrays@@ -69,8 +72,8 @@    -- * Sparsify : remove almost-0 elements (|x| < eps) -- | Sparsify an SpVector-sparsifySV :: SpVector Double -> SpVector Double-sparsifySV (SV d im) = SV d $ IM.filter (\x -> abs x >= eps) im+sparsifySV :: Epsilon a => SpVector a -> SpVector a+sparsifySV (SV d im) = SV d $ IM.filter isNz im   @@ -79,7 +82,7 @@ -- * Matrix condition number  -- |uses the R matrix from the QR factorization-conditionNumberSM :: SpMatrix Double -> Double+conditionNumberSM :: (Epsilon a, RealFloat a) => SpMatrix a -> a conditionNumberSM m | isInfinite kappa = error "Infinite condition number : rank-deficient system"                     | otherwise = kappa where   kappa = lmax / lmin@@ -104,8 +107,8 @@ {-| 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+hhRefl :: Num a => SpVector a -> SpMatrix a+hhRefl = hhMat (fromInteger 2)   @@ -151,7 +154,7 @@ non-zero but A has zeros in row k for all columns less than j. -} -givens :: SpMatrix Double -> IxRow -> IxCol -> SpMatrix Double+givens :: (Floating a, Epsilon a, Ord a) => SpMatrix a -> IxRow -> IxCol -> SpMatrix a givens mm i j    | isValidIxSM mm (i,j) && isSquareSM mm =        sparsifySM $ fromListSM' [(i,i,c),(j,j,c),(j,i,-s),(i,j,s)] (eye (nrows mm))@@ -182,15 +185,15 @@ -- * QR decomposition  --- | Applies Givens rotation iteratively to zero out sub-diagonal elements-qr :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double)+-- | Given a matrix A, returns a pair of matrices (Q, R) such that Q R = A, Q is orthogonal and R is upper triangular. Applies Givens rotation iteratively to zero out sub-diagonal elements+qr :: (Epsilon a, Floating a, Real a) => SpMatrix a -> (SpMatrix a, SpMatrix a) 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 ]-gmats :: SpMatrix Double -> [SpMatrix Double]+gmats :: (Epsilon a, Real a, Floating a) => SpMatrix a -> [SpMatrix a] gmats mm = gm mm (subdiagIndicesSM mm) where  gm m ((i,j):is) = let g = givens m i j                    in g : gm (g #~# m) is@@ -218,13 +221,13 @@ -- ** QR algorithm  -- | `eigsQR n mm` performs `n` iterations of the QR algorithm on matrix `mm`, and returns a SpVector containing all eigenvalues-eigsQR :: Int -> SpMatrix Double -> SpVector Double+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 norm2 (dm1 ^-^ dm2) <= eps+                 in nearZero $ norm2 (dm1 ^-^ dm2)   @@ -234,13 +237,13 @@ -- ** 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 :: 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, _)] = norm2 (b2 ^-^ b1) <= eps +    f [(b1, _), (b2, _)] = nearZero $ norm2 (b2 ^-^ b1)   rayleighStep aa (b, mu) = (b', mu') where       ii = eye (nrows aa)       nom = (aa ^-^ (mu `matScale` ii)) <\> b@@ -253,13 +256,13 @@ -- * Householder vector   -- (Golub & Van Loan, Alg. 5.1.1, function `house`)-hhV :: SpVector Double -> (SpVector Double, Double)+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) | sigma <= eps = (vtemp, 0)+  (v, beta) | nearZero sigma = (vtemp, 0)             | otherwise = let mu = sqrt (headSV x**2 + sigma)                               xh = headSV x                               vh | xh <= 1 = xh - mu@@ -298,109 +301,92 @@   +-- * 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 -} --- | LU factors-lu :: SpMatrix Double -> (SpMatrix Double, SpMatrix Double)-lu aa = (lfin, ufin) where-  (ixf,lf,uf) = execState (modifyUntil q (luUpd aa)) (luInit aa)-  lfin = lf-  ufin = uUpd aa (ixf, lf, uf)-  q (i, _, _) = i == (nrows aa - 1) --- | First iteration of LU-luInit ::-  (Num t, Fractional a) => SpMatrix a -> (t, SpMatrix a, SpMatrix a)-luInit aa = (1, l0, u0) where-  n = nrows aa-  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 --- | LU update step-luUpd :: SpMatrix Double-     -> (Int, SpMatrix Double, SpMatrix Double)-     -> (Int, SpMatrix Double, SpMatrix Double)-luUpd aa (i, l, u) = (i', l', u') where-  n = nrows aa  -  u' = uUpdSparse aa (i, l, u)  -- update U-  l' = lUpdSparse aa (i, l, u') -- update L-  i' = i + 1     -- increment i  -uUpd' ::-  Num a =>-  ([(Int, a)] -> [(Int, a)]) ->-  SpMatrix a ->-  (Rows, SpMatrix a, SpMatrix a) ->-  SpMatrix a-uUpd' ff amat (ix, lmat, umat) = insertRow umat uv ix where-  n = nrows amat-  colsix = [ix .. n - 1]-  us = ff $ zip colsix $ map (solveForUij amat lmat umat ix) colsix-  uv = fromListSV n us -uUpd :: Num a => SpMatrix a -> (Rows, SpMatrix a, SpMatrix a) -> SpMatrix a-uUpd = uUpd' id --- update U while sparsifying-uUpdSparse ::-  SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double-uUpdSparse = uUpd' (filter (isNz . snd))    --- solve for element Uij-solveForUij ::-  Num a => SpMatrix a -> SpMatrix a -> SpMatrix a -> IxRow -> IxCol -> a-solveForUij amat lmat umat i j = a - p where-  a = amat @@! (i, j)-  p = contractSub lmat umat i j (i - 1) +-- * 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 -} --- solve for element Lij-solveForLij ::-  SpMatrix Double -> SpMatrix Double -> SpMatrix Double -> IxRow -> IxCol -> Double-solveForLij amat lmat umat 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 = amat @@! (i, j)-   ujj = umat @@! (j , j)   -- NB this must be /= 0-   p = contractSub lmat umat i j (i - 1)+-- | Given a matrix A, returns a pair of matrices (L, U) 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, Real b) => (Int -> b) -> [Int] -> [(Int, b)]+onRangeSparse f ixs = filter (isNz . snd) $ zip ixs $ map f ixs -lUpd' :: ([(Rows, Double)] -> [(Int, Double)])-     -> SpMatrix Double-     -> (Rows, SpMatrix Double, SpMatrix Double)-     -> SpMatrix Double-lUpd' ff amat (ix, lmat, umat) = insertCol lmat lv ix where-  n = nrows amat-  rowsix = [ix + 1 .. n - 1]-  ls = ff $ zip rowsix $ map (\i -> solveForLij amat lmat umat i ix) rowsix-  lv = fromListSV n ls -lUpd :: SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double-lUpd = lUpd' id -lUpdSparse ::-  SpMatrix Double -> (Rows, SpMatrix Double, SpMatrix Double) -> SpMatrix Double-lUpdSparse = lUpd' (filter (isNz . snd))   @@ -442,7 +428,7 @@ -- | 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@@ -471,13 +457,13 @@  -- ** SSOR --- | `mSsor aa omega` : if `omega = 1` it returns the diagonal of `aa`, -mSsor :: Fractional a => SpMatrix a -> a -> SpMatrix a-mSsor aa omega = l ## r where+-- | `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 = d ^-^ scale omega e-  r = eye n ^-^ scale omega (reciprocal d ## f)+  l = (eye n ^-^ scale omega e) ## reciprocal d+  r = d ^-^ scale omega f    @@ -486,11 +472,27 @@   +-- Linear solver, LU-based+++++++++ -- * Iterative linear solvers  +-- ** GMRES +-- *** Left-preconditioning ++++ -- ** CGNE  cgneStep :: SpMatrix Double -> CGNE -> CGNE@@ -840,7 +842,7 @@ -- | convergence check (FIXME) normDiffConverged :: (Foldable t, Functor t) =>      (a -> SpVector Double) -> t a -> Bool-normDiffConverged fp xx = normSq (foldrMap fp (^-^) (zeroSV 0) xx) <= eps+normDiffConverged fp xx = nearZero $ normSq (foldrMap fp (^-^) (zeroSV 0) xx)     

@@ -1,4 +1,13 @@ {-# language ScopedTypeVariables #-}+-----------------------------------------------------------------------------+-- |+-- Copyright   :  (C) 2016 Marco Zocca+-- License     :  GPL-3 (see LICENSE)+-- Maintainer  :  zocca.marco gmail+-- Stability   :  provisional+-- Portability :  portable+--+----------------------------------------------------------------------------- module LibSpec where  import Numeric.LinearAlgebra.Sparse@@ -40,9 +49,9 @@     it "transposeSM : sparse matrix transpose" $       transposeSM m1 `shouldBe` m1t     it "matVec : matrix-vector product" $-      normSq ((aa0 #> x0true) ^-^ b0 ) <= eps `shouldBe` True+      nearZero ( normSq ((aa0 #> x0true) ^-^ b0 )) `shouldBe` True     it "vecMat : vector-matrix product" $-      normSq ((x0true <# aa0) ^-^ aa0tx0 ) <= eps `shouldBe` True  +      nearZero ( normSq ((x0true <# aa0) ^-^ aa0tx0 ))`shouldBe` True       it "matMat : matrix-matrix product" $       (m1 `matMat` m2) `shouldBe` m1m2     it "eye : identity matrix" $@@ -58,22 +67,22 @@     it "countSubdiagonalNZ : # of nonzero elements below the diagonal" $       countSubdiagonalNZSM m3 `shouldBe` 1     it "permutPairsSM : permutation matrices are orthogonal" $ do-      let pm0 = permutPairsSM 3 [(0,2), (1,2)]+      let pm0 = permutPairsSM 3 [(0,2), (1,2)] :: SpMatrix Double       pm0 ##^ pm0 `shouldBe` eye 3       pm0 #^# pm0 `shouldBe` eye 3              it "modifyInspectN : early termination by iteration count" $-      execState (modifyInspectN 2 ((< eps) . diffSqL) (/2)) 1 `shouldBe` 1/8+      execState (modifyInspectN 2 (nearZero . diffSqL) (/2)) (1 :: Double) `shouldBe` 1/8     it "modifyInspectN : termination by value convergence" $-      execState (modifyInspectN (2^16) ((< eps) . head) (/2)) 1 < eps `shouldBe` True +      nearZero (execState (modifyInspectN (2^16) (nearZero . head) (/2)) (1 :: Double)) `shouldBe` True    describe "Numeric.LinearAlgebra.Sparse : Linear solvers" $ do     -- it "TFQMR (2 x 2 dense)" $     --   normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True     it "BCG (2 x 2 dense)" $-      normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True+      nearZero (normSq (_xBcg (bcg aa0 b0 x0) ^-^ x0true)) `shouldBe` True     it "BiCGSTAB (2 x 2 dense)" $ -      normSq (aa0 <\> b0 ^-^ x0true) <= eps `shouldBe` True+      nearZero (normSq (aa0 <\> b0 ^-^ x0true)) `shouldBe` True     it "CGS (2 x 2 dense)" $ -      normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true) <= eps `shouldBe` True+      nearZero (normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true)) `shouldBe` True   describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do         it "QR (4 x 4 sparse)" $       checkQr tm4 `shouldBe` True@@ -84,6 +93,9 @@       checkLu tm6 `shouldBe` True     it "LU (10 x 10 sparse)" $       checkLu tm7 `shouldBe` True+  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ do+    it "chol (5 x 5 sparse)" $+      checkChol tm7 `shouldBe` True   {-@@ -252,10 +264,11 @@  {- QR-} -checkQr :: SpMatrix Double -> Bool++checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool checkQr a = c1 && c2 where   (q, r) = qr a-  c1 = normFrobenius ((q #~# r) ^-^ a) <= eps+  c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)   c2 = isOrthogonalSM q  @@ -266,13 +279,24 @@  {- LU -} -checkLu :: SpMatrix Double -> Bool+checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool checkLu a = lup == a where   (l, u) = lu a   lup = l #~# u   +{- Cholesky -}++checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool+checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where+  l = chol a++++++ {- eigenvalues -}  @@ -345,9 +369,9 @@ tm7 :: SpMatrix Double tm7 = a ^+^ b ^+^ c where   n = 5-  a = mkSubDiagonal n 1 $ replicate n 1-  b = mkSubDiagonal n 0 $ replicate n (-2)-  c = mkSubDiagonal n (-1) $ replicate n 1+  a = mkSubDiagonal n 1 $ replicate n (-1)+  b = mkSubDiagonal n 0 $ replicate n 2+  c = mkSubDiagonal n (-1) $ replicate n (-1)  -- -- run N iterations