sparse-linear-algebra 0.2.1.0 → 0.2.1.1

raw patch · 10 files changed

+381/−111 lines, 10 filesdep +vectorPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: vector

API changes (from Hackage documentation)

Files

• README.md +38/−16
• sparse-linear-algebra.cabal +2/−1
• src/Data/Sparse/Common.hs +57/−7
• src/Data/Sparse/IntMap2/IntMap2.hs +10/−2
• src/Data/Sparse/SpMatrix.hs +22/−18
• src/Data/Sparse/SpVector.hs +44/−0
• src/Data/Sparse/Utils.hs +22/−1
• src/Numeric/LinearAlgebra/Class.hs +9/−2
• src/Numeric/LinearAlgebra/Sparse.hs +105/−23
• test/LibSpec.hs +72/−41

@@ -6,32 +6,38 @@  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. -Algorithms :+Contents :  * Iterative linear solvers -    * BiConjugate Gradient (BCG)+    * BiConjugate Gradient (`bcg`) -    * Conjugate Gradient Squared (CGS)+    * Conjugate Gradient Squared (`cgs`) -    * BiConjugate Gradient Stabilized (BiCGSTAB) (non-Hermitian systems)+    * BiConjugate Gradient Stabilized (`bicgstab`) (non-Hermitian systems) -    * Transpose-Free Quasi-Minimal Residual (TFQMR)+    * Transpose-Free Quasi-Minimal Residual (`tfqmr`) +* Direct linear solvers++    * LU-based (`luSolve`)+ * Matrix factorization algorithms -    * QR+    * QR (`qr`) -    * LU+    * LU (`lu`) -    * Cholesky+    * Cholesky (`chol`)  * Eigenvalue algorithms -    * QR algorithm+    * Arnoldi iteration (`arnoldi`) -    * Rayleigh quotient iteration+    * QR (`eigsQR`) +    * Rayleigh quotient iteration (`eigRayleigh`)+ * Utilities : Vector and matrix norms, matrix condition number, Givens rotation, Householder reflection  * Predicates : Matrix orthogonality test (A^T A ~= I)@@ -43,9 +49,9 @@  The module `Numeric.LinearAlgebra.Sparse` contains the user interface. -### Creation and pretty-printing+### Creation of sparse data -To create a sparse matrix from an array of its entries we use `fromListSM` :+The `fromListSM` function creates a sparse matrix from an array of its entries we use :      fromListSM :: Foldable t => (Int, Int) -> t (IxRow, IxCol, a) -> SpMatrix a @@ -53,8 +59,14 @@      > amat = fromListSM (3,3) [(0,0,2),(1,0,4),(1,1,3),(1,2,2),(2,2,5)] -And similarly for sparse vectors : `fromListSV :: Int -> [(Int, a)] -> SpVector a`.+and similarly +    fromListSV :: Int -> [(Int, a)] -> SpVector a++can be used to create sparse vectors.++### Displaying sparse data+ Both sparse vectors and matrices can be pretty-printed using `prd`:      > prd amat@@ -64,7 +76,7 @@     [4,3,2]     [0,0,5] -The zeros are just added at pretty printing time; sparse vectors and matrices should only contain non-zero entries.+The zeros are just added at printing time; sparse vectors and matrices should only contain non-zero entries.  ### Matrix operations @@ -79,7 +91,7 @@     [0.0,0.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`) is fixed at 10^-8.+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     ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )@@ -91,7 +103,7 @@  ### Linear systems -Linear systems can be solved with either `linSolve` (which also requires choosing a method) or with `<\>` (which uses BiCGSTAB as default) :+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 BiCGSTAB as default) :      > b = fromListSV 3 [(0,3),(1,2),(2,5)]     > x = amat <\> b@@ -106,6 +118,16 @@     ( 3 elements ) ,  3 NZ ( sparsity 1.0 )      [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'++    ( 3 elements ) ,  3 NZ ( sparsity 1.0 )++    [1.5,-2.0,1.0]+   

