sparse-linear-algebra 0.2.9.7 → 0.2.9.8

raw patch · 23 files changed

+2012/−863 lines, 23 filesdep −vector-spacedep ~primitivePVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies removed: vector-space

Dependency ranges changed: primitive

API changes (from Hackage documentation)

Files

• CHANGELOG.markdown +12/−3
• README.md +66/−35
• sparse-linear-algebra.cabal +29/−13
• src/Control/Iterative.hs +9/−2
• src/Data/Sparse/Common.hs +69/−47
• src/Data/Sparse/Internal/CSB.hs +35/−77
• src/Data/Sparse/Internal/CSC.hs +125/−0
• src/Data/Sparse/Internal/IntM.hs +42/−32
• src/Data/Sparse/Internal/SList.hs +97/−0
• src/Data/Sparse/Internal/SVector.hs +310/−0
• src/Data/Sparse/Internal/SVector/Mutable.hs +119/−0
• src/Data/Sparse/Internal/Stream.hs +58/−0
• src/Data/Sparse/Internal/TriMatrix.hs +282/−0
• src/Data/Sparse/Internal/Utils.hs +3/−3
• src/Data/Sparse/PPrint.hs +89/−65
• src/Data/Sparse/SpMatrix.hs +49/−33
• src/Data/Sparse/SpVector.hs +29/−20
• src/Data/Sparse/Utils.hs +87/−18
• src/Numeric/LinearAlgebra/Class.hs +83/−129
• src/Numeric/LinearAlgebra/EigenSolvers/Experimental.hs +0/−1
• src/Numeric/LinearAlgebra/LinearSolvers/Experimental.hs +1/−1
• src/Numeric/LinearAlgebra/Sparse.hs +107/−105
• test/LibSpec.hs +311/−279

@@ -1,4 +1,13 @@ +	0.2.9.7+	------++	Improved pretty printer:++	* Fixed display precision (e.g. 2 decimal digits), fixed column width output for vectors and matrices+	+	* Small and large values (wrt fixed precision) switch to scientific notation+	 	0.2.9.4 	------- @@ -6,10 +15,10 @@   	0.2.9.1- 	-------  	* Uses classes from `vector-space` such as AdditiveGroup, VectorSpace and InnerSpace+	 	* QuickCheck tests for algebraic properties, such as matrix-vector products and soon more abstract ones e.g. positive semi-definite matrices-	* Getting rid of `error` in favor of MonadThrow exceptions for high-level operations-such as matrix algorithms+	+	* Getting rid of `error` in favor of MonadThrow exceptions for high-level operations such as matrix algorithms

@@ -4,10 +4,22 @@  [![Hackage](https://img.shields.io/hackage/v/sparse-linear-algebra.svg)](https://hackage.haskell.org/package/sparse-linear-algebra)  [![Build Status](https://travis-ci.org/ocramz/sparse-linear-algebra.png)](https://travis-ci.org/ocramz/sparse-linear-algebra) -This library provides common numerical analysis functionality, without requiring any external bindings. It is not optimized for performance (yet), but it serves as an experimental platform for scientific computation in a purely functional setting.+This library provides common numerical analysis functionality, without requiring any external bindings. It aims to serve as an experimental platform for scientific computation in a purely functional setting. -Contents :+## News +Oct 7., 2017: The library is evolving in a number of ways, to reflect performance observations and user requests:++* typeclasses and instances for primitive types will become `sparse-linear-algebra-core`, along with a typeclass-oriented reformulation of the numerical algorithms that used to depend on the nested IntMap representation.+This will let other developers build on top of this library, in the spirit of `vector-space` and `linear`.++* The `vector`-based backend is being reworked.++* An `accelerate`-based backend is under development.+++## Contents+ * Iterative linear solvers (`<\>`)      * Generalized Minimal Residual (GMRES) (non-Hermitian systems) @@ -96,56 +108,68 @@ Both sparse vectors and matrices can be pretty-printed using `prd`:      λ> prd amat-    ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 ) -    2.0  _   _ -    4.0 3.0 2.0-     _   _  5.0+    ( 3 rows, 3 columns ) , 5 NZ ( density 55.556 % ) -Note: sparse data should only contain non-zero entries not to waste memory and computation.+    2.00   , _      , _      +    4.00   , 3.00   , 2.00   +    _      , _      , 5.00        -### Matrix operations+*Note (sparse storage)*: sparse data should only contain non-zero entries not to waste memory and computation. +*Note (approximate output)*: `prd` rounds the results to two significant digits, and switches to scientific notation for large or small values. Moreover, values which are indistinguishable from 0 (see the `Numeric.Eps` module) are printed as `_`. +++### Matrix factorizations, matrix product+ There are a few common matrix factorizations available; in the following example we compute the LU factorization of matrix `amat` and verify it with the matrix-matrix product `##` of its factors :      λ> (l, u) <- lu amat     λ> prd $ l ## u-    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )+    +    ( 3 rows, 3 columns ) , 9 NZ ( density 100.000 % ) -    2.0  _   _ -    4.0 3.0 2.0-     _   _  5.0+    2.00   , _      , _      +    4.00   , 3.00   , 2.00   +    _      , _      , 5.00        + Notice that the result is _dense_, i.e. certain entries are numerically zero but have been inserted into the result along with all the others (thus taking up memory!). To preserve sparsity, we can use a sparsifying matrix-matrix product `#~#`, which filters out all the elements x for which `|x| <= eps`, where `eps` (defined in `Numeric.Eps`) depends on the numerical type used (e.g. it is 10^-6 for `Float`s and 10^-12 for `Double`s).      λ> prd $ l #~# u-    ( 3 rows, 3 columns ) , 5 NZ ( sparsity 0.5555555555555556 )+    +    ( 3 rows, 3 columns ) , 5 NZ ( density 55.556 % ) -    2.0  _   _ -    4.0 3.0 2.0-     _   _  5.0    +    2.00   , _      , _      +    4.00   , 3.00   , 2.00   +    _      , _      , 5.00  + A matrix is transposed using the `transpose` function.  Sometimes we need to compute matrix-matrix transpose products, which is why the library offers the infix operators `#^#` (i.e. matrix transpose * matrix) and `##^` (matrix * matrix transpose):      λ> amat' = amat #^# amat     λ> prd amat'-    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )+    +    ( 3 rows, 3 columns ) , 9 NZ ( density 100.000 % ) -    20.0 12.0 8.0-    12.0 9.0 6.0-    8.0 6.0 29.0+    20.00  , 12.00  , 8.00   +    12.00  , 9.00   , 6.00   +    8.00   , 6.00   , 29.00      +     -    λ> l <- chol amat'-    λ> prd $ l ##^ l-    ( 3 rows, 3 columns ) , 9 NZ ( sparsity 1.0 )+    λ> lc <- chol amat'+    λ> prd $ lc ##^ lc+    +    ( 3 rows, 3 columns ) , 9 NZ ( density 100.000 % ) -    20.000000000000004 12.0 8.0-    12.0 9.0 10.8-    8.0 10.8 29.0+    20.00  , 12.00  , 8.00   +    12.00  , 9.00   , 10.80  +    8.00   , 10.80  , 29.00       + In the last example we have also shown the Cholesky decomposition (M = L L^T where L is a lower-triangular matrix), which is only defined for symmetric positive-definite matrices.  ### Linear systems@@ -155,25 +179,29 @@     λ> b = fromListDenseSV 3 [3,2,5] :: SpVector Double     λ> x <- amat <\> b     λ> prd x-    ( 3 elements ) ,  3 NZ ( sparsity 1.0 ) -    1.4999999999999998 -1.9999999999999998 0.9999999999999998+    ( 3 elements ) ,  3 NZ ( density 100.000 % ) +    1.50   , -2.00  , 1.00      ++ The result can be verified by computing the matrix-vector action `amat #> x`, which should (ideally) be very close to the right-hand side `b` :      λ> prd $ amat #> x-    ( 3 elements ) ,  3 NZ ( sparsity 1.0 ) -    2.9999999999999996 1.9999999999999996 4.999999999999999+    ( 3 elements ) ,  3 NZ ( density 100.000 % ) -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:+    3.00   , 2.00   , 5.00       +     +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 LU factors defined previously we can cross-verify the two solution methods:+     λ> x' <- luSolve l u b     λ> prd x' -    ( 3 elements ) ,  3 NZ ( sparsity 1.0 )+    ( 3 elements ) ,  3 NZ ( density 100.000 % ) -    1.5 -2.0 1.0+    1.50   , -2.00  , 1.00        @@ -182,6 +210,7 @@   + ## License  GPL3, see LICENSE@@ -200,6 +229,8 @@  [2] G.H. Golub and C.F. Van Loan, Matrix Computations, 3rd ed., 1996 -[3] L.N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997+[3] T.A. Davis, Direct Methods for Sparse Linear Systems, 2006 -[4] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77, 2nd ed., 1992 +[4] L.N. Trefethen, D. Bau, Numerical Linear Algebra, SIAM, 1997++[5] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in Fortran 77, 2nd ed., 1992 

@@ -1,5 +1,5 @@ name:                sparse-linear-algebra-version:             0.2.9.7+version:             0.2.9.8 synopsis:            Numerical computation in native Haskell description:   /Overview/@@ -26,11 +26,23 @@ data-files:          e05r0000.mtx                      e05r0000_rhs1.mtx cabal-version:       >=1.10-tested-with:         GHC == 8.0.1+tested-with:         GHC == 8.0.2 +source-repository head+  type:     git+  location: https://github.com/ocramz/sparse-linear-algebra+                     +flag dump+  manual: True+  default: False+ library   default-language:    Haskell2010-  ghc-options:         -Wall+  ghc-options:         -Wall -O2++  if flag(dump)+    ghc-options: -ddump-simpl -ddump-stg -ddump-to-file+     hs-source-dirs:      src   exposed-modules:     Numeric.LinearAlgebra.Sparse                        Numeric.LinearAlgebra.Class@@ -42,7 +54,14 @@   other-modules:       Data.Sparse.Internal.IntMap2                        Data.Sparse.Internal.IntM                        Data.Sparse.Internal.Utils-                       Data.Sparse.Internal.CSB                       +                       Data.Sparse.Internal.SVector+                       Data.Sparse.Internal.SVector.Mutable+                       Data.Sparse.Internal.Stream                   +                       -- Data.Sparse.Internal.SHVector+                       Data.Sparse.Internal.SList+                       Data.Sparse.Internal.TriMatrix+                       Data.Sparse.Internal.CSB+                       Data.Sparse.Internal.CSC                                               Data.Sparse.Utils                        Data.Sparse.PPrint                        Data.Sparse.Types@@ -50,18 +69,20 @@                        Numeric.LinearAlgebra.LinearSolvers.Experimental                        Numeric.LinearAlgebra.EigenSolvers.Experimental   build-depends:       base >= 4.7 && < 5-                     , vector-space+                     , primitive+                     -- , deepseq                      , containers                      , vector                      , vector-algorithms+                     -- , hybrid-vectors                      , exceptions                      , mtl >= 2.2.1                                           , transformers                      -- , writer-cps-transformers                      -- , monad-log == 0.1.1.0                      -- , monad-par-                     -- , QuickCheck >= 2.8.2-                     -- , hspec+                     , QuickCheck >= 2.8.2+                     , hspec                      -- , primitive >= 0.6.1.0                      -- , mwc-random                      @@ -76,8 +97,6 @@ --                      , matrix-market-attoparsec --                      , scientific   -- test-suite spec   default-language:    Haskell2010   ghc-options:         -Wall@@ -95,13 +114,10 @@                      , mwc-random                      , primitive >= 0.6.1.0                      , sparse-linear-algebra-                     , vector-space                      , scientific                      , matrix-market-attoparsec -source-repository head-  type:     git-  location: https://github.com/ocramz/sparse-linear-algebra+   

@@ -29,9 +29,9 @@ import qualified Control.Monad.Trans.State.Strict  as MTS -- (runStateT) import Data.Foldable (foldrM) -import Data.VectorSpace  + data ConvergenceStatus a = BufferNotReady                          | Converging                          | Converged a@@ -78,7 +78,7 @@    x <- get    y <- lift $ f x    when (printDebugInfo config) $ liftIO $ do-     putStrLn $ unwords ["Iteration", show i]+     putStrLn $ unwords ["Iteration", show i, "\n"]      printDebugIO config (pf y)     put y    if q y@@ -247,6 +247,13 @@   -- | Helpers++-- | Relative residual+relRes :: (Normed t, LinearVectorSpace t) =>+     MatrixType t -> t -> t -> Magnitude t+relRes aa b x = n / d where+  n = norm2 $ (aa #> x) ^-^ b+  d = norm2 b   -- meanl :: (Foldable t, Fractional a) => t a -> a

@@ -51,8 +51,13 @@ import Data.Maybe (fromMaybe, maybe) import qualified Data.Vector as V -import Data.VectorSpace ++++++ -- withBoundsSM m ij e f --   | isValidIxSM m ij = f m ij --   | otherwise = error e@@ -230,15 +235,15 @@   :: (Normed f1, Num b, Functor f) => f (f1 b) -> f1 b -> f b -} -instance LinearVectorSpace (SpVector Double) where-  type MatrixType (SpVector Double) = SpMatrix Double-  (#>) = matVecSD-  (<#) = vecMatSD+-- instance LinearVectorSpace (SpVector Double) where+--   type MatrixType (SpVector Double) = SpMatrix Double+--   (#>) = matVecSD+--   (<#) = vecMatSD -instance LinearVectorSpace (SpVector (Complex Double)) where-  type MatrixType (SpVector (Complex Double)) = SpMatrix (Complex Double)-  (#>) = matVecSD-  (<#) = vecMatSD+-- instance LinearVectorSpace (SpVector (Complex Double)) where+--   type MatrixType (SpVector (Complex Double)) = SpMatrix (Complex Double)+--   (#>) = matVecSD+--   (<#) = vecMatSD     @@ -342,75 +347,92 @@  toDenseRow :: Num a => SpMatrix a -> IM.Key -> [a] toDenseRow sm irow =-  fmap (\icol -> sm @@ (irow,icol)) [0..ncol-1] where (_, ncol) = (nrows sm, ncols sm)+  fmap (\icol -> sm @@ (irow,icol)) [0..ncols sm-1]   -toDenseRowClip :: (Show a, Epsilon a) => SpMatrix a -> IM.Key -> Int -> String-toDenseRowClip sm irow ncomax-  | nco > ncomax = unwords (map showNz h) ++  " ... " ++ showNz t-  | otherwise = unwords $ showNz <$> dr-     where dr = toDenseRow sm irow-           h = take (ncomax - 2) dr-           t = last dr-           (_, nco) = dim sm+prdR, prdC :: PPrintOptions+prdR = PPOpts 4 2 7   -- reals+prdC = PPOpts 4 2 16  -- complex values  --- printDenseSM :: (Show t, Num t) => SpMatrix t -> IO ()++-- -- printDenseSM :: (Show t, Num t) => SpMatrix t -> IO () -- printDenseSM :: (ScIx c ~ (Int, Int), FDSize c ~ (Int, Int), SpContainer c a, Show a, Epsilon a) => c a -> IO () printDenseSM sm = do   newline-  putStrLn $ sizeStr sm+  putStrLn $ sizeStrSM sm   newline-  printDenseSM0 sm+  prd0 sm -printDenseSM0 sm = do+++printDenseSM0+  :: (SpMatrix a -> IxRow -> Int -> String) -> SpMatrix a -> IO ()+printDenseSM0 f sm = do   printDenseSM' sm 5 5   newline-  where    -    -- printDenseSM' :: (Show t, Num t) => SpMatrix t -> Int -> Int -> IO ()+  where     printDenseSM' sm' nromax ncomax = mapM_ putStrLn rr_' where       (nr, _) = (nrows sm, ncols sm)-      rr_ = map (\i -> toDenseRowClip sm' i ncomax) [0..nr - 1]+      rr_ = map (\i -> f sm' i ncomax) [0..nr-1]       rr_' | nr > nromax = take (nromax - 2) rr_ ++ [" ... "] ++[last rr_]-           | otherwise = rr_+           | otherwise = rr_  +printDenseSM0r :: SpMatrix Double -> IO ()+printDenseSM0r sm = printDenseSM0 g sm+  where+    g sm' irow ncolmax = printDN (ncols sm') ncolmax prdR $ toDenseRow sm' irow+printDenseSM0c :: SpMatrix (Complex Double) -> IO ()+printDenseSM0c sm = printDenseSM0 g sm+  where+    g sm' irow ncolmax = printCN (ncols sm') ncolmax prdC $ toDenseRow sm' irow  -toDenseListClip :: (Show a, Epsilon a) => SpVector a -> Int -> String-toDenseListClip sv ncomax-  | dim sv > ncomax = unwords (map showNz h) ++  " ... " ++ showNz t-  | otherwise = unwords $ showNz <$> dr-     where dr = toDenseListSV sv-           h = take (ncomax - 2) dr-           t = last dr--printDenseSV :: (Show t, Epsilon t) => SpVector t -> IO ()+-- printDenseSV :: (Show t, Epsilon t) => SpVector t -> IO ()+printDenseSV :: PrintDense (SpVector a) => SpVector a -> IO () printDenseSV sv = do   newline   putStrLn $ sizeStrSV sv   newline-  printDenseSV0 sv+  prd0 sv -printDenseSV0 :: (Show t, Epsilon t) => SpVector t -> IO ()-printDenseSV0 sv = do-  printDenseSV' sv 5+printDenseSV0r :: SpVector Double -> IO ()+printDenseSV0r = printDenseSV0 g where+  g l n = printDN l n prdR+printDenseSV0c :: SpVector (Complex Double) -> IO ()+printDenseSV0c = printDenseSV0 g where+  g l n = printCN l n prdC++printDenseSV0 :: Num a =>+  (Int -> Int -> [a] -> String) -> SpVector a -> IO ()+printDenseSV0 f sv = do+  printDenseSV' (svDim sv) 5   newline where-    printDenseSV' v nco = putStrLn rr_' where-      rr_ = toDenseListClip v nco :: String-      rr_' | dim sv > nco = unwords [take (nco - 2) rr_ , " ... " , [last rr_]]-           | otherwise = rr_+    printDenseSV' l n = putStrLn (f l n vd)+    vd = toDenseListSV sv ++     -- ** Pretty printer typeclass +instance PrintDense (SpVector Double) where+  prd = printDenseSV+  prd0 = printDenseSV0r -instance (Show a, Epsilon a) => PrintDense (SpVector a) where+instance PrintDense (SpVector (Complex Double)) where   prd = printDenseSV-  prd0 = printDenseSV0+  prd0 = printDenseSV0c -instance (Show a, Epsilon a) => PrintDense (SpMatrix a) where+instance PrintDense (SpMatrix Double) where   prd = printDenseSM-  prd0 = printDenseSM0+  prd0 = printDenseSM0r++instance PrintDense (SpMatrix (Complex Double)) where+  prd = printDenseSM+  prd0 = printDenseSM0c++   

@@ -5,7 +5,7 @@ -- import Control.Monad.Primitive import Control.Monad.ST -import qualified Data.Foldable as F -- (foldl')+import Data.Foldable (foldl') -- import Data.List (group, groupBy)  import qualified Data.Vector as V @@ -14,12 +14,14 @@ import qualified Data.Vector.Algorithms.Merge as VA (sortBy) -- import qualified Data.Vector.Generic as VG (convert) +import qualified Data.IntMap as IM+ import Control.Monad import Data.Maybe import Data.Ord (comparing) -import Data.Complex-import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)+-- import Data.Complex+-- import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)  import Control.Concurrent -- import qualified Control.Monad.Par as Par@@ -61,93 +63,49 @@ -- 2) " have at most (bRows x bCols) NZ -- data CsbMatrix a = CSB {-  csbNrows, csbNcols :: {-# UNPACK #-} !Int,+  csbNrows, csbNcols, csbBeta :: {-# UNPACK #-} !Int,   csbVal :: V.Vector a,-  csbRowPtr, csbRowIx, csbColIx :: V.Vector Int+  csbBlkPtr, csbRowIx, csbColIx :: V.Vector Int                        } deriving (Eq, Functor)  -{- |-encode:--Given:-* Block size-* Matrix size (#rows, #cols)-* (row, column) index pair--Compute:-* Block index-* (Row, column) index (relative to block)--}------ data Block a = Block {---   bRows, bCols :: {-# UNPACK #-} !Int,---   rowIx :: V.Vector Int,---   colIx :: V.Vector Int,---   val :: V.Vector a      -- invariant: ---   } deriving (Eq, Show)---- instance Functor Block where---   fmap f (Block br bc ri ci x) = Block br bc ri ci (fmap f x)---- instance Foldable Block where---   foldr f z bl = foldr f z (val bl)---- instance Traversable Block where---   traverse f (Block br bc ri ci v) = Block br bc ri ci <$> traverse f v----------- -- * CSB Matrix---- -- | Invariants:--- -- 1) csbNrows = sum bRows, csbNcols = sum bCols---- data CsbMatrix a =---   CB { csbNrows :: {-# UNPACK #-} !Int,---        csbNcols :: {-# UNPACK #-} !Int,---        -- csbBlkPtr :: V.Vector Int,---        csbBlocks :: V.Vector (Block a) -- block ordering function f(i,j)---      } deriving Eq---- instance Functor CsbMatrix where---   fmap f (CB m n bb) = CB m n $ fmap (fmap f) bb---- -- instance Foldable CsbMatrix where--- --   foldr f z cm = foldr f z (csbBlocks cm)---- -- instance Traversable CsbMatrix where--- --   traverse f (CB m n bp x) = CB m n bp <$> traverse f x----+csbParams :: (Int, Int)  -- ^ Matrix size+          -> Int         -- ^ Block parameter (i.e. block edge length)+          -> (Int, Int)  -- ^ # of blocks per side+csbParams (m,n) beta = (ceiling m', ceiling n')+  where m' = fromIntegral m / fromIntegral beta+        n' = fromIntegral n / fromIntegral beta --- * fromList--- m n [(i, j, x)] -> C+-- | Which block (in row, col format) do matrix indices (i,j) fall in ?+bBin, bCoords :: Integral t => t -> t -> t -> (t, t)+bBin beta i j = (i `div` beta, j `div` beta) +-- | Block coordinates+bCoords beta i j = (i `mod` beta, j `mod` beta) +-- | Block index (row-major order)+blockIx :: (Int, Int) -> Int -> Int -> Int -> Int+blockIx dims beta i j = bx + by*nbx where+  (bx, by) = bBin beta i j+  (nbx, _) = csbParams dims beta   +type Block i a = [(i,i,a)] --- * Mat-vec (y = A x)--- y[i] = < A[i,:] , x >+{-# inline consBlockElem #-}+consBlockElem ::+  IM.Key -> (i, i, a) -> IM.IntMap (Block i a) -> IM.IntMap (Block i a)+consBlockElem ib x bb = IM.insert ib (x : blocki) bb where+  blocki = fromMaybe [] (IM.lookup ib bb)  --- * MatT-vec (y = A^T x)--- y[i] = < A[:, j] , x>-+consBlocks :: Foldable t =>+     (Int, Int) -> Int -> t (Int, Int, a) -> IM.IntMap (Block Int a)+consBlocks dims beta ijx = foldl' ins IM.empty ijx where+  ins accb (i,j,x) = consBlockElem ib (i,j,x) accb where+    ib = blockIx dims beta i j +     --- * Mat-vec (y = A x, parallel over CSB and block-partitioned x)--- | Invariants:--- 1) # cols block Aij == # entries j'th partition of x  

@@ -0,0 +1,125 @@+{-# language TypeFamilies, FlexibleInstances, MultiParamTypeClasses #-}+module Data.Sparse.Internal.CSC where++import Control.Monad (forM_, when)+import Control.Monad.Primitive (PrimMonad(..), PrimState(..))++import qualified Data.Graph as G+import qualified Data.Vector as V+import qualified Data.Vector.Mutable as VM++-- import Data.Sparse.Types+import Data.Sparse.Internal.SVector+import Data.Sparse.Internal.Utils+import Numeric.LinearAlgebra.Class+++data CscMatrix a =+  CscM {+      cscNrows :: {-# UNPACK #-} !Int,+      cscNcols :: {-# UNPACK #-} !Int,+      cscNz :: {-# UNPACK #-} !Int,+      cscRowIx :: V.Vector Int,+      cscColPtr :: V.Vector Int,+      cscVal :: V.Vector a} deriving Eq++instance Functor CscMatrix where+  fmap f (CscM m n nz cc rp x) = CscM m n nz cc rp (fmap f x)++instance Foldable CscMatrix where+  foldr f z cm = foldr f z (cscVal cm)++instance Show a => Show (CscMatrix a) where+  show m'@(CscM m n nz cix rp x) = szs where+    szs = unwords ["CSC (",show m, "x", show n,"),",show nz, "NZ ( sparsity",show (spy m'),"), row indices:",show cix,", column pointers:", show rp,", data:",show x]++instance FiniteDim (CscMatrix a) where+  type FDSize (CscMatrix a) = (Int, Int)+  dim m = (cscNrows m, cscNcols m)++instance HasData (CscMatrix a) where+  nnz = cscNz+  +instance Sparse (CscMatrix a) where+  spy m = fromIntegral (nnz m) / fromIntegral (cscNrows m * cscNcols m)++  +++-- * Creation+-- | O(N log N) : Copy a Vector containing (row index, column index, entry) into a CSC structure. Sorts the Vector by row columns ( O(log N) ), unzips row indices and data ( O(N) ) and generates the column pointer vector ( 2 O(N) passes )+toCSC :: Int -> Int -> V.Vector (Int, Int, a) -> CscMatrix a+toCSC m n ijxv = CscM m n nz rix crp x where+  nz = V.length x+  (rix, cix, x) = V.unzip3 (sortWith snd3 ijxv)  -- sort by columns+  crp = csPtrV (==) m cix++-- | O(N) : Rebuilds the (row, column, entry) Vector from the CSC representation. +fromCSC :: CscMatrix a -> V.Vector (Int, Int, a)+fromCSC mc = V.zip3 ii jj xx where (ii,jj,xx) = fromCSC0 mc++fromCSC0 :: CscMatrix a -> (V.Vector Int, V.Vector Int, V.Vector a)+fromCSC0 mc = (rowIx, cols, cscVal mc) where+  (_, n) = dim mc+  rowIx = cscRowIx mc+  l = length rowIx+  cp = cscColPtr mc+  cols = V.create $ do+    colv <- VM.replicate l 0+    forM_ [0 .. n-1] (\i -> go colv i 0)+    return colv+  go :: PrimMonad m => VM.MVector (PrimState m) Int -> Int -> Int -> m ()  +  go vm irow j = when (j < nj) $ do+                          VM.write vm (j + jmin) irow+                          go vm irow (succ j) where+    jmin = cp V.! irow+    jmax = cp V.! (irow + 1)+    nj = jmax - jmin+++-- ** Extract a column+-- | O(1) : extract a column from the CSC matrix.+extractColCSC :: CscMatrix a -> Int -> SVector a+extractColCSC (CscM m _ _ rix cp x) jcol = SV m ixs vals where+  jmin = cp V.! jcol+  jmax = cp V.! (jcol + 1)+  ixs = V.slice jmin (jmax - jmin) rix+  vals = V.slice jmin (jmax - jmin) x++-- ** Extract the diagonal element, if this exists (NB this only holds for a lower triangular matrix)+extractDiagSubdiagCSC :: CscMatrix a -> Int -> Maybe (a, SVector a)+extractDiagSubdiagCSC cscm j+  | V.length ixs > 1 && j == V.head ixs = Just (diagEl, subdiagVec)+  | otherwise = Nothing where+      (SV m ixs vals) = extractColCSC cscm j+      diagEl = V.head vals+      subdiagVec = SV (m-1) (V.tail ixs) (V.tail vals) +++++++-- | O(N log N) : Transpose CSC matrix+transposeCSC :: CscMatrix a -> CscMatrix a+transposeCSC mm = toCSC n m $ V.zip3 jj ii xx where+  (m,n) = dim mm+  (ii, jj, xx) = fromCSC0 mm+++-- * Helpers++cscToGraph :: CscMatrix a -> G.Graph+cscToGraph ll = G.buildG (0, m-1) $ V.toList (V.zip ill jll) -- graph of L+  where+   (m, _) = dim ll -- dimensions of L = bounds of G(L)+   (ill, jll, _) = fromCSC0 ll+++++-- -- example data++-- -- row = np.array([0, 0, 1, 2, 2, 2])+-- -- col = np.array([0, 2, 2, 0, 1, 2])+-- -- data = np.array([1, 2, 3, 4, 5, 6])

@@ -6,7 +6,7 @@  import GHC.Exts import Data.Complex-import Data.VectorSpace+ import qualified Data.IntMap as IM  @@ -16,38 +16,48 @@ -- | A synonym for IntMap  newtype IntM a = IntM {unIM :: IM.IntMap a} deriving (Eq, Show, Functor, Foldable) -+empty :: IntM a empty = IntM IM.empty +size :: IntM a -> Int size (IntM x) = IM.size x +singleton :: IM.Key -> a -> IntM a singleton i x = IntM $ IM.singleton i x  filterWithKey f im = IntM $ IM.filterWithKey f (unIM im) +insert :: IM.Key -> a -> IntM a -> IntM a insert k x (IntM im) = IntM $ IM.insert k x im +filterI :: (a -> Bool) -> IntM a -> IntM a filterI f (IntM im) = IntM $ IM.filter f im -lookup x (IntM im) = IM.lookup x im+lookup :: IM.Key -> IntM a -> Maybe a+lookup i (IntM im) = IM.lookup i im  lookupLT x (IntM im) = IM.lookupLT x im +foldlWithKey :: (a -> IM.Key -> b -> a) -> a -> IntM b -> a foldlWithKey f z (IntM im) = IM.foldlWithKey f z im +foldlWithKey' :: (a -> IM.Key -> b -> a) -> a -> IntM b -> a foldlWithKey' f z (IntM im) = IM.foldlWithKey' f z im  mapWithKey f (IntM im) = IntM $ IM.mapWithKey f im +keys :: IntM a -> [IM.Key] keys (IntM im) = IM.keys im  mapKeys f (IntM im) = IntM $ IM.mapKeys f im +union :: IntM a -> IntM a -> IntM a union (IntM a) (IntM b) = IntM $ IM.union a b  findMin (IntM im) = IM.findMin im findMax (IntM im) = IM.findMax im +(!) :: IntM a -> IM.Key -> a (IntM im) ! i = im IM.! i  @@ -75,45 +85,45 @@   --- | ParamInstance can be used with all types that are instances of Set (which are by construction also instances of Functor)-#define ParamInstance(f, t) \-  instance VectorSpace (f t) where {type (Scalar (f (t))) = (t); n *^ im = fmap (* n) im};\-  instance VectorSpace (f (Complex t)) where {type (Scalar (f (Complex t))) = Complex (t); n *^ im = fmap (* n) im};\-  instance InnerSpace (f t) where {a <.> b = sum $ liftI2 (*) a b};\-  instance InnerSpace (f (Complex t)) where {a <.> b = sum $ liftI2 (*) (conjugate <$> a) b};\-  -- instance Normed (f t) where {type RealScalar (f t) = t ; type Magnitude (f t) = t ; norm1 a = sum (abs <$> a) ; norm2Sq a = sum $ liftI2 (*) a a; normP p v = sum u**(1/p) where u = fmap (**p) v; normalize = normzPR ; normalize2 = normz2R}; \-  -- instance Normed (f (Complex t)) where {type RealScalar (f (Complex t)) = t; type Magnitude (f (Complex t)) = t; norm1 a = realPart $ sum (abs <$> a); norm2Sq a = realPart $ sum $ liftI2 (*) (conjugate <$> a) a; normP p v = realPart $ sum u**(1/(p :+ 0)) where u = fmap (**(p :+ 0)) v; normalize = normzPC; normalize2 = normz2C }+-- -- | ParamInstance can be used with all types that are instances of Set (which are by construction also instances of Functor)+-- #define ParamInstance(f, t) \+--   instance VectorSpace (f t) where {type (Scalar (f (t))) = (t); n .* im = fmap (* n) im};\+--   instance VectorSpace (f (Complex t)) where {type (Scalar (f (Complex t))) = Complex (t); n .* im = fmap (* n) im};\+--   instance InnerSpace (f t) where {a <.> b = sum $ liftI2 (*) a b};\+--   instance InnerSpace (f (Complex t)) where {a <.> b = sum $ liftI2 (*) (conjugate <$> a) b};\+--   -- instance Normed (f t) where {type RealScalar (f t) = t ; type Magnitude (f t) = t ; norm1 a = sum (abs <$> a) ; norm2Sq a = sum $ liftI2 (*) a a; normP p v = sum u**(1/p) where u = fmap (**p) v; normalize = normzPR ; normalize2 = normz2R}; \+--   -- instance Normed (f (Complex t)) where {type RealScalar (f (Complex t)) = t; type Magnitude (f (Complex t)) = t; norm1 a = realPart $ sum (abs <$> a); norm2Sq a = realPart $ sum $ liftI2 (*) (conjugate <$> a) a; normP p v = realPart $ sum u**(1/(p :+ 0)) where u = fmap (**(p :+ 0)) v; normalize = normzPC; normalize2 = normz2C }  -instance Normed (IntM Double) where-  type RealScalar (IntM Double) = Double-  type Magnitude (IntM Double) = Double-  norm1 a = sum (abs <$> a)-  norm2Sq a = sum $ liftI2 (*) a a-  normP p v = sum u**(1/p) where u = fmap (**p) v-  normalize p v = v ./ normP p v -  normalize2 v = v ./ norm2 v +-- instance Normed (IntM Double) where+--   type RealScalar (IntM Double) = Double+--   type Magnitude (IntM Double) = Double+--   norm1 a = sum (abs <$> a)+--   norm2Sq a = sum $ liftI2 (*) a a+--   normP p v = sum u**(1/p) where u = fmap (**p) v+--   normalize p v = v ./ normP p v +--   normalize2 v = v ./ norm2 v    -instance Normed (IntM (Complex Double)) where-  type RealScalar (IntM (Complex Double)) = Double-  type Magnitude (IntM (Complex Double)) = Double-  norm1 a = realPart $ sum (abs <$> a)-  norm2Sq a = realPart $ sum $ liftI2 (*) (conjugate <$> a) a-  normP p v = realPart $ sum u**(1/(p :+ 0)) where u = fmap (**(p :+ 0)) v-  normalize p v = v ./ toC (normP p v)-  normalize2 v = v ./ toC (norm2 v)+-- instance Normed (IntM (Complex Double)) where+--   type RealScalar (IntM (Complex Double)) = Double+--   type Magnitude (IntM (Complex Double)) = Double+--   norm1 a = realPart $ sum (abs <$> a)+--   norm2Sq a = realPart $ sum $ liftI2 (*) (conjugate <$> a) a+--   normP p v = realPart $ sum u**(1/(p :+ 0)) where u = fmap (**(p :+ 0)) v+--   normalize p v = v ./ toC (normP p v)+--   normalize2 v = v ./ toC (norm2 v)      --- | IntMap instances-#define IntMapInstance(t) \-  ParamInstance( IntM, t )+-- -- | IntMap instances+-- #define IntMapInstance(t) \+--   ParamInstance( IntM, t ) -IntMapInstance(Double)--- IntMapInstance(Float)+-- IntMapInstance(Double)+-- -- IntMapInstance(Float)   

@@ -0,0 +1,97 @@+{-# language TypeFamilies, FlexibleInstances, DeriveFunctor #-}+module Data.Sparse.Internal.SList where+++import Data.Sparse.Utils+import Numeric.LinearAlgebra.Class+++++-- | Sparse list+newtype SList a = SL {unSL :: [(Int, a)]} deriving (Eq, Show, Functor)++emptySL :: SList a+emptySL = SL []++consSL :: (Int, a) -> SList a -> SList a+consSL x (SL xs) = SL (x : xs)++headSL :: SList a -> Maybe (Int, a)+headSL (SL (x:_)) = Just x+headSL (SL []) = Nothing++fromList :: [(Int, a)] -> SList a+fromList = SL++toList :: SList a -> [(Int, a)]+toList = unSL+++{-|+NB : unionWith and intersectWith work only if the indices are _sorted_++NB2 : we use the _descending_ order comparison+-}+++-- | Inner product between sparse lists+sldot :: (Elt a, Ord i) => [(i, a)] -> [(i, a)] -> a+sldot u v = sum $ intersectWithD pf u v where+  pf x y = conj x * y++-- | Vector sum of sparse lists+slsum :: (Ord i, Elt a) => [(i, a)] -> [(i, a)] -> [(i, a)]+slsum = unionWithD (+) 0++++++++++-- | `vector-space` instances++instance Elt a => AdditiveGroup (SList a) where+  zeroV = SL []+  negateV = fmap (* (-1))+  u ^+^ v = SL $ slsum (unSL u) (unSL v)+++instance Elt a => VectorSpace (SList a) where+  type Scalar (SList a) = a+  a .* v = fmap (* a) v+++instance (AdditiveGroup a, Elt a) => InnerSpace (SList a) where+  u <.> v = sldot (unSL u) (unSL v)  +++-- instance InnerSpace (SList Double) where+--   u <.> v = inner (unSV u) (unSV v)++-- instance InnerSpace (SList (Complex Double)) where+--   u <.> v = innerC (unSV u) (unSV v)+++++-- test data++l1, l2 :: [(Int, Double)]+l1 = [(0, pi), (2, pi),  (3, 5.4) ]++l2 = [(1, exp 1), (2, 3.4)]++-- l1c :: [(Int, Complex Double)]+-- l1c = zip ii $ zipWith (:+) [1..3] [3,2,1] where+--   ii = [0, 2, 5]++-- sl1c = SL l1c+++-- helpers+-- sortIndices :: [(IM.Key, a)] -> [(IM.Key, a)]+-- sortIndices = IM.toList . IM.fromList

@@ -0,0 +1,310 @@+{-# language TypeFamilies, FlexibleInstances, MultiParamTypeClasses, CPP #-}+module Data.Sparse.Internal.SVector where++import Control.Arrow+import Control.Monad (unless)+import Control.Monad.IO.Class+import qualified Data.Foldable as F -- (foldl')+-- import Data.List (group, groupBy)++import qualified Data.Vector as V +-- import qualified Data.Vector.Unboxed as VU+import qualified Data.Vector.Mutable as VM+-- import qualified Data.Vector.Algorithms.Merge as VA (sortBy)+-- import qualified Data.Vector.Generic as VG (convert)++import Data.Complex+import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)+import Control.Monad.Primitive++-- import Data.Sparse.Utils+-- import Data.Sparse.Types+import Data.Sparse.Internal.SVector.Mutable hiding (fromList)+import qualified Data.Sparse.Internal.SVector.Mutable as SMV (fromList)+++import Numeric.LinearAlgebra.Class++++++data SVector a = SV { svDim :: {-# UNPACK #-} !Int,+                      svIx :: V.Vector Int,+                      svVal :: V.Vector a } deriving Eq++instance Show a => Show (SVector a) where+  show (SV n ix v) = unwords ["SV (",show n,"),",show nz,"NZ:",show (V.zip ix v)]+    where nz = V.length ix++instance Functor SVector where+  fmap f (SV n ix v) = SV n ix (fmap f v)++instance Foldable SVector where+  foldr f z (SV _ _ v) = foldr f z v++instance Traversable SVector where+  traverse f (SV n ix v) = SV n ix <$> traverse f v++instance HasData (SVector a) where+  nnz = length . svIx+  -- dat (SV _ _ x) = x++++-- ** Construction ++fromDense :: V.Vector a -> SVector a+fromDense xs = SV n (V.enumFromTo 0 (n-1)) xs where+  n = V.length xs++fromVector :: Int -> V.Vector (Int, a) -> SVector a+fromVector n ixs = SV n ix xs where+  (ix, xs) = V.unzip ixs++fromList :: Int -> [(Int, a)] -> SVector a+fromList n = fromVector n . V.fromList++-- * Query++-- | O(N) : Lookup an index in a SVector (based on `find` from Data.Foldable)+index :: SVector a -> Int -> Maybe a+index cv i =+  case F.find (== i) (svIx cv) of+    Just i' -> Just $ svVal cv V.! i'+    Nothing -> Nothing+      +++-- * Lifted binary functions++-- | O(N) : Applies a function to the index _intersection_ of two CsrVector s. Useful e.g. to compute the inner product of two sparse vectors.+intersectWith :: (a -> b -> c) -> SVector a -> SVector b -> SVector c+intersectWith g v1 v2 = SV nfin ixf vf where+   nfin = V.length vf+   (ixf, vf) = V.unzip $ V.create (intersectWithM g v1 v2)++intersectWithM :: PrimMonad m =>+     (a -> b -> c)+     -> SVector a+     -> SVector b+     -> m (VM.MVector (PrimState m) (Int, c))+intersectWithM g (SV n1 ixu u) (SV n2 ixv v) = do+  vm <- VM.new n+  (vm', nfin) <- go ixu u ixv v 0 vm+  return $ VM.take nfin vm'+  where+    n = min n1 n2+    go ixu_ u_ ixv_ v_ i vm+      | V.null u_ || V.null v_ || i == n = return (vm, i)+      | otherwise =  do+           let (u, us) = headTail u_ +               (v, vs) = headTail v_ +               (ix1, ix1s) = headTail ixu_+               (ix2, ix2s) = headTail ixv_+           if ix1 == ix2 then do VM.write vm i (ix1, g u v)+                                 go ix1s us ix2s vs (i + 1) vm+                     else if ix1 < ix2 then go ix1s us ixv  v_ i vm+                                       else go ixu  u_ ix2s vs i vm++      +-- | O(N) : Applies a function to the index _union_ of two CsrVector s. Useful e.g. to compute the vector sum of two sparse vectors.+unionWith :: (t -> t -> a) -> t -> SVector t -> SVector t -> SVector a+unionWith g z v1 v2 = SV n ixf vf where+   n = max (svDim v1) (svDim v2)+   (ixf, vf) = V.unzip $ V.create (unionWithM g z v1 v2)++unionWithM :: PrimMonad m =>+     (a -> a -> b)+     -> a+     -> SVector a+     -> SVector a+     -> m (VM.MVector (PrimState m) (Int, b))+unionWithM g z (SV n1 ixu u) (SV n2 ixv v) = do+  vm <- VM.new n+  (vm', nfin) <- go ixu u ixv v 0 vm+  let vm'' = VM.take nfin vm'+  return vm''+  where+    n = min n1 n2+    go iu u_ iv v_ i vm+          | (V.null u_ && V.null v_) || i == n = return (vm , i)+          | V.null u_ = do+                 VM.write vm i (V.head iv, g z (V.head v_))+                 go iu u_ (V.tail iv) (V.tail v_) (i + 1) vm+          | V.null v_ = do+                 VM.write vm i (V.head iu, g (V.head u_) z)+                 go (V.tail iu) (V.tail u_) iv v_ (i + 1) vm+          | otherwise =  do+           let (u, us) = headTail u_+               (v, vs) = headTail v_+               (iu1, ius) = headTail iu+               (iv1, ivs) = headTail iv+           if iu1 == iv1 then do VM.write vm i (iu1, g u v)+                                 go ius us ivs vs (i + 1) vm+                         else if iu1 < iv1 then do VM.write vm i (iu1, g u z)+                                                   go ius us iv v_ (i + 1) vm+                                           else do VM.write vm i (iv1, g z v)+                                                   go iu u_ ivs vs (i + 1) vm+  +++#define SVType(t) \+  instance AdditiveGroup (SVector t) where {zeroV = SV 0 V.empty V.empty; (^+^) = unionWith (+) 0 ; negateV = fmap negate};\+  instance VectorSpace (SVector t) where {type Scalar (SVector t) = (t); n .* x = fmap (* n) x };\+  instance InnerSpace (SVector t) where {a <.> b = sum $ intersectWith (*) a b}+  -- instance Normed (CsrVector t) where {norm p v = norm' p v}+  -- instance Hilbert (CsrVector t) where {x `dot` y = sum $ intersectWithCV (*) x y };\+  ++#define SVTypeC(t) \+  instance AdditiveGroup (SVector (Complex t)) where {zeroV = SV 0 V.empty V.empty; (^+^) = unionWith (+) (0 :+ 0) ; negateV = fmap negate};\+  instance VectorSpace (SVector (Complex t)) where {type Scalar (SVector (Complex t)) = (Complex t); n .* x = fmap (* n) x };\+  instance InnerSpace (SVector (Complex t)) where {x <.> y = sum $ intersectWith (*) (conjugate <$> x) y};\++-- #define NormedType(t) \+--   instance Normed (CsrVector t) where { norm p v | p==1 = norm1 v | otherwise = norm2 v ; normalize p v = v ./ norm p v};\+--   instance Normed (CsrVector (Complex t)) where { norm p v | p==1 = norm1 v | otherwise = norm2 v ; normalize p v = v ./ norm p v}+++-- -- CVType(Int)+-- -- CVType(Integer)+-- SVType(Float)+-- SVType(Double)+-- -- CVType(CSChar)+-- -- CVType(CInt)+-- -- CVType(CShort)+-- -- CVType(CLong)+-- -- CVType(CLLong)+-- -- CVType(CIntMax)+-- SVType(CFloat)+-- SVType(CDouble)++-- SVTypeC(Float)+-- SVTypeC(Double)  +-- SVTypeC(CFloat)+-- SVTypeC(CDouble)++-- -- NormedType(Float)+-- -- NormedType(Double)+-- -- NormedType(CFloat)+-- -- NormedType(CDouble)    +++++++++{-| Modify the mutable vector operand by applying a binary function over the index union of the two.++e.g.++g = (+)+z = 0+u = [(1, a), (2, b)]+v = [(0, d), (2, e)]++unionWithSMV g z u v = [(0, d), (1, a), (2, b + e)]++invariants :++* uu, vv nonzero values+* ixu, ixv nonzero indices+* ixu, ixv sorted in ascending order+* n1 == n2+* length ixu == length uu+* length ixv == length vv++-}+unionWithSMV :: PrimMonad m =>+     (a -> a -> a)+     -> a -> SVector a -> SMVector m a -> m (SMVector m a)+unionWithSMV g z (SV n ixu uu) (SMV n2 ixm_ vm_) = do+  ixmnew <- VM.new nnzero  -- create new mutable vectors+  vmnew <- VM.new nnzero+  unless (n == n2) (error "unionWithSMV : operand vectors must have the same length")+  (ixm, vm, nfin) <- go 0 ixmnew vmnew+  let ixm' = VM.take nfin ixm+      vm' = VM.take nfin vm+  return $ SMV n ixm' vm'+  where+    nnzero = lu + lv+    lu = V.length ixu+    lv = VM.length ixm_+    go i ixm vm+      | i == lu && i == lv || i == n = return (ixm, vm, i)+      | i == lu = do+          v0 <- VM.read vm i+          VM.write ixm i i+          VM.write vm i (g z v0)+          go (i + 1) ixm vm+      | i == lv = do+          let u0 = uu V.! i+          VM.write ixm_ i i+          VM.write vm_ i (g u0 z)+          go (i + 1) ixm vm+      | otherwise = do+          let u = uu V.! i    -- read head elements and indices+              iu = ixu V.! i+          v <- VM.read vm i+          iv <- VM.read ixm i+          if iu == iv then+            do+              VM.write ixm i iu         -- write `iu` at position `i` +              VM.write vm i (g u v)+              go (i + 1) ixm vm+          else if iu < iv then+            do+              VM.write ixm i iu         -- write `iu` at position `i` +              VM.write vm i (g u z)+              go (i + 1) ixm vm+              else+                do+                  VM.write ixm i iv     -- write `iv` at position `i` +                  VM.write vm i (g z v)+                  go (i + 1) ixm vm+                               +                             ++-- -- test data+-- testUnionWithSMV :: IO (SVector Double)+-- testUnionWithSMV = do +--   let v = fromList 4 [(1, 1), (2, 1)]+--   vm <- SMV.fromList 4 [(0, pi), (2, pi)]+--   (vmres, nfin) <- unionWithSMV (+) 0 v vm+--   liftIO $ print nfin+--   freeze vmres+++++            +          ++-- * To/from SMVector++thaw :: PrimMonad m => SVector a -> m (SMVector m a)+thaw (SV n ix v) = do+  vm <- V.thaw v+  ixm <- V.thaw ix+  return $ SMV n ixm vm++freeze :: PrimMonad m => SMVector m a -> m (SVector a)+freeze (SMV n ixm vm) = do+  v <- V.freeze vm+  ix <- V.freeze ixm+  return $ SV n ix v++++++-- * helpers++both :: Arrow arr => arr b c -> arr (b, b) (c, c)+both f = f *** f++headTail :: V.Vector a -> (a, V.Vector a)+headTail = V.head &&& V.tail

@@ -0,0 +1,119 @@+{-# language FlexibleContexts, TypeFamilies #-}+module Data.Sparse.Internal.SVector.Mutable where++import Data.Foldable+import qualified Data.Vector as V +import qualified Data.Vector.Mutable as VM+import qualified Data.Vector.Generic as VG++import Control.Arrow ((&&&))+import Control.Monad.Primitive++data SMVector m a = SMV { smvDim :: {-# UNPACK #-} !Int,+                          smvIx :: VM.MVector (PrimState m) Int,+                          smvVal :: VM.MVector (PrimState m) a }+++fromList :: (PrimMonad m) => Int -> [(Int, a)] -> m (SMVector m a)+fromList n ixv = do+  let ns = length ixv+  vm <- VM.new ns+  ixm <- VM.new ns+  fromListOverwrite ixm vm n ixv+  ++fromListOverwrite :: PrimMonad m =>+     VM.MVector (PrimState m) Int+     -> VM.MVector (PrimState m) a+     -> Int+     -> [(Int, a)]+     -> m (SMVector m a)+fromListOverwrite ixm vm n ixv = go ixm vm ixv 0+    where+      go ixm_ vm_ [] _ = return $ SMV n ixm_ vm_+      go ixm_ vm_ ((i, x) : xs) iwrite = do+        VM.write vm_ iwrite x+        VM.write ixm_ iwrite i+        go ixm_ vm_ xs (iwrite + 1)+      +++-- instance Foldable (SMVector m) where+--   foldr f z (SMV _ _ v) = foldr f z v++-- instance Traversable (SMVector m) where+--   traverse f (SMV n ix v) = SMV n ix <$> traverse f v++++-- | traverse the sparse mutable vector and operate on the union set between the indices of the immutable vector `v` and those of the mutable one `vm`, _overwriting_ the values in vm.+-- Invariant: the index set of `v` is a strict subset of that of `vm` (i.e. we assume that `vm` is preallocated properly, or, we assume there won't be any out-of-bound writes attempted)+unionWithM_prealloc g z v vm@(SMV nvm vmix vmv) = undefined+++++++++-- -- | unionWithSMV takes the union of two sparse mutable vectors given a binary function and a neutral element.++-- -- unionWithSMV :: PrimMonad m =>+-- --      (a -> a -> b)+-- --      -> a+-- --      -> SMVector m a+-- --      -> SMVector m a+-- --      -> m (VM.MVector (PrimState m) (Int, b))++-- unionWithSMV g z (SMV n1 ixu uu) (SMV n2 ixv vv) = do+--   vm0 <- VM.new n+--   (vm, nfin) <- go ixu uu ixv vv 0 vm0  -- populate+--   -- let (ixm, vm') = VG.unzip (VM.take nfin vm')+--   let vm_trim = VM.take nfin vm+--       (ixm, vm') = VG.unzip $ VG.convert vm_trim+--   vOut <- V.freeze ixm+--   return (vOut, vm')+--   where+--     headTail = V.head &&& V.tail+--     n = min n1 n2+--     go iu u_ iv v_ i vm+--           | (VM.null u_ && VM.null v_) || i == n = return (vm , i)+--           | VM.null u_ = do+--               v0 <- VM.read v_ 0+--               VM.write vm i (V.head iv, g z v0)+--               go iu u_ (V.tail iv) (VM.tail v_) (i + 1) vm+--           | VM.null v_ = do+--               u0 <- VM.read u_ 0+--               VM.write vm i (V.head iu, g u0 z)+--               go (V.tail iu) (VM.tail u_) iv v_ (i + 1) vm+--           | otherwise =  do+--              u <- VM.read u_ 0+--              v <- VM.read v_ 0+--              let us = VM.tail u_+--                  vs = VM.tail v_+--              let (iu1, ius) = headTail iu+--                  (iv1, ivs) = headTail iv+--              if iu1 == iv1 then do VM.write vm i (iu1, g u v)+--                                    go ius us ivs vs (i + 1) vm+--                            else if iu1 < iv1 then do VM.write vm i (iu1, g u z)+--                                                      go ius us iv v_ (i + 1) vm+--                                              else do VM.write vm i (iv1, g z v)+--                                                      go iu u_ ivs vs (i + 1) vm+++++++++++-- -- * helpers++-- -- both :: Arrow arr => arr b c -> arr (b, b) (c, c)+-- -- both f = f *** f+++

@@ -0,0 +1,58 @@+{-# language ExistentialQuantification #-}+module Data.Sparse.Internal.Stream where++import Data.Vector.Fusion.Stream.Monadic (Stream(..), Step(..))+import Data.Vector.Fusion.Bundle.Size (Size(..))++-- from Data.Vector.Fusion.Stream.Monadic++-- -- | Result of taking a single step in a stream+-- data Step s a = Yield a s  -- ^ a new element and a new seed+--               | Skip    s  -- ^ just a new seed+--               | Done       -- ^ end of stream++-- -- | Result of taking a single step in a stream+-- data Step s a where+--   Yield :: a -> s -> Step s a+--   Skip  :: s -> Step s a+--   Done  :: Step s a++-- instance Functor (Step s) where+--   {-# INLINE fmap #-}+--   fmap f (Yield x s) = Yield (f x) s+--   fmap _ (Skip s) = Skip s+--   fmap _ Done = Done++-- -- | Monadic streams+-- data Stream m a = forall s. Stream (s -> m (Step s a)) s+++-- . -- ++data MergeState sa sb i a+  = MergeL sa sb i a+  | MergeR sa sb i a+  | MergeLeftEnded sb+  | MergeRightEnded sa+  | MergeStart sa sb++mergeStreamsWith0 f (Stream stepa sa0) (Stream stepb sb0) =+  Stream step (MergeStart sa0 sb0) -- (toMax na + toMax nb)+  where+    step (MergeStart sa sb) = do+      r <- stepa sa+      return $ case r of+        Yield (i, a) sa' -> Skip (MergeL sa' sb i a)+        Skip sa'         -> Skip (MergeStart sa' sb)+        Done             -> Skip (MergeLeftEnded sb)+    step (MergeL sa sb i a) = do+      r <- stepb sb+      return $ case r of+        Yield (j, b) sb' -> case compare i j of+          LT -> Yield (i, a)     (MergeR sa sb' j b)+          EQ -> case f a b of+            Just c  -> Yield (i, c) (MergeStart sa sb')+            Nothing -> Skip (MergeStart sa sb')+          GT -> Yield (j, b)     (MergeL sa sb' i a)+        Skip sb' -> Skip (MergeL sa sb' i a)+        Done     -> Yield (i, a) (MergeRightEnded sa)          

@@ -0,0 +1,282 @@+{-# LANGUAGE FlexibleContexts #-}+{-# language TypeFamilies, FlexibleInstances, DeriveFunctor #-}+{-# language DeriveFoldable, DeriveTraversable #-}+module Data.Sparse.Internal.TriMatrix where++-- import qualified Data.Map.Strict as M+import qualified Data.IntMap.Strict as IM+import Data.IntMap.Strict ((!))+-- import qualified Data.Set as S+-- import qualified Data.Vector as V++import Data.Foldable (foldrM)+import Data.List (sort)+import Data.Maybe (fromMaybe)+-- import Data.Monoid+import Data.Complex+import qualified Data.Graph as G+import qualified Data.Tree as T++import Numeric.Eps+import Data.Sparse.Types+-- import Data.Sparse.Utils+import Data.Sparse.Internal.SList++import Numeric.LinearAlgebra.Class+import Data.Sparse.Internal.CSC+import Data.Sparse.Internal.SVector+import qualified Data.Sparse.Internal.SVector.Mutable as SMV+import Data.Sparse.SpMatrix (fromListSM, fromListDenseSM, insertSpMatrix, zeroSM, transposeSM, sparsifySM)+import Data.Sparse.Common (prd, prd0, (@@!), nrows, ncols, lookupSM, extractRow, extractCol, SpVector, SpMatrix, foldlWithKeySV, (##), (#~#))+import Control.Iterative (IterationConfig(IterConf), modifyUntilM, modifyUntilM')+import Data.Sparse.PPrint++import Control.Monad.Catch (MonadThrow, throwM)+import Control.Exception.Common++import Control.Monad (when)+import Control.Monad.IO.Class+import Control.Monad.Trans.State -- (execStateT, get, put, modify)+import Control.Monad.Primitive+-- import Control.Monad.State (MonadState())++import qualified Data.Vector as V (Vector, freeze, thaw, fromList, toList, zip)+import qualified Data.Vector.Mutable as VM (MVector, new, set, write, read, clone)++++flattenForest :: T.Forest a -> [a]+flattenForest = concatMap T.flatten++-- | Given a lower triangular system L x = b, finds the nonzero entries of the solution vector x as the set of reachable nodes from the r.h.s. via the graph G(L^T). Node indices are _sorted_ afterwards, for a total complexity of O(N)+reachableFromRHS :: G.Graph -> V.Vector Int -> V.Vector Int+reachableFromRHS g vs = V.fromList . sort . flattenForest $ G.dfs (G.transposeG g) (V.toList vs)++++-- triLowerSolve ll b = do+--   initializeSoln xinz b+--   where+--   -- xinz : nonzeros of solution vector x obtained from reachable nodes of b via G(L^T)+--   xinz = reachableFromRHS (cscToGraph ll) (svIx b)++-- -- tlUpdateColumn lldiag llsubdiag x j = undefined where+-- --   xj' = xj / lldiag++-- | Build the initial solution vector `x` for the triangular solve, given a vector of nonzero indices `ixnz` and the right hand side of the linear system `bb`:+-- -- 1. Initialize an empty mutable vector of length `length ixnz` to 0+-- -- 2. Copy the entries from `bb` to `x`.+--+-- Note: this assumes that the index set of `bb` is strictly contained within that of `x` (which is true for the case of the triangular solve, see `reachableFromRHS`)+initializeSoln ::+  (PrimMonad m, Num a) => V.Vector Int -> SVector a -> m (SMV.SMVector m a)+initializeSoln ixnz (SV n ixb b) = do+  let nnzx = length ixnz +  xm <- VM.new nnzx+  ixnzm <- V.thaw ixnz+  VM.set xm 0+  SMV.fromListOverwrite ixnzm xm n $ V.toList (V.zip ixb b)+++++tlUpdateSubdiag :: VectorSpace v => v -> Scalar v -> v -> v+tlUpdateSubdiag lsubdiag xj x = x ^-^ (xj .* lsubdiag)+   +++  ++++g0 = G.buildG (0,2) [(0,0),(2,0),(1,1),(2,2)]++t0 :: T.Forest Int+t0 = G.dfs (G.transposeG g0) [0]++g1 = G.buildG (0, 4) [(0,0), (1,0), (1,1), (2,0), (2,2), (3,1), (3,2), (3,3), (4,0), (4,4)]++++{- | triangular sparse matrix, row-major order++Intmap-of-sparse lists+* fast random access of rows+* fast consing of row elements+-}++newtype TriMatrix a = TM { unTM :: IM.IntMap (SList a) } deriving (Show, Functor)++emptyIMSL :: Int -> IM.IntMap (SList a)+emptyIMSL n = IM.fromList [(i, emptySL) | i <- [0 .. n-1]]++emptyTM :: Int -> TriMatrix a+emptyTM n = TM (emptyIMSL n)++-- | `appendIM i x im` appends an element `x` to the i'th SList in an IntMap-of-SLists structure+appendIM :: IM.Key -> (Int, a) -> IM.IntMap (SList a) -> IM.IntMap (SList a)+appendIM i x im = IM.insert i (x `consSL` e) im where+  e = fromMaybe emptySL (IM.lookup i im)+++-- | Nested lookup with default value = 0+lookupWD :: Num a =>+     (irow -> mat -> Maybe row)    -- ^ row lookup+     -> (jcol -> row -> Maybe a)   -- ^ in-row lookup+     -> mat                +     -> irow+     -> jcol+     -> a+lookupWD rlu clu aa i j = fromMaybe 0 (rlu i aa >>= clu j)++ ++++++++{- | LU factorization : store L and U^T in TriMatrix format -}+++lu :: (Scalar (SpVector t) ~ t, Elt t, AdditiveGroup t, VectorSpace (SpVector t),+      MonadThrow m, MonadIO m, PrintDense (SpMatrix t),+      Epsilon t) =>+     SpMatrix t -> m (SpMatrix t, SpMatrix t) -- ^ L, U+lu amat = do+  let d@(m,n) = (nrows amat, ncols amat)+      q (_, _, i) = i == m    -- stopping criterion+      luInit = (lmat0, umat0, 1) where+         urow0 = extractRow amat 0                 -- first row of U+         lcol0 = extractCol amat 0 ./ (urow0 @@ 0) -- first col of L, div by U00+         umat0 = foldlWithKeySV ins (emptyIMSL n) urow0 -- populate umat0+         lmat0 = IM.insert 0 (SL [(0, 1)]) l0 where     -- populate lmat0+           l0 = foldlWithKeySV ins (emptyIMSL m) lcol0 +         ins acc i x = appendIM i (0, x) acc+      luStep (lmat, umat, i) = do+          let (umat', uii) = uStep amat lmat umat i  -- new U+          when (nearZero uii) $+             throwM (NeedsPivoting "LU" (unwords ["U", show (i,i)]) :: MatrixException Double)+          let lmat' = lStep amat lmat umat' uii i  -- new L+          return (lmat', umat', i + 1)             +  -- (lfin, ufin, _) <- execStateT (modifyUntilM q luStep) luInit+  (lfin, ufin, _) <- modifyUntilM' (luConfig d) q luStep luInit  +  let uu = fillSM d True ufin+      ll = fillSM d False lfin+  return (ll, uu)+++luConfig d = IterConf 0 True vf prf where+        vf (l, u, _) = (l, u)+        prf (l, u) = do+          prd0 $ fillSM d False l+          prd0 $ fillSM d True u+++++uStep :: (Elt a, Epsilon a, AdditiveGroup a) =>+     SpMatrix a+     -> IM.IntMap (SList a)+     -> IM.IntMap (SList a)+     -> IM.Key+     -> (IM.IntMap (SList a), a)   -- ^ updated U, i'th diagonal element Uii+uStep amat lmat umat i = (foldr ins umat [i .. n-1], udiag) where+  n = ncols amat+  udiag = amat@@!(i,i) - (li <.> (umat ! i)) -- i'th diag element of U+  li = lmat ! i                            -- i'th row of L+  ins j acc+      | i == j   = appendIM j (i, udiag) acc+      | isNz uij = appendIM j (i, uij) acc+      | otherwise = acc where+    uij = aij - li <.> uj +    aij = amat @@! (i,j)+    uj = umat ! j+  ++lStep :: (Elt a, Epsilon a, AdditiveGroup a) =>+     SpMatrix a+     -> IM.IntMap (SList a)+     -> IM.IntMap (SList a)+     -> a                   -- ^ diagonal element of U (must be nonzero)+     -> IM.Key+     -> IM.IntMap (SList a) -- ^ updated L+lStep amat lmat umat udiag j = foldr ins lmat [j .. m-1] where+  m = nrows amat+  uj = umat ! j+  ins i acc+    | i == j   = appendIM i (j, 1) acc  -- write 1 on the diagonal +    | isNz lij = appendIM i (j, lij) acc+    | otherwise = acc where+    lij = (aij - li <.> uj)/udiag+    aij = amat @@! (i,j)+    li = lmat ! i+++++fillSM :: (Rows, Cols) -> Bool -> IM.IntMap (SList a) -> SpMatrix a+fillSM (m,n) transpq tm = IM.foldlWithKey rowIns (zeroSM m n) tm where+  rowIns accRow i row = foldr ins accRow (unSL row) where+    ins (j, x) acc | transpq = insertSpMatrix j i x acc   -- transposed fill+                   | otherwise = insertSpMatrix i j x acc++++-- test data++test mm = do+  (l, u) <- lu mm+  prd l+  prd u+  prd mm+  prd $ l #~# u++++tm0, tm2, tm4, tm9 :: SpMatrix Double++tm0 = fromListDenseSM 2 [1,3,2,4]++tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41]+tm4 = sparsifySM $ fromListDenseSM 4 [1,0,0,0,2,5,0,10,3,6,8,11,4,7,9,12]+tm9 = fromListSM (4, 4) [(0,0,pi), (1,1, 3), (3, 0, 23), (1,3, 45), (2,2,4), (3,2, 1), (3,1, 5), (3,3, exp 1)]++-- -- complex+tmc4 :: SpMatrix (Complex Double)+tmc4 = fromListDenseSM 3 [3:+1, 4:+(-1), (-5):+3, 2:+2, 3:+(-2), 5:+0.2, 7:+(-2), 9:+(-1), 2:+3]+++++-- λ> test tmc4++-- 1.00            , _               , _               +-- 1.10 - 0.70i    , 1.00            , _               +-- -1.20 + 1.40i   , -0.56 + 1.04i   , 1.00            ++-- 3.00 + 1.00i    , 2.00 + 2.00i    , 7.00 - 2.00i    +-- _               , 2.20 - 5.60i    , -0.10 - 3.70i   +-- _               , _               , 16.99 + 8.24i   ++++++++-- playground++data T x+  = Leaf+  | Branch (T x) x (T x)+  deriving (Functor, Foldable, Traversable, Show)++next :: Monad m => StateT Int m Int+next = do x <- get ; modify (+ 1) ; return x++next' f = do x <- get;  modify f; return x++labelElt :: Monad m => x -> StateT Int m (Int, x)+labelElt x = (,) <$> next <*> pure x

@@ -24,10 +24,10 @@  -- | Given a number of rows(resp. columns) `n` and a _sorted_ Vector of Integers in increasing order (containing the row(col) indices of nonzero entries), return the cumulative vector of nonzero entries of length `n + 1` (the "row(col) pointer" of the CSR(CSC) format). NB: Fused count-and-accumulate -- E.g.:--- > csrPtrV (==) 4 (V.fromList [1,1,2,3])+-- > csPtrV (==) 4 (V.fromList [1,1,2,3]) -- [0,0,2,3,4]-csrPtrV :: (a -> Int -> Bool) -> Int -> V.Vector a -> V.Vector Int-csrPtrV eqf n xs = V.create createf where+csPtrV :: (a -> Int -> Bool) -> Int -> V.Vector a -> V.Vector Int+csPtrV eqf n xs = V.create createf where   createf :: ST s (VM.MVector s Int)   createf = do     let c = 0

@@ -1,3 +1,4 @@+{-# language FlexibleContexts, RecordWildCards #-} ----------------------------------------------------------------------------- -- | -- Copyright   :  (C) 2016 Marco Zocca@@ -6,8 +7,14 @@ -- Stability   :  provisional -- Portability :  portable --+--+-- Pretty printing helper functions ------------------------------------------------------------------------------module Data.Sparse.PPrint (prd, prd0, PrintDense, newline, PPrintOptions, prdef, prepD, prepC, printDN, printCN) where+module Data.Sparse.PPrint (prd, prd0, PrintDense, newline+                           , PPrintOptions(..)+                             -- , prdefR, prdefC+                           , printDN , printCN+                          ) where  import Data.Complex import Data.List@@ -27,92 +34,109 @@ newline = putStrLn ""    --- | Pretty printing options: # of integer digits, # of decimal digits+-- | Pretty printing options: total length in # digits (including the decimal point), # of decimal digits data PPrintOptions =   PPOpts {-     pprintInt :: Int,-     pprintDec :: Int+     pprintLen, pprintDec, pprintColWidth :: Int      } deriving (Eq, Show) -prdef :: PPrintOptions-prdef = PPOpts 1 2+-- | Some defaults+prdefR, prdefC :: PPrintOptions+prdefR = PPOpts 4 2 7   -- reals+prdefC = PPOpts 4 2 16  -- complex values  +-- | Pretty print an array of real numbers+printDN :: (PrintfArg a, Epsilon a, Ord a) =>+     Int -> Int -> PPrintOptions -> [a] -> String+printDN l n = printNpad l n f where+  f o x | isNz x = printf (prepD o x) x+        | otherwise = printf "_"  +-- | Pretty print an array of complex numbers+printCN :: (PrintfArg a, Epsilon a, Epsilon (Complex a), Ord a) =>+     Int -> Int -> PPrintOptions -> [Complex a] -> String+printCN l n = printNpad l n f where+  f o x | nearZero (re x) && isNz (imagPart x) =+               printf (prepD o (imagPart x)++ "i") (aim x)+        | nearZero (imagPart x) && isNz (re x) =+               printf (prepD o (realPart x)) (re x)+        | isNz x = printf (prepC o x) (re x) (aim x)+        | otherwise = printf "_"+++++++-- | printf an array of items with padding space to render a fixed column width+printNpad ::+     Int     -- ^ Length of list to be provided+     -> Int  -- ^ +     -> (PPrintOptions -> a -> String)+     -> PPrintOptions -> [a] -> String+printNpad llen nmax f o@PPOpts{..} xxl = commas [h,l] where+  h = commas $ take hlen ll+  l = last ll+  hlen = min (llen-1) (nmax-1)+  ll = unfoldr g (0, xxl) +  g (i, x:xs) | i<nmax-2 || llen>=nmax-1 = Just (s', sxs)+              | i==nmax-2 = Just (dots', sxs)+              | null xs = Just (s', sxs)+              | otherwise = Just ("", sxs) where                  +                  s = f o x+                  sxs = (succ i, xs)+                  s' = s ++ spaces (n - length s) +                  dots' = dots ++ spaces (n - length dots)+  g (_, []) = Nothing+  n = pprintColWidth+  dots = " ... "++++++ -- | printf format string--- prepD :: (Ord t, Epsilon t) => PPrintOptions -> t -> String-prepD opts x = pstr -- printf pstr x-  where-  pstr | nearZero x = "_"-       | abs x >= 10 || abs x < 0.1 = s ++ "e"-       | otherwise = s ++ "f"-    where-      s = concat ["%" , show ni, ".", show nd]-  nd = pprintDec opts-  ni = pprintInt opts+prepD :: (Ord t, Epsilon t) => PPrintOptions -> t -> String+prepD PPOpts{..} x = s where+  s | nearZero x  = "_"+    | abs x > magHi || abs x < magLo = s0 ++ "e"+    | otherwise = s0 ++ "f"+  s0 = concat ["%" , show nl, ".", show nd]+  -- s0 = "%1." ++ show nd+  nl = pprintLen  +  nd = pprintDec +  nint = nl - nd  -- # of integer digits+  magLo = 10 ** (- fromIntegral nd)+  magHi = 10 ** fromIntegral nint+             + -- | printf format string for a Complex prepC :: (Epsilon t, Ord t) => PPrintOptions -> Complex t -> String prepC opts (r :+ i) = prepD opts r ++ oi where-    oi | isNz i = concat [s, prepD opts i', "i"]-       | otherwise = []-    s | signum i >= 0 = "+"-      | otherwise = "-"-    i' = abs i+  oi = concat [s, prepD opts i', "i"]+  s | signum i >= 0 = " + "+    | otherwise = " - "+  i' = abs i  --- | printf for a list of values------ > printDN prepD (PPOpts 1 3) [1,pi]-printDN-  :: (PrintfArg t1, PrintfType t, Ord t1, Epsilon t1) =>-     PPrintOptions -> [t1] -> t-printDN opts xl-  | null xl = printf "\n"-  | null xs = printf (prepD opts x) x-  | n==1 = let [x1]=xs in printf (commas (pr <$> xl)++"\n") x x1-  | n==2 = let [x1,x2]=xs in printf (commas (pr <$> xl)++"\n") x x1 x2-  | n==3 = let [x1,x2,x3]=xs in printf (commas (pr <$> xl)++"\n") x x1 x2 x3-  | n==4 = let [x1,x2,x3,x4]=xs in printf (commas (pr <$> xl)++"\n") x x1 x2 x3 x4-  | otherwise = let xs@[x,x1,x2,x3]=take 4 xl-                    xfin=last xl-                in printf (commas (pr <$> xs) ++ ", ... , " ++ pr xfin++"\n") x x1 x2 x3 xfin-  where-    (x:xs) = xl-    n = length xs-    pr = prepD opts  --- | printf for list of complex values--- --- > printCN pdef [(1:+pi), (3.5:+4.3), (pi:+(-3.4))]-printCN-  :: (PrintfArg t1, PrintfType t, Epsilon t1, Ord t1) =>-     PPrintOptions -> [Complex t1] -> t-printCN opts xl-  | null xl = printf "\n"-  | null xs = printf (pr x) (re x) (aim x)-  | n==1 = let [x1]=xs-           in printf (commas (pr <$> xl)++"\n") (re x) (aim x) (re x1) (aim x1)-  | n==2 = let [x1,x2]=xs-           in printf (commas (pr <$> xl) ++"\n") (re x) (aim x) (re x1) (aim x1) (re x2) (aim x2)-  | n==3 = let [x1,x2,x3]=xs-           in printf (commas (pr <$> xl)++"\n") (re x) (aim x) (re x1) (aim x1) (re x2) (aim x2) (re x3) (aim x3)-  | n==4 = let [x1,x2,x3,x4]=xs-           in printf (commas (pr <$> xl)++"\n") (re x) (aim x) (re x1) (aim x1) (re x2) (aim x2) (re x3) (aim x3) (re x4) (aim x4)-  | otherwise = let xs@[x,x1,x2,x3]=take 4 xl-                    xfin = last xl-                in printf (commas (pr <$> xs) ++ ", ... , " ++ pr xfin++"\n") (re x) (aim x) (re x1) (aim x1) (re x2) (aim x2) (re x3) (aim x3) (re xfin) (aim xfin)-  where-    (x:xs) = xl-    pr = prepC opts-    n = length xs+     +-- | Helpers++spaces :: Int -> String+spaces n = replicate n ' '+ commas :: [String] -> String     commas = intercalate ", "  re :: Complex a -> a re = realPart+ aim :: Num a => Complex a -> a aim = abs . imagPart

@@ -11,8 +11,8 @@ ----------------------------------------------------------------------------- module Data.Sparse.SpMatrix where -import Control.Exception.Common-import Data.Sparse.SpVector+-- import Control.Exception.Common+-- import Data.Sparse.SpVector import Data.Sparse.Utils import Data.Sparse.Types @@ -23,8 +23,8 @@ import qualified Data.Sparse.Internal.IntM as I import Data.Sparse.Internal.IntMap2 - import GHC.Exts+import Text.Printf  import qualified Data.IntMap.Strict as IM @@ -32,7 +32,7 @@ import Data.Foldable (foldl') import Data.Maybe -import Data.VectorSpace hiding (magnitude)+-- import Data.VectorSpace hiding (magnitude)   -- *@@ -56,16 +56,18 @@                 deriving (Eq, Functor, Foldable)  --sizeStr :: (FDSize f ~ (a1, a2), Sparse f a, Show a2, Show a1) => f a -> String-sizeStr sm =-  unwords ["(",show nr,"rows,",show nc,"columns ) ,",show nz,"NZ ( sparsity",show sy,")"] where-  (nr, nc) = dim sm+-- sizeStr :: (FDSize f ~ (a1, a2), Sparse f a, Show a2, Show a1) => f a -> String+sizeStrSM :: SpMatrix a -> String+sizeStrSM sm@(SM (nr, nc) _) =+  unwords ["(",show nr,"rows,",show nc,"columns ) ,",show nz,"NZ ( density",sys,")"] where+  -- (nr, nc) = dim sm   nz = nnz sm   sy = spy sm :: Double+  sys = printf "%1.3f %%" (sy * 100) :: String + instance Show a => Show (SpMatrix a) where-  show sm@(SM _ x) = unwords ["SM",sizeStr sm,show (toList $ toList <$> x)]+  show sm@(SM _ x) = unwords ["SM",sizeStrSM sm,show (toList $ toList <$> x)]   -- show sm@(SM _ x) = show x  -- instance Functor SpMatrix where@@ -90,21 +92,22 @@   -- | 'SpMatrix'es are maps between finite-dimensional spaces-instance FiniteDim SpMatrix where-  type FDSize SpMatrix = (Rows, Cols)+instance FiniteDim (SpMatrix a) where+  type FDSize (SpMatrix a) = (Rows, Cols)   dim = smDim -instance HasData SpMatrix a where-  type HDData SpMatrix a = IntM (IntM a)+instance HasData (SpMatrix a) where+  type HDData (SpMatrix a) = IntM (IntM a)   nnz = nzSM   dat = smData -instance Sparse SpMatrix a where+instance Sparse (SpMatrix a) where   spy = spySM  -- | 'SpMatrix'es are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0)-instance Num a => SpContainer SpMatrix a where-  type ScIx SpMatrix = (Rows, Cols)+instance Num a => SpContainer (SpMatrix a) where+  type ScIx (SpMatrix a) = (Rows, Cols)+  type ScElem (SpMatrix a) = a   scInsert (i,j) = insertSpMatrix i j   scLookup m (i, j) = lookupSM m i j   m @@ d | isValidIxSM m d = m @@! d@@ -126,6 +129,7 @@   -- *** Diagonal matrix+-- | `mkDiagonal n ll` : create a diagonal matrix of size `n` from a list `ll` of elements mkDiagonal :: Int -> [a] -> SpMatrix a mkDiagonal n = mkSubDiagonal n 0 @@ -140,22 +144,34 @@ -- | Permutation matrix from a (possibly incomplete) list of row swaps starting from row 0 -- e.g. `permutationSM 5 [1,3]` first swaps rows (0, 1) and then rows (1, 3) : -- --- [0,1,0,0,0]--- [0,0,0,1,0]--- [0,0,1,0,0]--- [1,0,0,0,0]--- [0,0,0,0,1]+-- >>> prd (permutationSM 5 [1,3] :: SpMatrix Double)+-- +-- @+-- ( 5 rows, 5 columns ) , 5 NZ ( density 20.000 % )+-- +-- _      , 1.00   , _      , _      , _      +-- _      , _      , _      , 1.00   , _      +-- _      , _      , 1.00   , _      , _      +-- 1.00   , _      , _      , _      , _      +-- _      , _      , _      , _      , 1.00+-- @ permutationSM :: Num a => Int -> [IxRow] -> SpMatrix a 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 and 4 :------ [1,0,0,0,0]--- [0,1,0,0,0]--- [0,0,0,0,1]--- [0,0,0,1,0]--- [0,0,1,0,0]+-- +-- >>> prd (permutPairsSM 5 [(2,4)] :: SpMatrix Double)+-- +-- @+-- ( 5 rows, 5 columns ) , 5 NZ ( density 20.000 % )+-- +-- 1.00   , _      , _      , _      , _      +-- _      , 1.00   , _      , _      , _      +-- _      , _      , _      , _      , 1.00   +-- _      , _      , _      , 1.00   , _      +-- _      , _      , 1.00   , _      , _+-- @ permutPairsSM :: Num a => Int -> [(IxRow, IxRow)] -> SpMatrix a permutPairsSM n iix = go iix (eye n) where   go ((i1, i2):iis) m = go iis (swapRows i1 i2 m)@@ -406,8 +422,8 @@ isUpperTriSM m = m == lm where   lm = ifilterSM (\i j _ -> i <= j) m --- -- -- |Is the matrix orthogonal? i.e. Q^t ## Q == I--- -- isOrthogonalSM :: (Eq a, Epsilon a) => SpMatrix a -> Bool+-- |Is the matrix orthogonal? i.e. Q^t ## Q == I+isOrthogonalSM :: (Eq a, Epsilon a, MatrixRing (SpMatrix a)) => SpMatrix a -> Bool isOrthogonalSM sm@(SM (_,n) _) = rsm == eye n where   rsm = roundZeroOneSM $ transpose sm ## sm @@ -539,7 +555,7 @@  -- ** Matrix stacking --- | Vertical stacking+-- | Vertical stacking of matrix blocks vertStackSM, (-=-) :: SpMatrix a -> SpMatrix a -> SpMatrix a vertStackSM mm1 mm2 = SM (m, n) $ I.union u1 u2 where   nro1 = nrows mm1@@ -550,7 +566,7 @@  (-=-) = vertStackSM --- | Horizontal stacking+-- | Horizontal stacking of matrix blocks horizStackSM, (-||-) :: SpMatrix a -> SpMatrix a -> SpMatrix a horizStackSM mm1 mm2 = t (t mm1 -=- t mm2) where   t = transposeSM@@ -559,7 +575,7 @@   --- | Assembles a square matrix from a list of square matrices, arranging these along the main diagonal+-- | Assemble a block-diagonal square matrix from a list of square matrices, arranging these along the main diagonal fromBlocksDiag :: [SpMatrix a] -> SpMatrix a fromBlocksDiag mml = fromListSM (n, n) lstot where   dims = map nrows mml

@@ -29,16 +29,18 @@ import Data.Complex import Data.Maybe +import Text.Printf+ import qualified Data.IntMap.Strict as IM import qualified Data.Foldable as F import qualified Data.Vector as V -import Data.VectorSpace hiding (magnitude)     + -- * Sparse Vector  data SpVector a = SV { svDim :: {-# UNPACK #-} !Int ,@@ -47,7 +49,7 @@ instance Show a => Show (SpVector a) where   show (SV d x) = "SV (" ++ show d ++ ") "++ show (toList x) --- | SpVector sparsity+-- | SpVector density spySV :: Fractional b => SpVector a -> b spySV s = fromIntegral (size (dat s)) / fromIntegral (dim s) @@ -57,7 +59,9 @@   sizeStrSV :: SpVector a -> String-sizeStrSV sv = unwords ["(",show (dim sv),"elements ) , ",show (nzSV sv),"NZ ( sparsity", show (spy sv),")"]+sizeStrSV sv = unwords ["(",show (dim sv),"elements ) , ",show (nzSV sv),"NZ ( density", sys,")"] where+  sy = spy sv :: Double+  sys = printf "%1.3f %%" (sy * 100) :: String   @@ -72,29 +76,34 @@     foldr f d v = F.foldr f d (svData v)  +foldlWithKeySV, foldlWithKeySV' :: (a -> IM.Key -> b -> a) -> a -> SpVector b -> a+foldlWithKeySV f d v = foldlWithKey f d (svData v) +foldlWithKeySV' f d v = foldlWithKey' f d (svData v)  + -- | 'SpVector's form a vector space because they can be multiplied by a scalar   -- | 'SpVector's are finite-dimensional vectors-instance FiniteDim SpVector where-  type FDSize SpVector = Int+instance FiniteDim (SpVector a) where+  type FDSize (SpVector a) = Int   dim = svDim   -instance HasData SpVector a where-  type HDData SpVector a = IntM a+instance HasData (SpVector a) where+  type HDData (SpVector a) = IntM a   dat = svData   nnz (SV _ x) = length x  -instance Sparse SpVector a where+instance Sparse (SpVector a) where   spy = spySV   -- | 'SpVector's are sparse containers too, i.e. any specific component may be missing (so it is assumed to be 0)-instance Elt a => SpContainer SpVector a where-  type ScIx SpVector = Int+instance Elt a => SpContainer (SpVector a) where+  type ScIx (SpVector a) = Int+  type ScElem (SpVector a) = a   scInsert = insertSpVector   scLookup v i = lookupSV i v   scToList = toListSV@@ -112,19 +121,19 @@   -#define SpVectorInstance(t) \-  instance AdditiveGroup (SpVector (t)) where { zeroV = SV 0 empty; (^+^) = liftU2 (+); negateV = fmap negate };\-  instance AdditiveGroup (SpVector (Complex t)) where { zeroV = SV 0 empty; (^+^) = liftU2 (+); negateV = fmap negate };\-  instance VectorSpace (SpVector t) where { type (Scalar (SpVector t)) = t; n *^ v = scale n v};\-  instance VectorSpace (SpVector (Complex t)) where { type (Scalar (SpVector (Complex t))) = Complex t; n *^ v = scale n v};\-  instance InnerSpace (SpVector (t)) where { (<.>) = dotS };\-  instance InnerSpace (SpVector (Complex (t))) where { (<.>) = dotS };\-  instance Normed (SpVector (t)) where {type RealScalar (SpVector (t)) = t; type Magnitude (SpVector (t)) = t; norm1 (SV _ v) = norm1 v; norm2Sq (SV _ v) = norm2Sq v ; normP p (SV _ v) = normP p v; normalize p (SV n v) = SV n (normalize p v); normalize2 (SV n v) = SV n (normalize2 v)};\-  instance Normed (SpVector (Complex t)) where {type RealScalar (SpVector (Complex t)) = t; type Magnitude (SpVector (Complex t)) = t; norm1 (SV _ v) = norm1 v; norm2Sq (SV _ v) = norm2Sq v ; normP p (SV _ v) = normP p v; normalize p (SV n v) = SV n (normalize p v); normalize2 (SV n v) = SV n (normalize2 v)}+-- #define SpVectorInstance(t) \+--   instance AdditiveGroup (SpVector (t)) where { zeroV = SV 0 empty; (^+^) = liftU2 (+); negateV = fmap negate };\+--   instance AdditiveGroup (SpVector (Complex t)) where { zeroV = SV 0 empty; (^+^) = liftU2 (+); negateV = fmap negate };\+--   instance VectorSpace (SpVector t) where { type (Scalar (SpVector t)) = t; n *^ v = scale n v};\+--   instance VectorSpace (SpVector (Complex t)) where { type (Scalar (SpVector (Complex t))) = Complex t; n *^ v = scale n v};\+--   instance InnerSpace (SpVector (t)) where { (<.>) = dotS };\+--   instance InnerSpace (SpVector (Complex (t))) where { (<.>) = dotS };\+--   instance Normed (SpVector (t)) where {type RealScalar (SpVector (t)) = t; type Magnitude (SpVector (t)) = t; norm1 (SV _ v) = norm1 v; norm2Sq (SV _ v) = norm2Sq v ; normP p (SV _ v) = normP p v; normalize p (SV n v) = SV n (normalize p v); normalize2 (SV n v) = SV n (normalize2 v)};\+--   instance Normed (SpVector (Complex t)) where {type RealScalar (SpVector (Complex t)) = t; type Magnitude (SpVector (Complex t)) = t; norm1 (SV _ v) = norm1 v; norm2Sq (SV _ v) = norm2Sq v ; normP p (SV _ v) = normP p v; normalize p (SV n v) = SV n (normalize p v); normalize2 (SV n v) = SV n (normalize2 v)}   -SpVectorInstance(Double)+-- SpVectorInstance(Double) -- SpVectorInstance(Float)  

@@ -1,13 +1,80 @@ module Data.Sparse.Utils where +-- import Control.Arrow (first, second)+import Data.Ord import qualified Data.Vector as V  -- * Misc. utilities  +++++-- | Intersection and union of sparse lists having indices in _ascending_ order+intersectWith :: Ord i => (a -> a -> b) -> [(i, a)] -> [(i, a)] -> [b]+intersectWith f = intersectWith0 (comparing fst) (lift2snd f)++unionWith :: Ord i =>+     (a -> a -> a) -> a -> [(i, a)] -> [(i, a)] -> [(i, a)]+unionWith = unionWith0 compare++++-- | Intersection of sparse lists having indices in _descending_ order+intersectWithD :: Ord i => (a -> a -> b) -> [(i, a)] -> [(i, a)] -> [b]+intersectWithD f = intersectWith0 (comparing (Down . fst)) (lift2snd f)++-- | Union of sparse lists having indices in _descending_ order+unionWithD :: Ord i =>+     (a -> a -> a) -> a -> [(i, a)] -> [(i, a)] -> [(i, a)]+unionWithD = unionWith0 (comparing Down)+++--+  +intersectWith0 :: (a -> b -> Ordering) -> (a -> b -> c) -> [a] -> [b] -> [c]+intersectWith0 q f = go [] where+  go acc ls@(x : xs) rs@(y : ys) =+    case q x y of EQ -> go (f x y : acc) xs ys+                  LT -> go acc xs rs+                  _  -> go acc ls ys+  go acc [] _ = acc+  go acc _ [] = acc++++++++unionWith0 :: (i -> i -> Ordering) -> (a -> a -> a) -> a -> [(i, a)] -> [(i, a)] -> [(i, a)]+unionWith0 q f z = go [] where+  go acc ls@((ix, x) : xs) rs@((iy, y) : ys) =+    case q ix iy of EQ -> go ((ix, f x y) : acc) xs ys+                    LT -> go ((ix, f x z) : acc) xs rs+                    _  -> go ((iy, f z y) : acc) ls ys+  go acc [] r = acc ++ r+  go acc l [] = acc ++ l++++++  +-- | Lift a binary function onto the second entry of a tuple+lift2snd :: (t -> t1 -> t2) -> (a, t) -> (a1, t1) -> t2+lift2snd f a b = f (snd a) (snd b)+++++++ -- | 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+safe :: (t -> Bool) -> (t -> a) -> t -> Maybe a+safe q f v | q v = Nothing               | otherwise = Just $ f v  @@ -90,30 +157,32 @@   head' :: V.Vector a -> Maybe a-head' = harness V.null V.head+head' = safe V.null V.head  tail' :: V.Vector a -> Maybe (V.Vector a)-tail' = harness V.null V.tail+tail' = safe V.null V.tail    --- | a cons-based moving-window datatype of length at least 3 -data W3 a = W3 a a a [a] deriving (Eq, Show)-initW3 :: a -> a -> a -> [a] -> W3 a-initW3 = W3 -pushW3 :: a -> W3 a -> W3 a-pushW3 i (W3 i0 i1 i2 is) = W3 i i0 i1 (i2 : initSafe is) -fstW3 (W3 i _ _ _) = i-sndW3 (W3 _ i _ _) = i-thirdW3 (W3 _ _ i _) = i -withInitW3, withTailW3 :: W3 t -> (t -> t -> t1) -> t1-withInitW3 (W3 a b _ _) f = f a b-withTailW3 (W3 _ b c _) f = f b c+-- unionWith0 :: (a -> a -> Ordering) -> (a -> a -> a) -> a -> [a] -> [a] -> [a]+-- unionWith0 q f z = go [] where+--   go acc ls@(x : xs) rs@(y : ys) =+--     case q x y of EQ -> go (f x y : acc) xs ys+--                   LT -> go (f x z : acc) xs rs+--                   _  -> go (f z y : acc) ls ys+--   go acc [] r = acc ++ r+--   go acc l [] = acc ++ l -initSafe (x:xs) = x : initSafe xs-initSafe [] = []+-- union :: Ord a => [a] -> [a] -> [a]+-- union u_ v_ = go u_ v_ where+--   go [] x = x+--   go y [] = y+--   go uu@(u:us) vv@(v:vs)+--     | u == v =    u : go us vs+--     | u < v =     u : go us vv +--     | otherwise = v : go uu vs

@@ -5,7 +5,7 @@ -- | -- Module      :  Numeric.LinearAlgebra.Class -- Copyright   :  (c) Marco Zocca 2017--- License     :  GPL-style (see the file LICENSE)+-- License     :  GPL-3 (see the file LICENSE) -- -- Maintainer  :  zocca marco gmail -- Stability   :  experimental@@ -18,8 +18,6 @@  -- import Control.Applicative import Data.Complex--- import Data.Ratio--- import Foreign.C.Types (CSChar, CInt, CShort, CLong, CLLong, CIntMax, CFloat, CDouble)  -- import Control.Exception -- import Control.Exception.Common@@ -30,8 +28,6 @@  import qualified Data.Vector as V (Vector) -import Data.VectorSpace hiding (magnitude)- import Data.Sparse.Types import Numeric.Eps @@ -48,7 +44,7 @@  instance Elt Double where {type EltMag Double = Double ; mag = id} instance Elt Float where {type EltMag Float = Float; mag = id}-instance RealFloat e => Elt (Complex e) where+instance (RealFloat e) => Elt (Complex e) where   type EltMag (Complex e) = e   conj = conjugate   mag = magnitude@@ -56,39 +52,54 @@   +infixl 6 ^+^, ^-^ +-- * Additive group+class AdditiveGroup v where+  -- | The zero element: identity for '(^+^)'+  zeroV :: v+  -- | Add vectors+  (^+^) :: v -> v -> v+  -- | Additive inverse+  negateV :: v -> v+  -- | Group subtraction+  (^-^) :: v -> v -> v  --- * Vector space---- | Scale a vector -(.*) :: VectorSpace v => Scalar v -> v -> v-(.*) = (*^)+infixr 7 .* --- | Scale a vector by the reciprocal of a number (e.g. for normalization)-(./) :: (VectorSpace v, Fractional (Scalar v)) => v -> Scalar v -> v-v ./ n = recip n .* v+-- * Vector space @v@.+class (AdditiveGroup v, Num (Scalar v)) => VectorSpace v where+  type Scalar v :: *+  -- | Scale a vector+  (.*) :: Scalar v -> v -> v --- | Convex combination of two vectors (NB: 0 <= `a` <= 1). -cvx :: (VectorSpace e, Num (Scalar e)) => Scalar e -> e -> e -> e-cvx a u v = a .* u ^+^ ((1-a) .* v)+-- | Adds inner (dot) products.+class (VectorSpace v, AdditiveGroup (Scalar v)) => InnerSpace v where+  -- | Inner/dot product+  (<.>) :: v -> v -> Scalar v +-- | Inner product+dot :: InnerSpace v => v -> v -> Scalar v+dot = (<.>)+   --- linearCombination :: (VectorSpace v , Foldable t) => t (Scalar v, v) -> v--- linearCombination  =  foldr (\(a, x) (b, y) -> (a .* x) ^+^ (b .* y)) +infixr 7 ./+infixl 7 *. --- linComb a v = a .* v+-- | Scale a vector by the reciprocal of a number (e.g. for normalization)+(./) :: (VectorSpace v, s ~ Scalar v, Fractional s) => v -> s -> v+v ./ s = (recip s) .* v +-- | Vector multiplied by scalar+(*.) :: (VectorSpace v, s ~ Scalar v) => v -> s -> v+(*.) = flip (.*)     --- * Hilbert space (inner product)---- | Inner product-dot :: InnerSpace v => v -> v -> Scalar v-dot = (<.>)-  -+-- | Convex combination of two vectors (NB: 0 <= `a` <= 1). +cvx :: VectorSpace v => Scalar v -> v -> v -> v+cvx a u v = a .* u ^+^ ((1-a) .* v)   @@ -152,27 +163,27 @@   -instance Normed Double where-  type Magnitude Double = Double-  type RealScalar Double = Double-  norm1 = abs-  norm2Sq = abs-  normP _ = abs-  normalize _ _ = 1-  normalize2 _ = 1-  norm2 = abs-  norm2' = abs+-- instance Normed Double where+--   type Magnitude Double = Double+--   type RealScalar Double = Double+--   norm1 = abs+--   norm2Sq = abs+--   normP _ = abs+--   normalize _ _ = 1+--   normalize2 _ = 1+--   norm2 = abs+--   norm2' = abs -instance Normed (Complex Double) where-  type Magnitude (Complex Double) = Double-  type RealScalar (Complex Double) = Double-  norm1 (r :+ i) = abs r + abs i-  norm2Sq = (**2) . magnitude-  normP p (r :+ i) = (r**p + i**p)**(1/p)-  normalize p c = c ./ normP p c-  normalize2 c = c ./ magnitude c-  norm2 = magnitude-  norm2' = magnitude+-- instance Normed (Complex Double) where+--   type Magnitude (Complex Double) = Double+--   type RealScalar (Complex Double) = Double+--   norm1 (r :+ i) = abs r + abs i+--   norm2Sq = (**2) . magnitude+--   normP p (r :+ i) = (r**p + i**p)**(1/p)+--   normalize p c = c ./ normP p c+--   normalize2 c = c ./ magnitude c+--   norm2 = magnitude+--   norm2' = magnitude     @@ -230,7 +241,6 @@   - -- a "sparse matrix ring" ?  -- class MatrixRing m a => SparseMatrixRing m a where@@ -282,66 +292,29 @@  -- * FiniteDim : finite-dimensional objects -class Functor f => FiniteDim f where-  type FDSize f :: *+class FiniteDim f where+  type FDSize f   -- | Dimension (i.e. Int for SpVector, (Int, Int) for SpMatrix)-  dim :: f a -> FDSize f--class FiniteDim' f where-  type FDSize' f :: *-  dim' :: f -> FDSize' f----- -- | unary dimension-checking bracket--- withDim :: (FiniteDim f, Show s) =>---      f e---      -> (FDSize f -> f e -> Bool)---      -> (f e -> c)---      -> String---      -> (f e -> s)---      -> c--- withDim x p f e ef | p (dim x) x = f x---                    | otherwise = error e' where e' = e ++ show (ef x)---- -- | binary dimension-checking bracket--- withDim2 :: (FiniteDim f, FiniteDim g, Show s) =>---      f e---      -> g e---      -> (FDSize f -> FDSize g -> f e -> g e -> Bool)---      -> (f e -> g e -> c)---      -> String---      -> (f e -> g e -> s)---      -> c--- withDim2 x y p f e ef | p (dim x) (dim y) x y = f x y---                       | otherwise = error e' where e' = e ++ show (ef x y)--+  dim :: f -> FDSize f     -- * HasData : accessing inner data (do not export) -class HasData f a where-  type HDData f a :: *+class HasData f where+  type HDData f    -- | Number of nonzeros-  nnz :: f a -> Int-  dat :: f a -> HDData f a+  nnz :: f -> Int+  dat :: f -> HDData f -class HasData' f where-  type HDD f :: *-  nnz' :: f -> Int-  dat' :: f -> HDD f   -- * Sparse : sparse datastructures -class (FiniteDim f, HasData f a) => Sparse f a where+class (FiniteDim f, HasData f) => Sparse f where   -- | Sparsity (fraction of nonzero elements)-  spy :: Fractional b => f a -> b--class (FiniteDim' f, HasData' f) => Sparse' f where-  spy' :: Fractional b => f -> b+  spy :: Fractional b => f -> b   @@ -358,35 +331,28 @@   -- * SpContainer : sparse container datastructures. Insertion, lookup, toList, lookup with 0 default-class Sparse c a => SpContainer c a where+class Sparse c => SpContainer c where   type ScIx c :: *-  scInsert :: ScIx c -> a -> c a -> c a-  scLookup :: c a -> ScIx c -> Maybe a-  scToList :: c a -> [(ScIx c, a)]+  type ScElem c+  scInsert :: ScIx c -> ScElem c -> c -> c+  scLookup :: c -> ScIx c -> Maybe (ScElem c)+  scToList :: c -> [(ScIx c, ScElem c)]   -- -- | Lookup with default, infix form ("safe" : should throw an exception if lookup is outside matrix bounds)-  (@@) :: c a -> ScIx c -> a+  (@@) :: c -> ScIx c -> ScElem c   -class SpContainer' c where-  type ScIx' c :: *-  scInsert' :: ScIx' c -> a -> c -> c-  scLookup' :: c -> ScIx' c -> Maybe a-  scToList' :: c -> [a]-  -- (@@')    -- -- * SparseVector -class SpContainer v e => SparseVector v e where+class SpContainer v => SparseVector v where   type SpvIx v :: *-  svFromList :: Int -> [(SpvIx v, e)] -> v e-  svFromListDense :: Int -> [e] -> v e-  svConcat :: Foldable t => t (v e) -> v e+  svFromList :: Int -> [(SpvIx v, ScElem v)] -> v+  svFromListDense :: Int -> [ScElem v] -> v+  svConcat :: Foldable t => t v -> v   -- svZipWith :: (e -> e -> e) -> v e -> v e -> v e  @@ -398,33 +364,26 @@  -- * SparseMatrix -class SpContainer m e => SparseMatrix m e where-  smFromVector :: LexOrd -> (Int, Int) -> V.Vector (IxRow, IxCol, e) -> m e+class SpContainer m => SparseMatrix m where+  smFromVector :: LexOrd -> (Int, Int) -> V.Vector (IxRow, IxCol, ScElem m) -> m   -- smFromFoldableDense :: Foldable t => t e -> m e  -  smTranspose :: m e -> m e+  smTranspose :: m -> m   -- smExtractSubmatrix ::   --   m e -> (IxRow, IxRow) -> (IxCol, IxCol) -> m e-  encodeIx :: m e -> LexOrd -> (IxRow, IxCol) -> LexIx-  decodeIx :: m e -> LexOrd -> LexIx -> (IxRow, IxCol)+  encodeIx :: m -> LexOrd -> (IxRow, IxCol) -> LexIx+  decodeIx :: m -> LexOrd -> LexIx -> (IxRow, IxCol)   -- data RowsFirst = RowsFirst -- data ColsFirst = ColsFirst --- class SpContainer m e => SparseMatrix m o e where---   smFromVector :: o -> (Int, Int) -> V.Vector (IxRow, IxCol, e) -> m e---   -- smFromFoldableDense :: Foldable t => t e -> m e  ---   smTranspose :: o -> m e -> m e---   -- smExtractSubmatrix ::---   --   m e -> (IxRow, IxRow) -> (IxCol, IxCol) -> m e---   encodeIx :: m e -> o -> (IxRow, IxCol) -> LexIx---   decodeIx :: m e -> o -> LexIx -> (IxRow, IxCol)      + -- * SparseMatVec  -- | Combining functions for relating (structurally) matrices and vectors, e.g. extracting/inserting rows/columns/submatrices@@ -434,11 +393,6 @@ --   smvInsertCol :: m e -> v e -> IxCol -> m e --   smvExtractRow :: m e -> IxRow -> v e --   smvExtractCol :: m e -> IxCol -> v e  -----   

@@ -8,7 +8,6 @@ import Control.Monad.Catch import Control.Monad.State.Strict -import Data.VectorSpace   -- | `eigRayleigh 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.

@@ -7,7 +7,7 @@  import Control.Monad.Catch -import Data.VectorSpace+   -- ** TFQMR

@@ -13,8 +13,8 @@          -- * Linear solvers          -- ** Iterative methods          (<\>),-         -- ** Preconditioners-         jacobiPre, ilu0Pre, mSsorPre,   +         -- -- ** Preconditioners+         -- jacobiPre, ilu0Pre, mSsorPre,             -- ** Moore-Penrose pseudoinverse          pinv,          -- ** Direct methods@@ -29,15 +29,13 @@          eigsArnoldi,          -- * Matrix factorization algorithms          -- ** QR-         qr,+         qr,           -- ** LU          lu,          -- ** Cholesky          chol,           -- ** Arnoldi iteration          arnoldi,    -         -- * Matrix partitioning-         diagPartitions,          -- * Utilities          -- ** Givens' rotation          givens,@@ -45,7 +43,9 @@          conditionNumberSM,          -- ** Householder reflection          hhRefl,-         -- -- * Householder bidiagonalization         +         -- -- * Householder bidiagonalization+         -- ** Matrix partitioning+         diagPartitions,          -- -- * Random arrays          -- randArray,          -- -- * Random matrices and vectors@@ -65,15 +65,23 @@          constv,          -- ** SpMatrix          fromListSM, toListSM,+          -- ** Packing and unpacking, rows or columns of a sparse matrix          -- *** ", using lists as container          fromRowsL, toRowsL,          fromColsL, toColsL,          -- *** ", using Vector as container          fromRowsV, fromColsV,+         -- ** Block operations+         (-=-), (-||-), fromBlocksDiag,+         -- ** Special matrices+         eye, mkDiagonal, mkSubDiagonal,+         permutationSM, permutPairsSM,+         -- * Predicates+         isOrthogonalSM, isDiagonalSM,          -- * Manipulation of sparse data          filterSV, ifilterSV,-         -- * Sparsity-related+         -- * Sparsity-related predicates          nearZero, nearOne, isNz,          -- * Operators          -- ** Scaling@@ -122,28 +130,12 @@  -- import Control.Applicative ((<|>)) --- import Control.Monad (replicateM) import Control.Monad.State.Strict--- import Control.Monad.Writer--- import Control.Monad.Trans.Class-import qualified Control.Monad.Trans.State  as MTS -- (runStateT)--- import Control.Monad.Trans.Writer (runWriterT)+import qualified Control.Monad.Trans.State  as MTS  import Data.Complex -import Data.VectorSpace hiding (magnitude)- import qualified Data.Sparse.Internal.IntM as I--- import Data.Utils.StrictFold (foldlStrict) -- hidden in `containers` --- import qualified System.Random.MWC as MWC--- import qualified System.Random.MWC.Distributions as MWC---- import Data.Monoid--- import qualified Data.Foldable as F--- import qualified Data.Traversable as T---- import qualified Data.List as L- import Data.Maybe  import qualified Data.Vector as V@@ -160,7 +152,8 @@ -- * Matrix condition number  -- | Matrix condition number: computes the QR factorization and extracts the extremal eigenvalues from the R factor-conditionNumberSM :: (MonadThrow m, MatrixRing (SpMatrix a), Num' a) =>+conditionNumberSM :: (MonadThrow m, MonadIO m, MatrixRing (SpMatrix a),+                      PrintDense (SpMatrix a), Num' a) =>      SpMatrix a -> m a conditionNumberSM m = do   (_, r) <- qr m@@ -228,24 +221,24 @@                     -{- $givens-Givens' method, row version: given a row, column pair (i,j), choose a row index i' distinct from i such that :--* i' is below the diagonal--* the corresponding element is nonzero--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 column indices < j.-NB: The Givens' matrix differs from Identity in 4 entries (geometrically, it is a planar rotation embedded in R^n)-NB2: The form of a Complex rotation matrix is as follows (`*` indicates complex conjugation):--@-    ( c    s )- G =(        )-    ( -s*  c*)-@---}+-- | Givens matrix : a planar rotation embedded in R^n. +--+-- >>> aa = fromListSM +-- >>> g <- givens aa 1 0+-- +-- Row version of the method: given a matrix element below the diagonal, indexed (i,j), choose a row index i' that is below the diagonal as well and distinct from i such that the corresponding element is nonzero.+--+-- To zero out entry A(i, j) we must find row i' such that A(i', j) is non-zero but A has zeros in row i' for all column indices < j.+-- +-- NB: The Givens' matrix differs from Identity in 4 entries+-- +-- NB2: The form of a Complex rotation matrix in R^2 is as follows (`*` indicates complex conjugation):+--+-- @+--     ( c    s )+--  G =(        )+--     ( -s*  c*)+-- @ {-# inline givens #-} givens :: (Elt a, MonadThrow m) => SpMatrix a -> IxRow -> IxCol -> m (SpMatrix a) givens aa i j @@ -285,6 +278,7 @@   +     @@ -298,32 +292,17 @@ However, we must also accumulate the `G_i` in order to build `Q`, and the present formulation follows this definition closely. -} {-# inline qr #-}-qr :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+qr :: (Elt a, MatrixRing (SpMatrix a), PrintDense (SpMatrix a),+      Epsilon a, MonadThrow m, MonadIO m) =>      SpMatrix a-     -> m (SpMatrix a, SpMatrix a)  -- ^ Q, R+     -> m (SpMatrix a, SpMatrix a) -- ^ Q, R qr mm = do -     (qt, r, _) <- MTS.execStateT (modifyUntilM haltf qrstepf) gminit-     return (transpose qt, r) -  where-    gminit = (eye (nrows mm), mm, subdiagIndicesSM mm)-    haltf (_, _, iis) = null iis-    qrstepf (qmatt, m, iis) = do-        let (i, j) = head iis-        g <- givens m i j-        let-          qmatt' = g #~# qmatt  -- update Q'-          m' = g #~# m          -- update R-        return (qmatt', m', tail iis)----qr' mm = do       (qt, r, _) <- modifyUntilM' config haltf qrstepf gminit      return (transpose qt, r)    where     gminit = (eye (nrows mm), mm, subdiagIndicesSM mm)     haltf (_, _, iis) = null iis-    config = IterConf 0 True fst2 prd2 where+    config = IterConf 0 False fst2 prd2 where       fst2 (x,y,_) = (x,y)       prd2 (x,y) = do         prd0 x@@ -336,19 +315,18 @@           m' = g #~# m          -- update R         return (qmatt', m', tail iis)     -tm2 :: SpMatrix Double-tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41]     + -- * Eigenvalue algorithms  -- ** QR algorithm  -- | `eigsQR n mm` performs at most `n` iterations of the QR algorithm on matrix `mm`, and returns a SpVector containing all eigenvalues.-eigsQR :: (MonadThrow m, MonadIO m, Num' a, Normed (SpVector a), MatrixRing (SpMatrix a), Typeable (Magnitude (SpVector a))) =>+eigsQR :: (MonadThrow m, MonadIO m, Num' a, Normed (SpVector a), MatrixRing (SpMatrix a), Typeable (Magnitude (SpVector a)), PrintDense (SpVector a), PrintDense (SpMatrix a)) =>         Int            -- ^ Max. # of iterations      -> Bool           -- ^ Print debug information              -> SpMatrix a     -- ^ Operand matrix@@ -369,7 +347,8 @@  -- | `eigsArnoldi n aa b` computes at most n iterations of the Arnoldi algorithm to find a Krylov subspace of (A, b), denoted Q, along with a Hessenberg matrix of coefficients H. After that, it computes the QR decomposition of H, denoted (O, R) and the eigenvalues {λ_i} of A are listed on the diagonal of the R factor. eigsArnoldi :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,-      Elt t, V (SpVector t), Epsilon t, MonadThrow m) =>+      Elt t, V (SpVector t), Epsilon t, PrintDense (SpMatrix t),+      MonadThrow m, MonadIO m) =>      Int      -> SpMatrix t      -> SpVector t@@ -429,14 +408,15 @@  -- | Given a positive semidefinite matrix A, returns a lower-triangular matrix L such that L L^T = A . This is an implementation of the Cholesky–Banachiewicz algorithm, i.e. proceeding row by row from the upper-left corner. -- | NB: The algorithm throws an exception if some diagonal element of A is zero.-chol :: (Elt a, Epsilon a, MonadThrow m) =>+chol :: (Elt a, Epsilon a, MonadThrow m, MonadIO m, PrintDense (SpMatrix a)) =>         SpMatrix a      -> m (SpMatrix a)  -- ^ L chol aa = do   let n = nrows aa       q (i, _) = i == n+      config = IterConf 0 False snd prd0   l0 <- cholUpd aa (0, zeroSM n n)-  (_, lfin) <- MTS.execStateT (modifyUntilM q (cholUpd aa)) l0+  (_, lfin) <- modifyUntilM' config q (cholUpd aa) l0   return lfin   where    oops i = throwM (NeedsPivoting "chol" (unwords ["L", show (i,i)]) :: MatrixException Double)@@ -577,14 +557,15 @@            -tmc4, tmc5 :: SpMatrix (Complex Double)+tmc4, tmc5, tmc6 :: SpMatrix (Complex Double) tmc4 = fromListDenseSM 3 [3:+1, 4:+(-1), (-5):+3, 2:+2, 3:+(-2), 5:+0.2, 7:+(-2), 9:+(-1), 2:+3]  tvc4 = vc [1:+3,2:+2,1:+9] -tmc5 = fromListDenseSM 4 $ zipWith (:+) [16..31] [15,14..0]+tmc5 = fromListDenseSM 4 $ zipWith (:+) [16..31] [17,14..]  +tmc6 = fromListDenseSM 2 $ zipWith (:+) [1,2,3,4] [4,3,2,1]    @@ -717,25 +698,38 @@  -- * Linear solver, LU-based +luSolveConfig :: PrintDense (SpVector t) => IterationConfig (SpVector t, IxRow) (SpVector t)+luSolveConfig = IterConf 0 False fst prd0+++ -- | Direct solver based on a triangular factorization of the system matrix.-luSolve :: (Scalar (SpVector t) ~ t, MonadThrow m, Elt t, InnerSpace (SpVector t), Epsilon t) =>+luSolve :: (Scalar (SpVector t) ~ t, MonadThrow m, Elt t, InnerSpace (SpVector t),+            Epsilon t, PrintDense (SpVector t), MonadIO m) =>      SpMatrix t    -- ^ Lower triangular      -> SpMatrix t -- ^ Upper triangular      -> SpVector t -- ^ r.h.s.      -> m (SpVector t) luSolve ll uu b   | isLowerTriSM ll && isUpperTriSM uu = do-      w <- triLowerSolve ll b-      triUpperSolve uu w+      w <- triLowerSolve0 luSolveConfig ll b+      triUpperSolve0 luSolveConfig uu w   | otherwise = throwM (NonTriangularException "luSolve")  -- | Forward substitution solver-triLowerSolve :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),-      Epsilon t, MonadThrow m) =>-     SpMatrix t    -- ^ Lower triangular-     -> SpVector t -- ^ r.h.s-     -> m (SpVector t)-triLowerSolve ll b = do+triLowerSolve+  :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),+      PrintDense (SpVector t), Epsilon t, MonadThrow m, MonadIO m) =>+      SpMatrix t -> SpVector t -> m (SpVector t)+triLowerSolve = triLowerSolve0 luSolveConfig++++triLowerSolve0 :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),+      Epsilon t, MonadThrow m, MonadIO m) =>+        IterationConfig (SpVector t, IxRow) b+     -> SpMatrix t -> SpVector t -> m (SpVector t)  +triLowerSolve0 config ll b = do   let q (_, i) = i == nb       nb = svDim b       oops i = throwM (NeedsPivoting "triLowerSolve" (unwords ["L", show (i, i)]) :: MatrixException Double)@@ -757,19 +751,24 @@           then return (insertSpVector 0 w0 $ zeroSV (dim b), 1)           else oops (0 :: Int)   l0 <- lInit             -  (v, _) <- MTS.execStateT (modifyUntilM q lStep) l0+  (v, _) <- modifyUntilM' config q lStep l0   return $ sparsifySV v    -- NB in the computation of `xi` we must rebalance the subrow indices (extractSubRow_RK) because `dropSV` does that too, in order to take the inner product with consistent index pairs -- | Backward substitution solver-triUpperSolve :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),-      Epsilon t, MonadThrow m) =>-     SpMatrix t    -- ^ Upper triangular-     -> SpVector t -- ^ r.h.s-     -> m (SpVector t)-triUpperSolve uu w = do +triUpperSolve+  :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),+      PrintDense (SpVector t), Epsilon t, MonadThrow m, MonadIO m) =>+      SpMatrix t -> SpVector t -> m (SpVector t)+triUpperSolve = triUpperSolve0 luSolveConfig++triUpperSolve0 :: (Scalar (SpVector t) ~ t, Elt t, InnerSpace (SpVector t),+      Epsilon t, MonadThrow m, MonadIO m) =>+     IterationConfig (SpVector t, IxRow) b+     -> SpMatrix t -> SpVector t -> m (SpVector t)  +triUpperSolve0 conf uu w = do    let q (_, i) = i == (- 1)       nw = svDim w       oops i = throwM (NeedsPivoting "triUpperSolve" (unwords ["U", show (i, i)]) :: MatrixException Double)@@ -790,7 +789,7 @@           then return (insertSpVector i x0 (zeroSV nw), i - 1)           else oops (0 :: Int)   u0 <- uInit             -  (x, _) <- MTS.execStateT (modifyUntilM q uStep) u0+  (x, _) <- modifyUntilM' conf q uStep u0   return $ sparsifySV x        @@ -817,10 +816,10 @@ -- -- A common optimization involves interleaving the QR factorization (and the subsequent triangular solve) with the Arnoldi process (and employing an updating QR factorization which only requires one Givens' rotation at every update).  -gmres :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,-      Elt t, Normed (SpVector t), LinearVectorSpace (SpVector t), Epsilon t,-      MonadThrow m) =>-     SpMatrix t -> SpVector t -> m (SpVector t)+-- gmres :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,+--       Elt t, Normed (SpVector t), LinearVectorSpace (SpVector t), Epsilon t,+--       MonadThrow m) =>+--      SpMatrix t -> SpVector t -> m (SpVector t) gmres aa b = do   let m = ncols aa   (qa, ha) <- arnoldi aa b m   -- at most m steps of Arnoldi (aa, b)@@ -999,7 +998,7 @@        BICGSTAB_ -> solver "BICGSTAB" nits _xBicgstab (bicgstabStep aa r0hat) (bicgsInit aa b x0)        BCG_ -> solver "BCG" nits _xBcg (bcgStep aa) (bcgInit aa b x0)        CGS_ -> solver "CGS" nits _x  (cgsStep aa r0hat) (cgsInit aa b x0)-       GMRES_ -> gmres aa b+       GMRES_ -> gmres aa b                    CGNE_ -> solver "CGNE" nits _xCgne (cgneStep aa) (cgneInit aa b x0)      r0hat = b ^-^ (aa #> x0)      nits = 200@@ -1008,21 +1007,23 @@      solver fname nitermax fproj stepf initf = do       xf <- untilConvergedG fname config (const True) stepf initf       return $ fproj xf-      where-        config = IterConf nitermax False fproj prd0+       where+         config =  IterConf nitermax False fproj prd0++               -- | <\> uses the GMRES method as default-instance LinearSystem (SpVector Double) where-  aa <\> b = linSolve0 GMRES_ aa b (mkSpVR n $ replicate n 0.1)-    where n = ncols aa -instance LinearSystem (SpVector (Complex Double)) where-  aa <\> b = linSolve0 GMRES_ aa b (mkSpVC n $ replicate n 0.1)-    where n = ncols aa+-- instance LinearSystem (SpVector Double) where+--   aa <\> b = linSolve0 GMRES_ aa b (mkSpVR n $ replicate n 0.1)+--     where n = ncols aa +-- instance LinearSystem (SpVector (Complex Double)) where+--   aa <\> b = linSolve0 GMRES_ aa b (mkSpVC n $ replicate n 0.1)+--     where n = ncols aa   @@ -1033,6 +1034,7 @@   + -- | TODO : if system is poorly conditioned, is it better to warn the user or just switch solvers (e.g. via the pseudoinverse) ?  -- linSolveQR aa b init f1 stepf@@ -1055,12 +1057,12 @@  -- test data --- aa4 : eigenvalues 1 (mult.=2) and -1-aa4 :: SpMatrix Double-aa4 = fromListDenseSM 3 [3,2,-2,2,2,-1,6,5,-4] +-- -- aa4 : eigenvalues 1 (mult.=2) and -1+-- aa4 :: SpMatrix Double+-- aa4 = fromListDenseSM 3 [3,2,-2,2,2,-1,6,5,-4]  -aa4c :: SpMatrix (Complex Double)-aa4c = toC <$> aa4+-- aa4c :: SpMatrix (Complex Double)+-- aa4c = toC <$> aa4   

@@ -20,14 +20,12 @@ -- import Control.Monad (replicateM) import Control.Monad.Catch -- (MonadThrow (..)) import Control.Monad.IO.Class-import Control.Monad.State+-- import Control.Monad.State -- import Data.Foldable (foldrM)  import Data.Complex -- import Data.Either (either)-import Data.Typeable--import Data.VectorSpace hiding (magnitude)+-- import Data.Typeable  -- import Control.Monad.State.Strict (execState) @@ -42,23 +40,22 @@ main :: IO () main = hspec spec --- niter = 5- spec :: Spec spec = do   describe "Numeric.LinearAlgebra.Sparse : Library" $ do-    prop "Subtraction is cancellative" $ \(x :: SpVector Double) ->-      norm2Sq (x ^-^ x) `shouldBe` zeroV-    it "<.> : inner product (Real)" $-      tv0 <.> tv0 `shouldBe` 61-    it "<.> : inner product (Complex)" $-      tvc2 <.> tvc3 `shouldBe` 2 :+ (-2)  +    -- prop "Subtraction is cancellative" $ \(x :: SpVector Double) ->+    --   norm2Sq (x ^-^ x) `shouldBe` zeroV+    -- it "<.> : inner product (Real)" $+    --   tv0 <.> tv0 `shouldBe` 61+    -- it "<.> : inner product (Complex)" $+    --   tvc2 <.> tvc3 `shouldBe` 2 :+ (-2)       it "transpose : sparse matrix transpose" $       transpose m1 `shouldBe` m1t-    it "(#>) : matrix-vector product (Real)" $-      nearZero ( norm2Sq ((aa0 #> x0true) ^-^ b0 )) `shouldBe` True-    it "(<#) : vector-matrix product (Real)" $-      nearZero ( norm2Sq ((x0true <# aa0) ^-^ aa0tx0 ))`shouldBe` True  +    -- it "(#>) : matrix-vector product (Real)" $+    --   nearZero ( norm2Sq ((aa0 #> x0true) ^-^ b0 )) `shouldBe` True+    -- it "(<#) : vector-matrix product (Real)" $+    --   nearZero ( norm2Sq ((x0true <# aa0) ^-^ aa0tx0 ))`shouldBe` True+       it "(##) : matrix-matrix product (Real, square)" $        (m1 ## m2) `shouldBe` m1m2     it "(##) : matrix-matrix product (Real, rectangular)" $ do@@ -85,80 +82,82 @@     it "isLowerTriSM : checks whether matrix is lower triangular" $       isLowerTriSM tm8' && isUpperTriSM tm8 `shouldBe` True       -    it "untilConvergedG0 : early termination by iteration count and termination by convergence" $ -     let-      n1 = 4-      nexp1 = fromIntegral n1 / fromIntegral (2^n1) -- 0.25-      f x = x/2-      mm1 = untilConvergedG0 "blah"-               (IterConf n1 False id print) (1/(2^n1)) f (fromIntegral n1 :: Double)-      n2 = 2^16-      mm2 = untilConvergedG0 "blah"-               (IterConf n2 False id print) (1/(2^n2)) f (fromIntegral n1 :: Double)-      eh (NotConvergedE _ _ x) = return x-      in-       do x1 <- mm1 `catch` eh-          x1 `shouldBe` nexp1-          x2 <- mm2 `catch` eh-          nearZero x2 `shouldBe` True+    -- it "untilConvergedG0 : early termination by iteration count and termination by convergence" $ +    --  let+    --   n1 = 4+    --   nexp1 = fromIntegral n1 / fromIntegral (2^n1) -- 0.25+    --   f x = x/2+    --   mm1 = untilConvergedG0 "blah"+    --            (IterConf n1 False id print) (1/(2^n1)) f (fromIntegral n1 :: Double)+    --   n2 = 2^16+    --   mm2 = untilConvergedG0 "blah"+    --            (IterConf n2 False id print) (1/(2^n2)) f (fromIntegral n1 :: Double)+    --   eh (NotConvergedE _ _ x) = return x+    --   in+    --    do x1 <- mm1 `catch` eh+    --       x1 `shouldBe` nexp1+    --       x2 <- mm2 `catch` eh+    --       nearZero x2 `shouldBe` True              -  describe "QuickCheck properties:" $ do-    -- prop "prop_matSPD_vec : (m #^# m) is symmetric positive definite" $-    --   \(PropMatSPDVec (m :: SpMatrix Double) v) -> prop_spd m v-    prop "prop_dot : (v <.> v) ~= 1 if ||v|| == 1" $-      \(v :: SpVector Double) -> prop_dot v-    prop "prop_matMat1 : (A ## B)^T == (B^T ## A^T)" $-      \p@(PropMatMat (_ :: SpMatrix Double) _) -> prop_matMat1 p-    prop "prop_matMat2 : M^T ##^ M == M #^# M^T" $-      \p@(PropMat (_ :: SpMatrix Double)) -> prop_matMat2 p-    -- prop "prop_matMat2 : M^T ##^ M == M #^# M^T , Complex" $-    --   \p@(PropMat (_ :: SpMatrix (Complex Double))) -> whenFail (prd $ unPropMat p) (prop_matMat2 p :: Bool)+  -- describe "QuickCheck properties:" $ do+  --   prop "prop_matSPD_vec : (m #^# m) is symmetric positive definite" $+  --     \(PropMatSPDVec (m :: SpMatrix Double) v) -> prop_spd m v+  --   -- prop "prop_matSPD_vec : (m #^# m) is symmetric positive definite" $+  --   --   \(PropMatSPDVec (m :: SpMatrix (Complex Double)) v) -> prop_spd m v  +  --   prop "prop_dot : (v <.> v) ~= 1 if ||v|| == 1" $+  --     \(v :: SpVector Double) -> prop_dot v+  --   prop "prop_matMat1 : (A ## B)^T == (B^T ## A^T)" $+  --     \p@(PropMatMat (_ :: SpMatrix Double) _) -> prop_matMat1 p+  --   prop "prop_matMat2 : M^T ##^ M == M #^# M^T" $+  --     \p@(PropMat (_ :: SpMatrix Double)) -> prop_matMat2 p+  --   -- prop "prop_matMat2 : M^T ##^ M == M #^# M^T , Complex" $+  --   --   \p@(PropMat (_ :: SpMatrix (Complex Double))) -> whenFail (prd $ unPropMat p) (prop_matMat2 p :: Bool)                                                        -    -- prop "prop_QR : Q R = A, Q is orthogonal, R is upper triangular" $-    --   \p@(PropMatI (_ :: SpMatrix Double)) -> prop_QR p-    -- prop "prop_Cholesky" $ \p@(PropMat_SPD (_ :: SpMatrix Double)) -> prop_Cholesky p-    -- prop "prop_linSolve GMRES" $ prop_linSolve GMRES_-      -- prop "aa2 is positive semidefinite" $ \(v :: SpVector Double) ->-    --   prop_psd aa2 v+  --   -- prop "prop_QR : Q R = A, Q is orthogonal, R is upper triangular" $+  --   --   \p@(PropMatI (_ :: SpMatrix Double)) -> prop_QR p+  --   -- prop "prop_Cholesky" $ \p@(PropMat_SPD (_ :: SpMatrix Double)) -> prop_Cholesky p+  --   -- prop "prop_linSolve GMRES" $ prop_linSolve GMRES_+  --     -- prop "aa2 is positive semidefinite" $ \(v :: SpVector Double) ->+  --   --   prop_psd aa2 v     -  describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers (Real)" $ do-    -- it "TFQMR (2 x 2 dense)" $-    it "GMRES (2 x 2 dense)" $-      -- checkLinSolveR GMRES_ aa0 b0 x0true >>= (`shouldBe` True)-      checkLinSolveR GMRES_ aa0 b0 x0true `shouldBeM` True-    it "GMRES (3 x 3 sparse, symmetric pos.def.)" $-      checkLinSolveR GMRES_ aa2 b2 x2 >>= (`shouldBe` True)-    it "GMRES (4 x 4 sparse)" $-      checkLinSolveR GMRES_ aa1 b1 x1 >>= (`shouldBe` True)-    it "GMRES (5 x 5 sparse)" $-      checkLinSolveR GMRES_ tm7 tvb7 tvx7 >>= (`shouldBe` True)  -    it "BCG (2 x 2 dense)" $-      checkLinSolveR BCG_ aa0 b0 x0true >>= (`shouldBe` True)-    it "BCG (3 x 3 sparse, symmetric pos.def.)" $-      checkLinSolveR BCG_ aa2 b2 x2 >>= (`shouldBe` True)-    -- it "BiCGSTAB (2 x 2 dense)" $ -    --   nearZero (normSq (linSolve BICGSTAB_ aa0 b0 ^-^ x0true)) `shouldBe` True-    it "BiCGSTAB (3 x 3 sparse, symmetric pos.def.)" $ -      checkLinSolveR BICGSTAB_ aa2 b2 x2 >>= (`shouldBe` True)-    it "CGS (2 x 2 dense)" $ -      checkLinSolveR CGS_ aa0 b0 x0true >>= (`shouldBe` True)-    it "CGS (3 x 3 sparse, SPD)" $ -      checkLinSolveR CGS_ aa2 b2 x2 >>= (`shouldBe` True)-    it "Moore-Penrose pseudoinverse (3 x 2 dense)" $-      checkPinv aa10 b10 x10 >>= (`shouldBe` True)+  -- describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers (Real)" $ do+  --   -- it "TFQMR (2 x 2 dense)" $+  --   it "GMRES (2 x 2 dense)" $+  --     -- checkLinSolveR GMRES_ aa0 b0 x0true >>= (`shouldBe` True)+  --     checkLinSolveR GMRES_ aa0 b0 x0true `shouldBeM` True+  --   it "GMRES (3 x 3 sparse, symmetric pos.def.)" $+  --     checkLinSolveR GMRES_ aa2 b2 x2 >>= (`shouldBe` True)+  --   it "GMRES (4 x 4 sparse)" $+  --     checkLinSolveR GMRES_ aa1 b1 x1 >>= (`shouldBe` True)+  --   it "GMRES (5 x 5 sparse)" $+  --     checkLinSolveR GMRES_ tm7 tvb7 tvx7 >>= (`shouldBe` True)  +  --   it "BCG (2 x 2 dense)" $+  --     checkLinSolveR BCG_ aa0 b0 x0true >>= (`shouldBe` True)+  --   it "BCG (3 x 3 sparse, symmetric pos.def.)" $+  --     checkLinSolveR BCG_ aa2 b2 x2 >>= (`shouldBe` True)+  --   -- it "BiCGSTAB (2 x 2 dense)" $ +  --   --   nearZero (normSq (linSolve BICGSTAB_ aa0 b0 ^-^ x0true)) `shouldBe` True+  --   it "BiCGSTAB (3 x 3 sparse, symmetric pos.def.)" $ +  --     checkLinSolveR BICGSTAB_ aa2 b2 x2 >>= (`shouldBe` True)+  --   it "CGS (2 x 2 dense)" $ +  --     checkLinSolveR CGS_ aa0 b0 x0true >>= (`shouldBe` True)+  --   it "CGS (3 x 3 sparse, SPD)" $ +  --     checkLinSolveR CGS_ aa2 b2 x2 >>= (`shouldBe` True)+  --   it "Moore-Penrose pseudoinverse (3 x 2 dense)" $+  --     checkPinv aa10 b10 x10 >>= (`shouldBe` True)          -- describe "Numeric.LinearAlgebra.Sparse : Iterative linear solvers (Complex)" $ do   --   it "<\\> (3 x 3 dense)" $   --     checkBackslash tmc4 tvc4 >>= (`shouldBe` True)       -  describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers (Real)" $ -    it "luSolve (4 x 4 sparse)" $ -      checkLuSolve aa1 b1 >>= (`shouldBe` (True, True, True))-  -- describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers (Complex)" $ -  --   it "luSolve (3 x 3 dense)" $ -  --     checkLuSolve tmc4 tvc4 >>= (`shouldBe` (True, True, True)) +  -- describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers (Real)" $ +  --   it "luSolve (4 x 4 sparse)" $ +  --     checkLuSolve aa1 b1 >>= (`shouldBe` (True, True, True))+  -- -- describe "Numeric.LinearAlgebra.Sparse : Direct linear solvers (Complex)" $ +  -- --   it "luSolve (3 x 3 dense)" $ +  -- --     checkLuSolve tmc4 tvc4 >>= (`shouldBe` (True, True, True))           describe "Numeric.LinearAlgebra.Sparse : QR factorization (Real)" $ do     it "qr (3 x 3 dense)" $ @@ -172,18 +171,16 @@       checkQr aa3cx >>= (`shouldBe` True)     it "qr (3 x 3 dense)" $       checkQr tmc4 >>= (`shouldBe` True)  -  describe "Numeric.LinearAlgebra.Sparse : LU factorization (Real)" $ do-    it "lu (3 x 3 dense)" $-      checkLu tm2 >>= (`shouldBe` True)-    it "lu (4 x 4 dense)" $-      checkLu tm6 >>= (`shouldBe` True)-    it "lu (5 x 5 sparse)" $-      checkLu tm7 >>= (`shouldBe` True)-  describe "Numeric.LinearAlgebra.Sparse : LU factorization (Complex)" $ do-    it "lu (3 x 3 dense)" $-      checkLu tmc4 >>= (`shouldBe` True)-    -- it "lu (4 x 4 dense)" $-    --   checkLu tmc5 >>= (`shouldBe` True)   +  -- describe "Numeric.LinearAlgebra.Sparse : LU factorization (Real)" $ do+  --   -- it "lu (3 x 3 dense)" $+  --   --   checkLu tm2 >>= (`shouldBe` True)+  --   it "lu (4 x 4 dense)" $+  --     checkLu tm6 >>= (`shouldBe` True)+  --   it "lu (5 x 5 sparse)" $+  --     checkLu tm7 >>= (`shouldBe` True)+  -- describe "Numeric.LinearAlgebra.Sparse : LU factorization (Complex)" $ +  --   it "lu (3 x 3 dense)" $+  --     checkLu tmc4 >>= (`shouldBe` True)      describe "Numeric.LinearAlgebra.Sparse : Cholesky factorization (Real, symmetric pos.def.)" $      it "chol (5 x 5 sparse)" $@@ -193,14 +190,14 @@   --     checkChol (tmc5 ##^ tmc5) >>= (`shouldBe` True)           -  describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration (Real)" $ do      -    it "arnoldi (4 x 4 dense)" $-      checkArnoldi tm6 4 >>= (`shouldBe` True)-    it "arnoldi (5 x 5 sparse)" $-      checkArnoldi tm7 5 >>= (`shouldBe` True)-  -- describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration (Complex)" $ do      +  -- describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration (Real)" $ do         --   it "arnoldi (4 x 4 dense)" $-  --     checkArnoldi tmc4 4 >>= (`shouldBe` True)      +  --     checkArnoldi tm6 4 >>= (`shouldBe` True)+  --   it "arnoldi (5 x 5 sparse)" $+  --     checkArnoldi tm7 5 >>= (`shouldBe` True)+  -- -- describe "Numeric.LinearAlgebra.Sparse : Arnoldi iteration (Complex)" $ do      +  -- --   it "arnoldi (4 x 4 dense)" $+  -- --     checkArnoldi tmc4 4 >>= (`shouldBe` True)         {- linear systems -}@@ -220,27 +217,27 @@   xhat <- linSolve0 method aa b x0r   return $ nearZero $ norm2 (x ^-^ xhat) -checkLinSolveR-  :: (MonadIO m, MonadCatch m) =>-     LinSolveMethod -     -> SpMatrix Double        -- ^ operator-     -> SpVector Double        -- ^ r.h.s-     -> SpVector Double        -- ^ candidate solution-     -> m Bool-checkLinSolveR method aa b x = checkLinSolve method aa b x x0r where-  x0r = mkSpVR n $ replicate n 0.1-  n = ncols aa+-- checkLinSolveR+--   :: (MonadIO m, MonadCatch m) =>+--      LinSolveMethod +--      -> SpMatrix Double        -- ^ operator+--      -> SpVector Double        -- ^ r.h.s+--      -> SpVector Double        -- ^ candidate solution+--      -> m Bool+-- checkLinSolveR method aa b x = checkLinSolve method aa b x x0r where+--   x0r = mkSpVR n $ replicate n 0.1+--   n = ncols aa -checkLinSolveC-  :: (MonadIO m, MonadCatch m) =>-     LinSolveMethod-     -> SpMatrix (Complex Double)-     -> SpVector (Complex Double)-     -> SpVector (Complex Double)-     -> m Bool-checkLinSolveC method aa b x = checkLinSolve method aa b x x0r where-  x0r = mkSpVC n $ replicate n (0.1 :+ 0.1)-  n = ncols aa+-- checkLinSolveC+--   :: (MonadIO m, MonadCatch m) =>+--      LinSolveMethod+--      -> SpMatrix (Complex Double)+--      -> SpVector (Complex Double)+--      -> SpVector (Complex Double)+--      -> m Bool+-- checkLinSolveC method aa b x = checkLinSolve method aa b x x0r where+--   x0r = mkSpVC n $ replicate n (0.1 :+ 0.1)+--   n = ncols aa   @@ -275,16 +272,25 @@  {- QR-} -checkQr :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>-     SpMatrix a -> m Bool-checkQr a = do-  (q, r) <- qr a+checkQr :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, PrintDense (SpMatrix a),+            MonadThrow m, MonadIO m) =>+     SpMatrix a+     -> m Bool+checkQr = checkQr0 qr++++checkQr0 :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+     (SpMatrix a -> m (SpMatrix a, SpMatrix a))+     -> SpMatrix a+     -> m Bool+checkQr0 mfqr a = do+  (q, r) <- mfqr a   let c1 = nearZero $ normFrobenius $ sparsifySM ((q ## r) ^-^ a)       c2 = isOrthogonalSM q       c3 = isUpperTriSM r   return $ c1 && c2 && c3 - -- stepQR a = do --   (q, r) <- qr a --   return $ r #~# q@@ -310,7 +316,8 @@  {- Cholesky -} -checkChol :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+checkChol :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, PrintDense (SpMatrix a),+              MonadThrow m, MonadIO m) =>      SpMatrix a -> m Bool checkChol a = do -- c1 && c2 where   l <- chol a@@ -323,7 +330,8 @@  checkLuSolve :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,       Elt t, Normed (SpVector t), LinearVectorSpace (SpVector t),-      Epsilon t, MonadThrow m) =>+      PrintDense (SpVector t),      +      Epsilon t, MonadThrow m, MonadIO m) =>      SpMatrix t -> SpVector t -> m (Bool, Bool, Bool) checkLuSolve amat rhs = do   (lmat, umat) <- lu amat@@ -347,9 +355,11 @@   {- triangular solvers -}+ checkTriUpperSolve :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,       Elt t, Normed (SpVector t), LinearVectorSpace (SpVector t), Epsilon t,-      MonadThrow m) =>+      PrintDense (SpVector t),+      MonadThrow m, MonadIO m) =>      SpMatrix t -> SpVector t -> m (SpVector t, Bool) checkTriUpperSolve umat rhs = do   xhat <- triUpperSolve umat rhs@@ -359,7 +369,8 @@  checkTriLowerSolve :: (Scalar (SpVector t) ~ t, MatrixType (SpVector t) ~ SpMatrix t,       Elt t, Normed (SpVector t), LinearVectorSpace (SpVector t), Epsilon t,-      MonadThrow m) =>+      PrintDense (SpVector t),      +      MonadThrow m, MonadIO m) =>      SpMatrix t -> SpVector t -> m (SpVector t, Bool) checkTriLowerSolve lmat rhs = do   xhat <- triLowerSolve lmat rhs@@ -497,27 +508,9 @@   -     ----------------- -- | Random sparse vector genSpV0 :: Arbitrary a => Int -> Int -> Gen (SpVector a) genSpV0 d n = do@@ -592,46 +585,46 @@   --- | An arbitrary Householder reflection matrix -newtype PropHhReflMat a = PropHhReflMat { unHh :: SpMatrix a} deriving Show-instance Arbitrary (PropHhReflMat Double) where-  arbitrary = sized (\n -> PropHhReflMat <$> genReflMat n) `suchThat` ((> 2) . nrows . unHh)+-- -- | An arbitrary Householder reflection matrix +-- newtype PropHhReflMat a = PropHhReflMat { unHh :: SpMatrix a} deriving Show+-- instance Arbitrary (PropHhReflMat Double) where+--   arbitrary = sized (\n -> PropHhReflMat <$> genReflMat n) `suchThat` ((> 2) . nrows . unHh) -genReflMat :: Int -> Gen (SpMatrix Double)-genReflMat n = do-  v <- normalize2 <$> (genSpVDense n :: Gen (SpVector Double))-  return $ hhRefl v+-- genReflMat :: Int -> Gen (SpMatrix Double)+-- genReflMat n = do+--   v <- normalize2 <$> (genSpVDense n :: Gen (SpVector Double))+--   return $ hhRefl v    --- a product of a "large" number of random Householder matrices-newtype PropMatUnitary a = PropMatUnitary {unUnitary :: SpMatrix a} deriving Show-instance Arbitrary (PropMatUnitary Double) where-  arbitrary = sized (\n -> PropMatUnitary <$> genReflMat n) `suchThat` ((> 2) . nrows . unUnitary)+-- -- a product of a "large" number of random Householder matrices+-- newtype PropMatUnitary a = PropMatUnitary {unUnitary :: SpMatrix a} deriving Show+-- instance Arbitrary (PropMatUnitary Double) where+--   arbitrary = sized (\n -> PropMatUnitary <$> genReflMat n) `suchThat` ((> 2) . nrows . unUnitary)   -genUnitary :: Int -> Gen (SpMatrix Double)-genUnitary n = do-  q <- vectorOf (10 * n) $ genReflMat n  -- random Householder matrix-  return $ foldr (##) (eye n) q +-- genUnitary :: Int -> Gen (SpMatrix Double)+-- genUnitary n = do+--   q <- vectorOf (10 * n) $ genReflMat n  -- random Householder matrix+--   return $ foldr (##) (eye n) q  -prop_unitary :: (MatrixRing (SpMatrix a), Epsilon a, Eq a) =>-     PropMatUnitary a -> Bool-prop_unitary (PropMatUnitary m) = isOrthogonalSM m+-- prop_unitary :: (MatrixRing (SpMatrix a), Epsilon a, Eq a) =>+--      PropMatUnitary a -> Bool+-- prop_unitary (PropMatUnitary m) = isOrthogonalSM m    --- | A symmetric, positive-definite matrix-newtype PropMatSPD a = PropMatSPD {unPropMatSPD :: SpMatrix a} deriving (Show)-instance Arbitrary (PropMatSPD Double) where-  arbitrary = sized genSpM_SPD `suchThat` ((> 2) . nrows . unPropMatSPD)+-- -- | A symmetric, positive-definite matrix+-- newtype PropMatSPD a = PropMatSPD {unPropMatSPD :: SpMatrix a} deriving (Show)+-- instance Arbitrary (PropMatSPD Double) where+--   arbitrary = sized genSpM_SPD `suchThat` ((> 2) . nrows . unPropMatSPD) -genSpM_SPD :: Int -> Gen (PropMatSPD Double)-genSpM_SPD n = do-  q <- genUnitary n -  d <- genSpMDiagonal (all (> 0)) n          -- positive diagonal-  return $ PropMatSPD ( q ## (d ##^ q) )+-- genSpM_SPD :: Int -> Gen (PropMatSPD Double)+-- genSpM_SPD n = do+--   q <- genUnitary n +--   d <- genSpMDiagonal (all (> 0)) n          -- positive diagonal+--   return $ PropMatSPD ( q ## (d ##^ q) )   @@ -660,12 +653,12 @@   --- -- | A symmetric positive definite matrix and vector of compatible size--- data PropMatSPDVec a = PropMatSPDVec (SpMatrix a) (SpVector a) deriving (Eq, Show)--- instance Arbitrary (PropMatSPDVec Double) where---   arbitrary = do---     PropMatVec m v <- arbitrary -- :: Gen (PropMatVec Double)---     return $ PropMatSPDVec (m #^# m) v+-- | A symmetric positive definite matrix and vector of compatible size+data PropMatSPDVec a = PropMatSPDVec (SpMatrix a) (SpVector a) deriving (Eq, Show)+instance Arbitrary (PropMatSPDVec Double) where+  arbitrary = do+    PropMatVec m v <- arbitrary -- :: Gen (PropMatVec Double)+    return $ PropMatSPDVec (m #^# m) v       @@ -715,14 +708,16 @@   --- | Cholesky factorization of a random SPD matrix -prop_Cholesky :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>-     PropMatSPD a -> m Bool-prop_Cholesky (PropMatSPD m) = checkChol m+-- -- | Cholesky factorization of a random SPD matrix +-- prop_Cholesky :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, PrintDense (SpMatrix a),+--                   MonadThrow m, MonadIO m) =>+--      PropMatSPD a -> m Bool+-- prop_Cholesky (PropMatSPD m) = checkChol m   -- | QR decomposition-prop_QR :: (Elt a, MatrixRing (SpMatrix a), Epsilon a, MonadThrow m) =>+prop_QR :: (Elt a, MatrixRing (SpMatrix a), PrintDense (SpMatrix a), Epsilon a,+            MonadThrow m, MonadIO m) =>      PropMatI a -> m Bool prop_QR (PropMatI m) = checkQr m @@ -735,13 +730,120 @@ --     b = aai #> x --   checkLinSolveR method aai b x --- -- test data    ++++++++++++-- --+++-- -- run N iterations ++-- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB+-- runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where+--    x0 = mkSpVectorD nd $ replicate nd 0.9+--    nd = dim r0+--    r0 = b ^-^ (aa #> x0)    +--    p0 = r0+--    bicgsInit = BICGSTAB x0 r0 p0++-- -- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS+-- runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where+--   x0 = mkSpVectorD nd $ replicate nd 0.1+--   nd = dim r0+--   r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution+--   p0 = r0+--   u0 = r0+--   cgsInit = CGS x0 r0 p0 u0++++-- solveRandomN ndim nsp niter = do+--   aa0 <- randSpMat ndim (nsp ^ 2)+--   let aa = aa0 ^+^ eye ndim+--   xtrue <- randSpVec ndim nsp+--   let b = aa #> xtrue+--       xhatB = head $ runNBiC niter aa b+--       xhatC = head $ runNCGS niter aa b+--   -- printDenseSM aa    +--   return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))+++ {-+random linear system +-}++++-- -- dense+-- solveRandom n = do+--   aa0 <- randMat n+--   let aa = aa0 ^+^ eye n+--   xtrue <- randVec n+--   -- x0 <- randVec n+--   let b = aa #> xtrue+--       dx = aa <\> b ^-^ xtrue+--   return $ normSq dx+--   -- let xhatB = _xBicgstab (bicgstab aa b x0 x0)+--   --     xhatC = _x (cgs aa b x0 x0)+--   -- return (aa, x, x0, b, xhatB, xhatC)++-- -- sparse+-- solveSpRandom :: Int -> Int -> IO Double+-- solveSpRandom n nsp = do+--   aa0 <- randSpMat n nsp+--   let aa = aa0 ^+^ eye n+--   xtrue <- randSpVec n nsp+--   let b = (aa ^+^ eye n) #> xtrue+--       dx = aa <\> b ^-^ xtrue+--   return $ normSq dx+++++-- solveRandomBanded n bw mu sig = do+--   let ndiags = 2*bw+--   bands <- replicateM (ndiags + 1) (randArray n mu sig)+--   xtrue <- randVec n+--   b <- randVec n+--   let+--     diags = [-bw .. bw - 1]++-- randDiagMat :: PrimMonad m =>+--      Rows -> Double -> Double -> Int -> m (SpMatrix Double)+-- randDiagMat n mu sig i = do+--   x <- randArray n mu sig+--   return $ mkSubDiagonal n i x++++++++++++++-- | Data++{-+ example 0 : 2x2 linear system  [1 2] [2] = [8]@@ -914,6 +1016,22 @@  tm2 = fromListDenseSM 3 [12, 6, -4, -51, 167, 24, 4, -68, -41] +-- λ> (l,u) <- lu tm2+-- λ> prd l++-- 1.00   , _      , _      +-- 0.50   , 1.00   , _      +-- -0.33  , 0.04   , 1.00   ++-- λ> prd u++-- 12.00  , -51.00 , 4.00   +-- _      , 1.92e2 , -70.00 +-- _      , _      , -37.12 ++++ tm3 = transposeSM $ fromListDenseSM 3 [1 .. 9]  @@ -937,7 +1055,7 @@  tvx7 = mkSpVR 5 [3,8,-12,4,9] -tvb7 = tm7 #> tvx7+-- tvb7 = tm7 #> tvx7   @@ -973,11 +1091,13 @@ -- Complex linear system  tmc4,tmc5, tmc6 :: SpMatrix (Complex Double)+-- tmc4: condition number = 4.4233 tmc4 = fromListDenseSM 3 [3:+1, 4:+(-1), (-5):+3, 2:+2, 3:+(-2), 5:+0.2, 7:+(-2), 9:+(-1), 2:+3] +-- tvc4 : unknown to be found tvc4 = vc [1:+3,2:+2,1:+9]--bc4 = tmc4 #> tvc4+-- bc4 : right-hand side+-- bc4 = tmc4 #> tvc4  tmc5 = fromListDenseSM 4 $ zipWith (:+) [16..31] [15,14..0] @@ -992,9 +1112,9 @@ aa10 :: SpMatrix Double aa10 = fromListDenseSM 3 [1,2,3,4,5,6] -x10, b10 :: SpVector Double+x10 :: SpVector Double x10 = fromListDenseSV 2 [2,3]-b10 = aa10 #> x10+-- b10 = aa10 #> x10  -- @@ -1035,91 +1155,3 @@ -- m2 = toCSR 4 4 $ V.fromList [(0,0,1), (0,2,5), (2,0,4), (2,3,1), (3,2,2)] -- m3 = toCSR 4 4 $ V.fromList [(1,0,5), (1,1,8), (2,2,3), (3,1,6)] --------- ------- -- run N iterations ---- -- runNBiC :: Int -> SpMatrix Double -> SpVector Double -> BICGSTAB--- runNBiC n aa b = map _xBicgstab $ runAppendN' (bicgstabStep aa x0) n bicgsInit where---    x0 = mkSpVectorD nd $ replicate nd 0.9---    nd = dim r0---    r0 = b ^-^ (aa #> x0)    ---    p0 = r0---    bicgsInit = BICGSTAB x0 r0 p0---- -- runNCGS :: Int -> SpMatrix Double -> SpVector Double -> CGS--- runNCGS n aa b = map _x $ runAppendN' (cgsStep aa x0) n cgsInit where---   x0 = mkSpVectorD nd $ replicate nd 0.1---   nd = dim r0---   r0 = b ^-^ (aa #> x0)    -- residual of initial guess solution---   p0 = r0---   u0 = r0---   cgsInit = CGS x0 r0 p0 u0------ solveRandomN ndim nsp niter = do---   aa0 <- randSpMat ndim (nsp ^ 2)---   let aa = aa0 ^+^ eye ndim---   xtrue <- randSpVec ndim nsp---   let b = aa #> xtrue---       xhatB = head $ runNBiC niter aa b---       xhatC = head $ runNCGS niter aa b---   -- printDenseSM aa    ---   return (normSq (xhatB ^-^ xtrue), normSq (xhatC ^-^ xtrue))----{--random linear system---}------ -- dense--- solveRandom n = do---   aa0 <- randMat n---   let aa = aa0 ^+^ eye n---   xtrue <- randVec n---   -- x0 <- randVec n---   let b = aa #> xtrue---       dx = aa <\> b ^-^ xtrue---   return $ normSq dx---   -- let xhatB = _xBicgstab (bicgstab aa b x0 x0)---   --     xhatC = _x (cgs aa b x0 x0)---   -- return (aa, x, x0, b, xhatB, xhatC)---- -- sparse--- solveSpRandom :: Int -> Int -> IO Double--- solveSpRandom n nsp = do---   aa0 <- randSpMat n nsp---   let aa = aa0 ^+^ eye n---   xtrue <- randSpVec n nsp---   let b = (aa ^+^ eye n) #> xtrue---       dx = aa <\> b ^-^ xtrue---   return $ normSq dx------- solveRandomBanded n bw mu sig = do---   let ndiags = 2*bw---   bands <- replicateM (ndiags + 1) (randArray n mu sig)---   xtrue <- randVec n---   b <- randVec n---   let---     diags = [-bw .. bw - 1]---- randDiagMat :: PrimMonad m =>---      Rows -> Double -> Double -> Int -> m (SpMatrix Double)--- randDiagMat n mu sig i = do---   x <- randArray n mu sig---   return $ mkSubDiagonal n i x