@@ -1,5 +1,5 @@ name:                sparse-linear-algebra-version:             0.2.1.0+version:             0.2.1.1 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@@ -34,6 +34,7 @@                      , primitive >= 0.6.1.0                      , mtl >= 2.2.1                      , mwc-random+                     , vector   executable sparse-linear-algebra

@@ -11,18 +11,20 @@        ( module X,          insertRowWith, insertRow, insertColWith, insertCol,          diagonalSM,-         outerProdSV, (><), toSV, svToSM, +         outerProdSV, (><), toSV, svToSM,+         lookupRowSM,           extractCol, extractRow,          extractVectorDenseWith, extractRowDense, extractColDense,          extractDiagDense,          extractSubRow, extractSubCol,          extractSubRow_RK, extractSubCol_RK,          matVec, (#>), vecMat, (<#),+         fromCols,          prd) where  import Data.Sparse.Utils as X import Data.Sparse.Types as X-import Data.Sparse.IntMap2.IntMap2 as X+import Data.Sparse.IntMap2.IntMap2 -- as X import Data.Sparse.SpMatrix as X import Data.Sparse.SpVector as X @@ -31,8 +33,20 @@  import qualified Data.IntMap as IM +import Data.Maybe (fromMaybe, maybe)+import qualified Data.Vector as V +-- withBoundsSM m ij e f+--   | isValidIxSM m ij = f m ij+--   | otherwise = error e +-- | 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+ -- * Insert row/column vector in matrix  -- | Insert row , using the provided row index transformation function@@ -107,12 +121,20 @@     | 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) ++ -- * Extract a SpVector from an SpMatrix -- ** Sparse extract  -- |Extract ith row extractRow :: SpMatrix a -> IxRow -> SpVector a-extractRow m i = toSV $ extractRowSM m i+extractRow m i+  | inBounds0 (nrows m) i = fromMaybe (zeroSV (ncols m)) (lookupRowSM m i)+  | otherwise = error $ unwords ["extractRow : index",show i,"out of bounds"]  -- |Extract jth column extractCol :: SpMatrix a -> IxCol -> SpVector a@@ -143,14 +165,32 @@   --- | extract row interval-extractSubRow :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a-extractSubRow m i (j1, j2)  = toSV $ extractSubRowSM m i (j1, j2)+-- | extract row interval (all entries between columns j1 and j2, INCLUDED, are returned)+-- extractSubRow :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a+-- extractSubRow m i (j1, j2) = case lookupRowSM m i of+--   Nothing -> zeroSV (ncols m)+--   Just rv -> ifilterSV (\j _ -> j >= j1 && j <= j2) rv +-- |", returning in Maybe+-- extractSubRow :: SpMatrix a -> IxRow -> (Int, Int) -> Maybe (SpVector a)+-- extractSubRow m i (j1, j2) =+--   resizeSV (j2 - j1) . ifilterSV (\j _ -> j >= j1 && j <= j2) <$> lookupRowSM m i++-- | Extract an interval of SpVector components, changing accordingly the resulting SpVector size. Keys are _not_ rebalanced, i.e. components are still labeled according with respect to the source matrix.+extractSubRow :: SpMatrix a -> IxRow -> (Int, Int) -> SpVector a+extractSubRow m i (j1, j2) = fromMaybe (zeroSV deltaj) vfilt where+  deltaj = j2 - j1 + 1+  vfilt = resizeSV deltaj .+          ifilterSV (\j _ -> j >= j1 && j <= j2) <$> lookupRowSM m i+ -- | extract row interval, rebalance keys by subtracting lowest one extractSubRow_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpVector a-extractSubRow_RK m i (j1, j2)  = toSV $ extractSubRowSM_RK m i (j1, j2)+extractSubRow_RK m i (j1, j2) = mapKeysSV (subtract j1) $ extractSubRow m i (j1, j2) +  -- toSV $ extractSubRowSM_RK m i (j1, j2)+++ -- | extract column interval extractSubCol :: SpMatrix a -> IxCol -> (IxRow, IxRow) -> SpVector a extractSubCol m j (i1, i2)  = toSV $ extractSubColSM m j (i1, i2)@@ -190,6 +230,16 @@ (<#) = vecMat    +++++-- | 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+  n = V.length qv+  m = dim $ V.head qv+  ins mm i c = insertCol mm c i   

@@ -3,8 +3,9 @@ import Numeric.LinearAlgebra.Class  import qualified Data.IntMap.Strict as IM-+import Data.Sparse.Types +import Data.Maybe   @@ -60,7 +61,14 @@ lookupIM2 i j imm = IM.lookup i imm >>= IM.lookup j {-# inline lookupIM2 #-}   --- |Ppopulate an IM2 from a list of (row index, column index, value)  +-- | Lookup with default 0+lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a+lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j)+++++-- |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)

@@ -75,7 +75,12 @@ instance Sparse SpMatrix a where   spy = spySM -+instance Num a => SpContainer SpMatrix a where+  type ScIx SpMatrix = (Rows, Cols)+  scInsert (i,j) = insertSpMatrix i j+  scLookup m (i, j) = lookupSM m i j+  m @@ d | isValidIxSM m d = m @@! d+         | otherwise = error $ "@@ : incompatible indices : matrix size is " ++ show (dim m) ++ ", but user looked up " ++ show d   @@ -200,28 +205,29 @@  -- ** Lookup +++ lookupSM :: SpMatrix a -> IxRow -> IxCol -> Maybe a lookupSM (SM _ im) i j = IM.lookup i im >>= IM.lookup j  -- | 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 -> (IxRow, IxCol) -> a+lookupWD_SM, (@@!):: Num a => SpMatrix a -> (IxRow, IxCol) -> a lookupWD_SM sm (i,j) =   fromMaybe 0 (lookupSM sm i j) -lookupWD_IM :: Num a => IM.IntMap (IM.IntMap a) -> (IxRow, IxCol) -> a-lookupWD_IM im (i,j) = fromMaybe 0 (IM.lookup i im >>= IM.lookup j) + -- | Zero-default lookup, infix form (no bound checking) (@@!) = lookupWD_SM --- | Zero-default lookup, infix form ("safe" : throws exception if lookup is outside matrix bounds)-m @@ d | isValidIxSM m d = m @@! d-       | otherwise = error $ "@@ : incompatible indices : matrix size is " ++ show (dim m) ++ ", but user looked up " ++ show d    ++ -- FIXME : to throw an exception or just ignore the out-of-bound access ?  @@ -287,21 +293,11 @@  -- *** Extract i'th row -- | Extract whole row-extractRowSM :: SpMatrix a -> IxRow -> SpMatrix a-extractRowSM sm i = extractSubmatrix sm (i, i) (0, ncols sm - 1)+-- -- moved to Data.Sparse.Common  --- | Extract column within a row range-extractSubRowSM :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a-extractSubRowSM sm i (j1, j2) = extractSubmatrix sm (i, i) (j1, j2) --- | Extract column within a row range, rebalance keys-extractSubRowSM_RK :: SpMatrix a -> IxRow -> (IxCol, IxCol) -> SpMatrix a-extractSubRowSM_RK sm i =-  extractSubmatrixRebalanceKeys sm (i, i)  -- -- *** Extract j'th column -- | Extract whole column extractColSM :: SpMatrix a -> IxCol -> SpMatrix a@@ -346,6 +342,14 @@   d = IM.filterWithKey ff (immSM m)   ff irow row = IM.size row == 1 &&                 IM.size (IM.filterWithKey (\j _ -> j == irow) row) == 1++-- | Is the matrix lower/upper triangular?+isLowerTriSM, isUpperTriSM :: Eq a => SpMatrix a -> Bool+isLowerTriSM m = m == lm where+  lm = ifilterSM (\i j _ -> i >= j) m++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

@@ -22,6 +22,7 @@  import qualified Data.IntMap as IM import qualified Data.Foldable as F+import qualified Data.Vector as V  -- * Sparse Vector @@ -74,6 +75,16 @@   spy = spySV  ++instance Num a => SpContainer SpVector a where+  type ScIx SpVector = Int+  scInsert = insertSpVector+  scLookup v i = lookupSV i v+  v @@ i = lookupDenseSV i v++++ instance Hilbert SpVector where   a `dot` b | dim a == dim b = dot (dat a) (dat b)             | otherwise =@@ -133,6 +144,10 @@   ++  ++ -- | 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 +167,28 @@   +++-- *** Vector-related++-- | Populate a SpVector with the contents of a Vector. +fromVector :: V.Vector a -> SpVector a+fromVector qv = V.ifoldl' ins (zeroSV n) qv where+  n = V.length qv+  ins vv i x = insertSpVector i x vv++-- | Populate a Vector with the entries of a SpVector, discarding the indices (NB: loses sparsity information).+toVector :: SpVector a -> V.Vector a+toVector = V.fromList . snd . unzip . toListSV++-- | -- | Populate a Vector with the entries of a SpVector, replacing the missing entries with 0+toVectorDense :: Num a => SpVector a -> V.Vector a+toVectorDense = V.fromList . toDenseListSV+++++ -- ** Element insertion  -- |insert element `x` at index `i` in a preexisting SpVector@@ -220,6 +257,13 @@ -- | Head element headSV :: Num a => SpVector a -> a headSV sv = fromMaybe 0 (IM.lookup 0 (dat sv))++-- | 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+-- | 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+   

@@ -1,8 +1,16 @@ module Data.Sparse.Utils where +import qualified Data.Vector as V+ -- * Misc. utilities  +-- | Wrap a function with a null check, returning in Maybe+harness :: (t -> Bool) -> (t -> a) -> t -> Maybe a+harness q f v | q v = Nothing+              | otherwise = Just $ f v++ -- | Componentwise tuple operations -- TODO : use semilattice properties instead maxTup, minTup :: Ord t => (t, t) -> (t, t) -> (t, t)@@ -46,7 +54,7 @@   go _ _ [] = mneutral  --- *** Bounds checking+-- ** Bounds checking type LB = Int type UB = Int @@ -63,3 +71,16 @@  inBounds02 :: (UB, UB) -> (Int, Int) -> Bool inBounds02 (bx,by) (i,j) = inBounds0 bx i && inBounds0 by j+++++-- ** Safe indexing++++head' :: V.Vector a -> Maybe a+head' = harness V.null V.head++tail' :: V.Vector a -> Maybe (V.Vector a)+tail' = harness V.null V.tail

@@ -181,9 +181,16 @@   --- class (Set f, Sparse f a) => SparseSet f a --- instance SparseSet SpVector a where++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+  -- -- | Lookup with default, infix form ("safe" : should throw an exception if lookup is outside matrix bounds)+  (@@) :: c a -> ScIx c -> a++   

@@ -5,19 +5,21 @@          -- * Matrix factorizations          qr, lu,          chol,-         -- * Incomplete LU-         ilu0,          -- * Condition number          conditionNumberSM,          -- * Householder reflection          hhMat, hhRefl,          -- * Givens' rotation          givens,+         -- * Arnoldi iteration+         arnoldi,          -- * Eigensolvers          eigsQR, eigRayleigh,          -- * Linear solvers          linSolve, LinSolveMethod, (<\>),-         -- ** Methods+         -- ** Direct methods+         luSolve,+         -- ** Iterative methods          cgne, tfqmr, bicgstab, cgs, bcg,          _xCgne, _xTfq, _xBicgstab, _x, _xBcg,          cgsStep, bicgstabStep,@@ -35,7 +37,7 @@          -- * Sparsify data          sparsifySV,          -- * Iteration combinators-         modifyInspectN, runAppendN',+         modifyInspectN, runAppendN', untilConverged,          diffSqL        )        where@@ -65,7 +67,7 @@ -- import qualified Data.List as L import Data.Maybe -+import qualified Data.Vector as V   @@ -73,7 +75,7 @@ -- * Sparsify : remove almost-0 elements (|x| < eps) -- | Sparsify an SpVector sparsifySV :: Epsilon a => SpVector a -> SpVector a-sparsifySV (SV d im) = SV d $ IM.filter isNz im+sparsifySV = filterSV isNz   @@ -152,8 +154,12 @@  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) && isSquareSM mm =@@ -165,11 +171,6 @@     a = mm @@ (i', j)     b = mm @@ (i, j)   -- element to zero out --- |Is the `k`th the first nonzero column in the row?-firstNonZeroColumn :: IM.IntMap a -> IxRow -> Bool-firstNonZeroColumn mm k = isJust (IM.lookup k mm) &&-                          isNothing (IM.lookupLT k mm)- -- |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@@ -177,6 +178,11 @@   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)   @@ -382,7 +388,7 @@   -- | 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 :: Epsilon b => (Int -> b) -> [Int] -> [(Int, b)] onRangeSparse f ixs = filter (isNz . snd) $ zip ixs $ map f ixs  @@ -404,14 +410,14 @@   --- 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+-- -- 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   @@ -433,7 +439,7 @@   (l, u) = lu aa   lh = sparsifyLU l aa   uh = sparsifyLU u aa-  sparsifyLU m m2 = SM (dim m) $ ifilterIM2 f (dat m) where+  sparsifyLU m m2 = ifilterSM f m where     f i j _ = isJust (lookupSM m2 i j)  @@ -441,10 +447,50 @@   +-- * Arnoldi iteration +-- | Given a matrix A 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 ::+  (Floating a, Eq a) => SpMatrix a -> Int -> (SpMatrix a, SpMatrix a)+arnoldi aa kn = (fromCols qvfin, hhfin) where+  (qvfin, hhfin, _) = execState (modifyUntil tf arnoldiStep) arnInit +  tf (_, _, ii) = ii == kn -- termination criterion+  (m, n) = dim aa+  arnInit = (qv1, hh1, 1) where      +      q0 = normalize 2 $ onesSV n -- starting basis vector+      aq0 = aa #> q0+      h11 = q0 `dot` aq0+      q1nn = (aq0 ^-^ (h11 .* q0))+      hh1 = fromListSM (m + 1, n) [(0, 0, h11), (1, 0, h21)] where        +        h21 = norm 2 q1nn+      q1 = normalize 2 q1nn+      qv1 = V.fromList [q0, q1]+  arnoldiStep (qv, hh, i) = (qv', hh', i + 1)+   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 = singletonSV $ norm 2 qipnn      -- normalization factor H_{i+1, i}+    qip = normalize 2 qipnn              -- q_{i + 1}+    hh' = insertCol hh (concatSV (fromVector hhcoli) qipnorm) i -- update H+    qv' = V.snoc qv qip        -- append q_{i+1} to Krylov basis Q_i  +   +++++++++ -- * Preconditioning  -- | Partition a matrix into strictly subdiagonal, diagonal and strictly superdiagonal parts@@ -472,10 +518,46 @@   --- Linear solver, LU-based+-- * 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 = lubwSolve uu (lufwSolve ll b)+  | otherwise = error "luSolve : factors must be triangular matrices"  +lufwSolve ll b = sparsifySV v where+  (v, _) = execState (modifyUntil q lStep) lInit where+  q (_, i) = i == dim b+  lStep (ww, i) = (wwi, 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+    wwi = 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+lubwSolve uu w = sparsifySV x where+  (x, _) = execState (modifyUntil q uStep) uInit+  q (_, i) = i == (- 1)+  uStep (xx, i) = (xxi, 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+    xxi = 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)   

@@ -20,8 +20,6 @@ import Data.Sparse.Common  -import qualified Data.IntMap as IM- import Control.Monad (replicateM) import Control.Monad.State.Strict (execState) @@ -69,12 +67,14 @@     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+    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 : Linear solvers" $ do+  describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers" $ do     -- it "TFQMR (2 x 2 dense)" $     --   normSq (_xTfq (tfqmr aa0 b0 x0) ^-^ x0true) <= eps `shouldBe` True     it "BCG (2 x 2 dense)" $@@ -83,6 +83,9 @@       nearZero (normSq (aa0 <\> b0 ^-^ x0true)) `shouldBe` True     it "CGS (2 x 2 dense)" $        nearZero (normSq (_x (cgs aa0 b0 x0 x0) ^-^ x0true)) `shouldBe` True+  describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers" $ do+    it "LU (unoptimized) (4 x 4 sparse)" $ +      checkLuSolve aa1 b1 `shouldBe` True            describe "Numeric.LinearAlgebra.Sparse : QR decomposition" $ do         it "QR (4 x 4 sparse)" $       checkQr tm4 `shouldBe` True@@ -93,11 +96,63 @@       checkLu tm6 `shouldBe` True     it "LU (10 x 10 sparse)" $       checkLu tm7 `shouldBe` True-  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $ do+  describe "Numeric.LinearAlgebra.Sparse : Cholesky decomposition (PSD matrices only)" $      it "chol (5 x 5 sparse)" $       checkChol tm7 `shouldBe` True+  describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration" $ do+    it "Arnoldi iteration (3 x 3 dense)" $+      checkArnoldi aa2 3 `shouldBe` True+    it "Arnoldi iteration (5 x 5 sparse)" $+      checkArnoldi tm7 5 `shouldBe` True      ++{- QR-}+++checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool+checkQr a = c1 && c2 where+  (q, r) = qr a+  c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)+  c2 = isOrthogonalSM q++++{- LU -}++checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool+checkLu a = lup == a where+  (l, u) = lu a+  lup = l #~# u++++{- Cholesky -}++checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool+checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where+  l = chol a+++{- direct linear solver -}++checkLuSolve :: (Epsilon a, Real a, Floating a) => SpMatrix a -> SpVector a -> Bool+checkLuSolve amat rhs = nearZero (normSq ( (lmat #> (umat #> xlu)) ^-^ rhs ))+  where+     (lmat, umat) = lu amat+     xlu = luSolve lmat umat rhs+      +  +{- Arnoldi iteration -}+checkArnoldi :: (Epsilon a, Floating a, Eq a) => SpMatrix a -> Int -> Bool+checkArnoldi aa kn = nearZero $ normFrobenius $ (aa #~# qvprev) ^-^ (qv #~# hh) where+  (qv, hh) = arnoldi aa kn+  (m, n) = dim qv+  qvprev = extractSubmatrix qv (0, m - 1) (0, n - 2)++++ {-  example 0 : 2x2 linear system@@ -112,13 +167,9 @@  -} -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)]+aa0 :: SpMatrix Double+aa0 = fromListDenseSM 2 [1,3,2,4]   -- b0, x0 : r.h.s and initial solution resp.@@ -161,7 +212,9 @@ b2 = mkSpVectorD 3 [4,-2,4]  +aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double + -- --  {-@@ -262,41 +315,13 @@   -{- QR-}  -checkQr :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkQr a = c1 && c2 where-  (q, r) = qr a-  c1 = nearZero $ normFrobenius ((q #~# r) ^-^ a)-  c2 = isOrthogonalSM q  -aa22 = fromListDenseSM 2 [2,1,1,2] :: SpMatrix Double   --{- LU -}--checkLu :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkLu a = lup == a where-  (l, u) = lu a-  lup = l #~# u----{- Cholesky -}--checkChol :: (Epsilon a, Real a, Floating a) => SpMatrix a -> Bool-checkChol a = nearZero $ normFrobenius ((l ##^ l) ^-^ a) where-  l = chol a------ {- eigenvalues -}  @@ -326,8 +351,7 @@ tv0, tv1 :: SpVector Double tv0 = mkSpVectorD 2 [5, 6] --tv1 = SV 2 $ IM.singleton 0 1+tv1 = fromListSV 2 [(0,1)]   -- wikipedia test matrix for Givens rotation @@ -403,3 +427,10 @@ --       xhatC = head $ runNCGS niter aa b --   -- printDenseSM aa     --   return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))+++tm8 :: SpMatrix Double+tm8 = fromListSM (2,2) [(0,0,1), (0,1,1), (1,1,1)]++tm8' :: SpMatrix Double+tm8' = fromListSM (2,2) [(0,0,1), (1,0,1), (1,1,1)]