eigen 3.3.4.1 → 3.3.4.2

raw patch · 23 files changed

+2590/−2549 lines, 23 filesdep +constraintsdep +ghc-primdep ~basedep ~binarydep ~bytestringPVP: major bump suggested

API removals or changes: PVP suggests a major version bump

Dependencies added: constraints, ghc-prim

Dependency ranges changed: base, binary, bytestring, primitive, transformers

API changes (from Hackage documentation)

Files

• Data/Eigen/Internal.hsc +0/−247
• Data/Eigen/LA.hs +0/−262
• Data/Eigen/Matrix.hs +0/−647
• Data/Eigen/Matrix/Mutable.hs +0/−116
• Data/Eigen/Parallel.hs +0/−25
• Data/Eigen/SparseLA.hs +0/−439
• Data/Eigen/SparseMatrix.hs +0/−538
• Data/Eigen/SparseMatrix/Mutable.hs +0/−144
• cbits/eigen-sparse-la.cpp +2/−2
• cbits/eigen-sparse.cpp +7/−1
• eigen.cabal +60/−46
• src/Eigen/Internal.hsc +369/−0
• src/Eigen/Matrix.hs +544/−0
• src/Eigen/Matrix/Mutable.hs +154/−0
• src/Eigen/Parallel.hs +28/−0
• src/Eigen/Solver/LA.hs +204/−0
• src/Eigen/Solver/SparseLA.hs +449/−0
• src/Eigen/SparseMatrix.hs +504/−0
• src/Eigen/SparseMatrix/Mutable.hs +163/−0
• test/rank.hs +12/−6
• test/regression.hs +57/−42
• test/solve-sparse.hs +25/−20
• test/solve.hs +12/−14

@@ -1,247 +0,0 @@-{-# OPTIONS_GHC -fno-warn-orphans #-}--{-# LANGUAGE CPP #-} -{-# LANGUAGE EmptyDataDecls  #-}-{-# LANGUAGE FlexibleInstances  #-}-{-# LANGUAGE ForeignFunctionInterface  #-}-{-# LANGUAGE FunctionalDependencies  #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE ScopedTypeVariables  #-}---- | This internal module is going to see a lot of refactoring.---   It is not recommended to import this, as the API is likely---   to experience heavy change.-module Data.Eigen.Internal where--import Control.Monad-import Data.Binary-import Data.Binary.Get-import Data.Binary.Put-import Data.Bits-import Data.Complex-import Foreign.C.String-import Foreign.C.Types-import Foreign.ForeignPtr-import Foreign.Ptr-import Foreign.Storable-import System.IO.Unsafe-import qualified Data.Vector.Storable as VS-import qualified Data.ByteString as BS-import qualified Data.ByteString.Internal as BSI--class (Num a, Cast a b, Cast b a, Storable b, Code b) => Elem a b | a -> b where--instance Elem Float CFloat where-instance Elem Double CDouble where-instance Elem (Complex Float) (CComplex CFloat) where-instance Elem (Complex Double) (CComplex CDouble) where--class Cast a b where-    cast :: a -> b--instance Storable a => Binary (VS.Vector a) where-    put vs = put (BS.length bs) >> putByteString bs where-        (fp,fs) = VS.unsafeToForeignPtr0 vs-        es = sizeOf (VS.head vs)-        bs = BSI.fromForeignPtr (castForeignPtr fp) 0 (fs * es)-        -    get = get >>= getByteString >>= \bs -> let-        (fp,fo,fs) = BSI.toForeignPtr bs-        es = sizeOf (VS.head vs)-        vs = VS.unsafeFromForeignPtr0 (Data.Eigen.Internal.plusForeignPtr fp fo) (fs `div` es)-        in return vs---- | Complex number for FFI with the same memory layout as std::complex\<T\>-data CComplex a = CComplex !a !a deriving Show--instance Storable a => Storable (CComplex a) where-    sizeOf _ = sizeOf (undefined :: a) * 2-    alignment _ = alignment (undefined :: a)-    poke p (CComplex x y) = do-        pokeElemOff (castPtr p) 0 x-        pokeElemOff (castPtr p) 1 y-    peek p = CComplex-        <$> peekElemOff (castPtr p) 0-        <*> peekElemOff (castPtr p) 1--data CTriplet a = CTriplet !CInt !CInt !a deriving Show--instance Storable a => Storable (CTriplet a) where-    sizeOf _ = sizeOf (undefined :: a) + sizeOf (undefined :: CInt) * 2-    alignment _ = alignment (undefined :: CInt)-    poke p (CTriplet row col val) = do-        pokeElemOff (castPtr p) 0 row-        pokeElemOff (castPtr p) 1 col-        pokeByteOff p (sizeOf (undefined :: CInt) * 2) val-    peek p = CTriplet-        <$> peekElemOff (castPtr p) 0-        <*> peekElemOff (castPtr p) 1-        <*> peekByteOff p (sizeOf (undefined :: CInt) * 2)--instance Cast CInt Int where; cast = fromIntegral-instance Cast Int CInt where; cast = fromIntegral-instance Cast CFloat Float where; cast (CFloat x) = x-instance Cast Float CFloat where; cast = CFloat-instance Cast CDouble Double where; cast (CDouble x) = x-instance Cast Double CDouble where; cast = CDouble-instance Cast (CComplex CFloat) (Complex Float) where; cast (CComplex x y) = cast x :+ cast y-instance Cast (Complex Float) (CComplex CFloat) where; cast (x :+ y) = CComplex (cast x) (cast y)-instance Cast (CComplex CDouble) (Complex Double) where; cast (CComplex x y) = cast x :+ cast y-instance Cast (Complex Double) (CComplex CDouble) where; cast (x :+ y) = CComplex (cast x) (cast y)--instance Cast a b => Cast (CTriplet a) (Int, Int, b) where; cast (CTriplet x y z) = (cast x, cast y, cast z)-instance Cast a b => Cast (Int, Int, a) (CTriplet b) where; cast (x,y,z) = CTriplet (cast x) (cast y) (cast z)--intSize :: Int-intSize = sizeOf (undefined :: CInt)--encodeInt :: CInt -> BS.ByteString-encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr--decodeInt :: BS.ByteString -> CInt-decodeInt (BSI.PS fp fo fs)-    | fs == sizeOf x = x-    | otherwise = error "decodeInt: wrong buffer size"-    where x = performIO $ withForeignPtr fp $ peek . (`plusPtr` fo)--data CSparseMatrix a b-type CSparseMatrixPtr a b = Ptr (CSparseMatrix a b)--data CSolver a b-type CSolverPtr a b = Ptr (CSolver a b)--performIO :: IO a -> a-performIO = unsafeDupablePerformIO--plusForeignPtr :: ForeignPtr a -> Int -> ForeignPtr b-plusForeignPtr fp fo = castForeignPtr fp' where-    vs :: VS.Vector CChar-    vs = VS.unsafeFromForeignPtr (castForeignPtr fp) fo 0-    (fp', _) = VS.unsafeToForeignPtr0 vs--foreign import ccall "eigen-proxy.h free" c_freeString :: CString -> IO ()--call :: IO CString -> IO ()-call func = func >>= \c_str -> when (c_str /= nullPtr) $-    peekCString c_str >>= \str -> c_freeString c_str >> fail str--foreign import ccall "eigen-proxy.h free" free :: Ptr a -> IO ()--foreign import ccall "eigen-proxy.h eigen_setNbThreads" c_setNbThreads :: CInt -> IO ()-foreign import ccall "eigen-proxy.h eigen_getNbThreads" c_getNbThreads :: IO CInt--class Code a where; code :: a -> CInt-instance Code CFloat where; code _ = 0-instance Code CDouble where; code _ = 1-instance Code (CComplex CFloat) where; code _ = 2-instance Code (CComplex CDouble) where; code _ = 3--newtype MagicCode = MagicCode CInt deriving Eq--instance Binary MagicCode where-    put (MagicCode code) = putWord32be $ fromIntegral code-    get = MagicCode . fromIntegral <$> getWord32be--magicCode :: Code a => a -> MagicCode-magicCode x = MagicCode (code x `xor` 0x45696730)--#let api1 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name--#api1 random,        "Ptr b -> CInt -> CInt -> IO CString"-#api1 identity,      "Ptr b -> CInt -> CInt -> IO CString"-#api1 add,           "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 sub,           "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 mul,           "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 diagonal,      "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 transpose,     "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 inverse,       "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 adjoint,       "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 conjugate,     "Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 normalize,     "Ptr b -> CInt -> CInt -> IO CString"-#api1 sum,           "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 prod,          "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 mean,          "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 norm,          "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 trace,         "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 squaredNorm,   "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 blueNorm,      "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 hypotNorm,     "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 determinant,   "Ptr b -> Ptr b -> CInt -> CInt -> IO CString"-#api1 rank,          "CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 image,         "CInt -> Ptr (Ptr b) -> Ptr CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 kernel,        "CInt -> Ptr (Ptr b) -> Ptr CInt -> Ptr CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 solve,         "CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"-#api1 relativeError, "Ptr b -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString"--#let api2 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall a b . Code b => %s\n%s = c_%s (code (undefined :: b))", #name, #name, args, #name, args, #name, #name--#api2 sparse_new,           "CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_clone,         "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_fromList,      "CInt -> CInt -> Ptr (CTriplet b) -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_toList,        "CSparseMatrixPtr a b -> Ptr (CTriplet b) -> CInt -> IO CString"-#api2 sparse_free,          "CSparseMatrixPtr a b -> IO CString"-#api2 sparse_makeCompressed,"CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_uncompress,    "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_isCompressed,  "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"-#api2 sparse_transpose,     "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_adjoint,       "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_pruned,        "CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_prunedRef,     "CSparseMatrixPtr a b -> Ptr b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_scale,         "CSparseMatrixPtr a b -> Ptr b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_nonZeros,      "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"-#api2 sparse_innerSize,     "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"-#api2 sparse_outerSize,     "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"-#api2 sparse_coeff,         "CSparseMatrixPtr a b -> CInt -> CInt -> Ptr b -> IO CString"-#api2 sparse_coeffRef,      "CSparseMatrixPtr a b -> CInt -> CInt -> Ptr (Ptr b) -> IO CString"-#api2 sparse_cols,          "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"-#api2 sparse_rows,          "CSparseMatrixPtr a b -> Ptr CInt -> IO CString"-#api2 sparse_norm,          "CSparseMatrixPtr a b -> Ptr b -> IO CString"-#api2 sparse_squaredNorm,   "CSparseMatrixPtr a b -> Ptr b -> IO CString"-#api2 sparse_blueNorm,      "CSparseMatrixPtr a b -> Ptr b -> IO CString"-#api2 sparse_add,           "CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_sub,           "CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_mul,           "CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_block,         "CSparseMatrixPtr a b -> CInt -> CInt -> CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_fromMatrix,    "Ptr b -> CInt -> CInt -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api2 sparse_toMatrix,      "CSparseMatrixPtr a b -> Ptr b -> CInt -> CInt -> IO CString"-#api2 sparse_values,        "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr b) -> IO CString"-#api2 sparse_outerStarts,   "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"-#api2 sparse_innerIndices,  "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"-#api2 sparse_innerNNZs,     "CSparseMatrixPtr a b -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"-#api2 sparse_setZero,       "CSparseMatrixPtr a b -> IO CString"-#api2 sparse_setIdentity,   "CSparseMatrixPtr a b -> IO CString"-#api2 sparse_reserve,       "CSparseMatrixPtr a b -> CInt -> IO CString"-#api2 sparse_resize,        "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString"--#api2 sparse_conservativeResize,    "CSparseMatrixPtr a b -> CInt -> CInt -> IO CString"-#api2 sparse_compressInplace,       "CSparseMatrixPtr a b -> IO CString"-#api2 sparse_uncompressInplace,     "CSparseMatrixPtr a b -> IO CString"---#let api3 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> CInt -> %s\n%s :: forall s a b . (Code s, Code b) => s -> %s\n%s s = c_%s (code (undefined :: b)) (code s)", #name, #name, args, #name, args, #name, #name--#api3 sparse_la_newSolver,          "Ptr (CSolverPtr a b) -> IO CString"-#api3 sparse_la_freeSolver,         "CSolverPtr a b -> IO CString"-#api3 sparse_la_factorize,          "CSolverPtr a b -> CSparseMatrixPtr a b -> IO CString"-#api3 sparse_la_analyzePattern,     "CSolverPtr a b -> CSparseMatrixPtr a b -> IO CString"-#api3 sparse_la_compute,            "CSolverPtr a b -> CSparseMatrixPtr a b -> IO CString"-#api3 sparse_la_tolerance,          "CSolverPtr a b -> Ptr CDouble -> IO CString"-#api3 sparse_la_setTolerance,       "CSolverPtr a b -> CDouble -> IO CString"-#api3 sparse_la_maxIterations,      "CSolverPtr a b -> Ptr CInt -> IO CString"-#api3 sparse_la_setMaxIterations,   "CSolverPtr a b -> CInt -> IO CString"-#api3 sparse_la_info,               "CSolverPtr a b -> Ptr CInt -> IO CString"-#api3 sparse_la_error,              "CSolverPtr a b -> Ptr CDouble -> IO CString"-#api3 sparse_la_iterations,         "CSolverPtr a b -> Ptr CInt -> IO CString"-#api3 sparse_la_solve,              "CSolverPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"--- #api3 sparse_la_solveWithGuess,     "CSolverPtr a b -> CSparseMatrixPtr a b -> CSparseMatrixPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api3 sparse_la_matrixQ,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api3 sparse_la_matrixR,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api3 sparse_la_setPivotThreshold,  "CSolverPtr a b -> CDouble -> IO CString"-#api3 sparse_la_rank,               "CSolverPtr a b -> Ptr CInt -> IO CString"-#api3 sparse_la_matrixL,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api3 sparse_la_matrixU,            "CSolverPtr a b -> Ptr (CSparseMatrixPtr a b) -> IO CString"-#api3 sparse_la_setSymmetric,       "CSolverPtr a b -> CInt -> IO CString"-#api3 sparse_la_determinant,        "CSolverPtr a b -> Ptr b -> IO CString"-#api3 sparse_la_logAbsDeterminant,  "CSolverPtr a b -> Ptr b -> IO CString"-#api3 sparse_la_absDeterminant,     "CSolverPtr a b -> Ptr b -> IO CString"-#api3 sparse_la_signDeterminant,    "CSolverPtr a b -> Ptr b -> IO CString"

@@ -1,262 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE ForeignFunctionInterface #-}-{-# LANGUAGE RecordWildCards #-}--{- |--The problem: You have a system of equations, that you have written as a single matrix equation--@Ax = b@--Where A and b are matrices (b could be a vector, as a special case). You want to find a solution x.--The solution: You can choose between various decompositions, depending on what your matrix A looks like, and depending on whether you favor speed or accuracy. However, let's start with an example that works in all cases, and is a good compromise:--@-import Data.Eigen.Matrix-import Data.Eigen.LA--main = do-    let-        a :: MatrixXd-        a = fromList [[1,2,3], [4,5,6], [7,8,10]]-        b = fromList [[3],[3],[4]]-        x = solve ColPivHouseholderQR a b-    putStrLn \"Here is the matrix A:\" >> print a-    putStrLn \"Here is the vector b:\" >> print b-    putStrLn \"The solution is:\" >> print x-@--produces the following output--@-Here is the matrix A:-Matrix 3x3-1.0 2.0 3.0-4.0 5.0 6.0-7.0 8.0 10.0--Here is the vector b:-Matrix 3x1-3.0-3.0-4.0--The solution is:-Matrix 3x1--2.0000000000000004-1.0000000000000018-0.9999999999999989-@--Checking if a solution really exists: Only you know what error margin you want to allow for a solution to be considered valid.--You can compute relative error using @'norm' (ax - b) / 'norm' b@ formula or use 'relativeError' function which provides the same calculation implemented slightly more efficient.---}--module Data.Eigen.LA (-    -- * Basic linear solving-    Decomposition(..),-    solve,-    relativeError,-    -- * Rank-revealing decompositions-    {- |-Certain decompositions are rank-revealing, i.e. are able to compute the 'rank' of a matrix. These are typically also the decompositions that behave best in the face of a non-full-rank matrix (which in the 'square' case means a singular matrix).--@-import Data.Eigen.Matrix-import Data.Eigen.LA--main = do-    let a = fromList [[1,2,5],[2,1,4],[3,0,3]] :: MatrixXd-    putStrLn "Here is the matrix A:" >> print a-    putStrLn "The rank of A is:" >> print (rank FullPivLU a)-    putStrLn "Here is a matrix whose columns form a basis of the null-space of A:" >> print (kernel FullPivLU a)-    putStrLn "Here is a matrix whose columns form a basis of the column-space of A:" >> print (image FullPivLU a)-@--produces the following output--@-Here is the matrix A:-Matrix 3x3-1.0 2.0 5.0-2.0 1.0 4.0-3.0 0.0 3.0--The rank of A is:-2-Here is a matrix whose columns form a basis of the null-space of A:-Matrix 3x1-0.5000000000000001-1.0--0.5--Here is a matrix whose columns form a basis of the column-space of A:-Matrix 3x2-5.0 1.0-4.0 2.0-3.0 3.0-@-    -}-    rank,-    kernel,-    image,-    -- * Multiple linear regression-    {- | A linear regression model that contains more than one predictor variable. -}-    linearRegression-) where--import Prelude as P-import Foreign.Storable-import Foreign.Marshal.Alloc-import qualified Foreign.Concurrent as FC-#if __GLASGOW_HASKELL__ >= 710-#else-import Control.Applicative-#endif-import Data.Eigen.Matrix-import qualified Data.Eigen.Internal as I-import qualified Data.Eigen.Matrix.Mutable as M-import qualified Data.Vector.Storable as VS--{- |-@-Decomposition           Requirements on the matrix          Speed   Accuracy  Rank  Kernel  Image--PartialPivLU            Invertible                          ++      +         -     -       --FullPivLU               None                                -       +++       +     +       +-HouseholderQR           None                                ++      +         -     -       --ColPivHouseholderQR     None                                +       ++        +     -       --FullPivHouseholderQR    None                                -       +++       +     -       --LLT                     Positive definite                   +++     +         -     -       --LDLT                    Positive or negative semidefinite   +++     ++        -     -       --JacobiSVD               None                                -       +++       +     -       --@-The best way to do least squares solving for square matrices is with a SVD decomposition ('JacobiSVD')--}--data Decomposition-    -- | LU decomposition of a matrix with partial pivoting.-    = PartialPivLU-    -- | LU decomposition of a matrix with complete pivoting.-    | FullPivLU-    -- | Householder QR decomposition of a matrix.-    | HouseholderQR-    -- | Householder rank-revealing QR decomposition of a matrix with column-pivoting.-    | ColPivHouseholderQR-    -- | Householder rank-revealing QR decomposition of a matrix with full pivoting.-    | FullPivHouseholderQR-    -- | Standard Cholesky decomposition (LL^T) of a matrix.-    | LLT-    -- | Robust Cholesky decomposition of a matrix with pivoting.-    | LDLT-    -- | Two-sided Jacobi SVD decomposition of a rectangular matrix.-    | JacobiSVD deriving (Eq, Enum, Show, Read)----- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@-solve :: I.Elem a b => Decomposition -> Matrix a b -> Matrix a b -> Matrix a b-solve d a b = I.performIO $ do-    x <- M.new (cols a) 1-    M.unsafeWith x $ \x_vals x_rows x_cols ->-        unsafeWith a $ \a_vals a_rows a_cols ->-            unsafeWith b $ \b_vals b_rows b_cols ->-                I.call $ I.solve (I.cast $ fromEnum d)-                    x_vals x_rows x_cols-                    a_vals a_rows a_cols-                    b_vals b_rows b_cols-    unsafeFreeze x---- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm-relativeError :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b -> a-relativeError x a b = I.performIO $-    unsafeWith x $ \x_vals x_rows x_cols ->-        unsafeWith a $ \a_vals a_rows a_cols ->-            unsafeWith b $ \b_vals b_rows b_cols ->-                alloca $ \pe -> do-                    I.call $ I.relativeError pe-                        x_vals x_rows x_cols-                        a_vals a_rows a_cols-                        b_vals b_rows b_cols-                    I.cast <$> peek pe---- | The rank of the matrix-rank :: I.Elem a b => Decomposition -> Matrix a b -> Int-rank d m = I.performIO $ alloca $ \pr -> do-    I.call $ unsafeWith m $ I.rank (I.cast $ fromEnum d) pr-    I.cast <$> peek pr---- | Return matrix whose columns form a basis of the null-space of @A@-kernel :: I.Elem a b => Decomposition -> Matrix a b -> Matrix a b-kernel d m1 = I.performIO $-    alloca $ \pvals ->-    alloca $ \prows ->-    alloca $ \pcols ->-        unsafeWith m1 $ \vals1 rows1 cols1 -> do-            I.call $ I.kernel (I.cast $ fromEnum d)-                pvals prows pcols-                vals1 rows1 cols1-            vals <- peek pvals-            rows <- I.cast <$> peek prows-            cols <- I.cast <$> peek pcols-            fp <- FC.newForeignPtr vals $ I.free vals-            return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols----- | Return a matrix whose columns form a basis of the column-space of @A@-image :: I.Elem a b => Decomposition -> Matrix a b -> Matrix a b-image d m1 = I.performIO $-    alloca $ \pvals ->-    alloca $ \prows ->-    alloca $ \pcols ->-        unsafeWith m1 $ \vals1 rows1 cols1 -> do-            I.call $ I.image (I.cast $ fromEnum d)-                pvals prows pcols-                vals1 rows1 cols1-            vals <- peek pvals-            rows <- I.cast <$> peek prows-            cols <- I.cast <$> peek pcols-            fp <- FC.newForeignPtr vals $ I.free vals-            return $ Matrix rows cols $ VS.unsafeFromForeignPtr0 fp $ rows * cols---{- |-[(coeffs, error) = linearRegression points] computes multiple linear regression @y = a1 x1 + a2 x2 + ... + an xn + b@ using 'ColPivHouseholderQR' decomposition--* point format is @[y, x1..xn]@--* coeffs format is @[b, a1..an]@--* error is calculated using 'relativeError'--@-import Data.Eigen.LA-main = print $ linearRegression [-    [-4.32, 3.02, 6.89],-    [-3.79, 2.01, 5.39],-    [-4.01, 2.41, 6.01],-    [-3.86, 2.09, 5.55],-    [-4.10, 2.58, 6.32]]-@-- produces the following output-- @- ([-2.3466569233817127,-0.2534897541434826,-0.1749653335680988],1.8905965120153139e-3)- @---}-linearRegression :: [[Double]] -> ([Double], Double)-linearRegression points = (coeffs, e) where-    a = fromList $ P.map ((1:).tail) points-    b = fromList $ P.map ((:[]).head) points-    x = solve ColPivHouseholderQR a b-    e = relativeError x a b-    coeffs = P.map head $ toList x-----

@@ -1,647 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE DataKinds #-}-{-# LANGUAGE EmptyDataDecls #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE KindSignatures #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE OverloadedStrings #-}-{-# LANGUAGE Rank2Types #-}-{-# LANGUAGE RecordWildCards #-}-{-# LANGUAGE ScopedTypeVariables #-}-{-# LANGUAGE TemplateHaskell #-}-{-# LANGUAGE TypeApplications #-}-{-# LANGUAGE TypeFamilies #-}-{-# LANGUAGE TypeOperators #-}-{-# LANGUAGE QuasiQuotes #-}-{-# LANGUAGE PolyKinds #-}-{-# LANGUAGE EmptyCase #-}--module Data.Eigen.Matrix-  ( -- * Matrix type-    -- | Matrix aliases follows Eigen naming convention-    Matrix(..)-  , MatrixXf-  , MatrixXd-  , MatrixXcf-  , MatrixXcd-  , I.Elem-  , I.CComplex-  , valid-    -- * Matrix conversions-  , fromList-  , toList-  , fromFlatList-  , toFlatList-  , generate-    -- * Standard matrices and special cases-  , empty-  , null-  , square-  , zero-  , ones-  , identity-  , constant-  , random-    -- * Accessing matrix data-  , cols-  , rows-  , dims-  , (!)-  , coeff-  , unsafeCoeff-  , col-  , row-  , block-  , topRows-  , bottomRows-  , leftCols-  , rightCols-    -- * Matrix properties-  , sum-  , prod-  , mean-  , minCoeff-  , maxCoeff-  , trace-  , norm-  , squaredNorm-  , blueNorm-  , hypotNorm-  , determinant-    -- * Generic reductions-  , fold-  , fold'-  , ifold-  , ifold'-  , fold1-  , fold1'-    -- * Boolean reductions-  , all-  , any-  , count-    -- * Basic matrix algebra-  , add-  , sub-  , mul-    -- * Mapping over elements-  , map-  , imap-  , filter-  , ifilter-    -- * Matrix transformations-  , diagonal-  , transpose-  , inverse-  , adjoint-  , conjugate-  , normalize-  , modify-  , convert-  , TriangularMode(..)-  , triangularView-  , lowerTriangle-  , upperTriangle-    -- * Matrix serialization-  , encode-  , decode-    -- * Mutable matrices-  , thaw-  , freeze-  , unsafeThaw-  , unsafeFreeze-    -- * Raw pointers-  , unsafeWith -) where--import Control.Monad-import Control.Monad.Primitive-import Control.Monad.ST-import Data.Binary hiding (encode, decode)-import Data.Complex hiding (conjugate)-import Data.Tuple-import Foreign.C.String-import Foreign.C.Types-import Foreign.Marshal.Alloc-import Foreign.Ptr-import Foreign.Storable-import Prelude hiding (null, sum, all, any, map, filter)-import Text.Printf-import qualified Data.Binary as B-import qualified Data.ByteString.Lazy as BSL-import qualified Data.Eigen.Internal as I-import qualified Data.Eigen.Matrix.Mutable as M-import qualified Data.List as L-import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Storable.Mutable as VSM-import qualified Prelude as P---- | Matrix to be used in pure computations, uses column major memory layout, features copy-free FFI with C++ <http://eigen.tuxfamily.org Eigen> library.-data Matrix a b where-    Matrix :: I.Elem a b => !Int -> !Int -> !(VS.Vector b) -> Matrix a b---- | Alias for single precision matrix-type MatrixXf = Matrix Float CFloat--- | Alias for double precision matrix-type MatrixXd = Matrix Double CDouble--- | Alias for single previsiom matrix of complex numbers-type MatrixXcf = Matrix (Complex Float) (I.CComplex CFloat)--- | Alias for double prevision matrix of complex numbers-type MatrixXcd = Matrix (Complex Double) (I.CComplex CDouble)---- | Pretty prints the matrix-instance (I.Elem a b, Show a) => Show (Matrix a b) where-    show m@(Matrix rows cols _) = concat [-        "Matrix ", show rows, "x", show cols,-        "\n", L.intercalate "\n" $ P.map (L.intercalate "\t" . P.map show) $ toList m, "\n"]---- | Basic matrix math exposed through Num instance: @(*)@, @(+)@, @(-)@, `fromInteger`, `signum`, `abs`, `negate`-instance I.Elem a b => Num (Matrix a b) where-    (*) = mul-    (+) = add-    (-) = sub-    fromInteger = constant 1 1 . fromInteger-    signum = map signum-    abs = map abs-    negate = map negate---- | Matrix binary serialization-instance I.Elem a b => Binary (Matrix a b) where-    put (Matrix rows cols vals) = do-        put $ I.magicCode (undefined :: b)-        put rows-        put cols-        put vals--    get = do-        get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))-        Matrix <$> get <*> get <*> get---- | Encode the matrix as a lazy byte string-encode :: I.Elem a b => Matrix a b -> BSL.ByteString-encode = B.encode---- | Decode matrix from the lazy byte string-decode :: I.Elem a b => BSL.ByteString -> Matrix a b-decode = B.decode---- | Empty 0x0 matrix-{-# INLINE empty #-}-empty :: I.Elem a b => Matrix a b-empty = Matrix 0 0 VS.empty---- | Is matrix empty?-{-# INLINE null #-}-null :: I.Elem a b => Matrix a b -> Bool-null (Matrix rows cols _) = rows == 0 && cols == 0---- | Is matrix square?-{-# INLINE square #-}-square :: I.Elem a b => Matrix a b -> Bool-square (Matrix rows cols _) = rows == cols---- | Matrix where all coeffs are filled with given value-{-# INLINE constant #-}-constant :: I.Elem a b => Int -> Int -> a -> Matrix a b-constant rows cols val = Matrix rows cols $ VS.replicate (rows * cols) (I.cast val)---- | Matrix where all coeff are 0-{-# INLINE zero #-}-zero :: I.Elem a b => Int -> Int -> Matrix a b-zero rows cols = constant rows cols 0---- | Matrix where all coeff are 1-{-# INLINE ones #-}-ones :: I.Elem a b => Int -> Int -> Matrix a b-ones rows cols = constant rows cols 1---- | The identity matrix (not necessarily square).-identity :: I.Elem a b => Int -> Int -> Matrix a b-identity rows cols = I.performIO $ do-    m <- M.new rows cols-    I.call $ M.unsafeWith m I.identity-    unsafeFreeze m---- | The random matrix of a given size-random :: I.Elem a b => Int -> Int -> IO (Matrix a b)-random rows cols = do-    m <- M.new rows cols-    I.call $ M.unsafeWith m I.random-    unsafeFreeze m---- | Number of rows for the matrix-{-# INLINE rows #-}-rows :: I.Elem a b => Matrix a b -> Int-rows (Matrix rows _ _) = rows---- | Number of columns for the matrix-{-# INLINE cols #-}-cols :: I.Elem a b => Matrix a b -> Int-cols (Matrix _ cols _) = cols---- | Mtrix size as (rows, cols) pair-{-# INLINE dims #-}-dims :: I.Elem a b => Matrix a b -> (Int, Int)-dims (Matrix rows cols _) = (rows, cols)---- | Matrix coefficient at specific row and col-{-# INLINE (!) #-}-(!) :: forall a b. (I.Elem a b) => Matrix a b -> (Int, Int) -> a-(!) m (row,col) = coeff row col m---- | Matrix coefficient at specific row and col-{-# INLINE coeff #-}-coeff :: I.Elem a b => Int -> Int -> Matrix a b -> a-coeff row col m@(Matrix rows cols _)-    | not (valid m) = error "matrix is not valid"-    | row < 0 || row >= rows = error $ printf "Matrix.coeff: row %d is out of bounds [0..%d)" row rows-    | col < 0 || col >= cols = error $ printf "Matrix.coeff: col %d is out of bounds [0..%d)" col cols-    | otherwise = unsafeCoeff row col m---- | Unsafe version of coeff function. No bounds check performed so SEGFAULT possible-{-# INLINE unsafeCoeff #-}-unsafeCoeff :: I.Elem a b => Int -> Int -> Matrix a b -> a-unsafeCoeff row col (Matrix rows _ vals) = I.cast $ VS.unsafeIndex vals $ col * rows + row---- | List of coefficients for the given col-{-# INLINE col #-}-col :: I.Elem a b => Int -> Matrix a b -> [a]-col c m@(Matrix rows _ _) = [coeff r c m | r <- [0..pred rows]]---- | List of coefficients for the given row-{-# INLINE row #-}-row :: I.Elem a b => Int -> Matrix a b -> [a]-row r m@(Matrix _ cols _) = [coeff r c m | c <- [0..pred cols]]---- | Extract rectangular block from matrix defined by startRow startCol blockRows blockCols-block :: I.Elem a b => Int -> Int -> Int -> Int -> Matrix a b -> Matrix a b-block startRow startCol blockRows blockCols m =-    generate blockRows blockCols $ \row col ->-        coeff (startRow + row) (startCol + col) m----valid :: I.Elem a b => Matrix a b -> Bool---valid (Matrix rows cols vals) =---  let rowsGood = rows >= 0---      colsGood = cols >= 0---      valsGood = VS.length vals == rows * cols---  in rowGood && colsGood && valsGood-   --- | Verify matrix dimensions and memory layout-{-# INLINE valid #-}-valid :: I.Elem a b => Matrix a b -> Bool-valid (Matrix rows cols vals) =-     rows >= 0-  && cols >= 0-  && VS.length vals == rows * cols---- | The maximum coefficient of the matrix-{-# INLINE maxCoeff #-}-maxCoeff :: (I.Elem a b, Ord a) => Matrix a b -> a-maxCoeff = fold1' max---- | The minimum coefficient of the matrix-{-# INLINE minCoeff #-}-minCoeff :: (I.Elem a b, Ord a) => Matrix a b -> a-minCoeff = fold1' min---- | Top @N@ rows of matrix-{-# INLINE topRows #-}-topRows :: I.Elem a b => Int -> Matrix a b -> Matrix a b-topRows n m@(Matrix _ cols _) = block 0 0 n cols m---- | Bottom @N@ rows of matrix-{-# INLINE bottomRows #-}-bottomRows :: I.Elem a b => Int -> Matrix a b -> Matrix a b-bottomRows n m@(Matrix rows cols _) = block (rows - n) 0 n cols m---- | Left @N@ columns of matrix-{-# INLINE leftCols #-}-leftCols :: I.Elem a b => Int -> Matrix a b -> Matrix a b-leftCols n m@(Matrix rows _ _) = block 0 0 rows n m---- | Right @N@ columns of matrix-{-# INLINE rightCols #-}-rightCols :: I.Elem a b => Int -> Matrix a b -> Matrix a b-rightCols n m@(Matrix rows cols _) = block 0 (cols - n) rows n m---- | Construct matrix from a list of rows, column count is detected as maximum row length. Missing values are filled with 0-fromList :: I.Elem a b => [[a]] -> Matrix a b-fromList list = Matrix rows cols vals where-    rows = length list-    cols = L.foldl' max 0 $ P.map length list-    vals = VS.create $ do-        vm <- VSM.replicate (rows * cols) (I.cast (0 `asTypeOf` (head (head list))))-        forM_ (zip [0..] list) $ \(row, vals) ->-            forM_ (zip [0..] vals) $ \(col, val) ->-                VSM.write vm (col * rows + row) (I.cast val)-        return vm---- | Convert matrix to a list of rows-toList :: I.Elem a b => Matrix a b -> [[a]]-toList m@(Matrix rows cols vals)-    | not (valid m) = error "matrix is not valid"-    | otherwise = [[I.cast $ vals `VS.unsafeIndex` (col * rows + row) | col <- [0..pred cols]] | row <- [0..pred rows]]---- | Build matrix of given dimensions and values from given list split on rows. Invalid list length results in error.-fromFlatList :: I.Elem a b => Int -> Int -> [a] -> Matrix a b-fromFlatList rows cols list-    | not (rows * cols == (length list)) = error $ concat ["cannot construct ", show rows, "x", show cols, " matrix from ", show $ length list, " values"]-    | otherwise = Matrix rows cols vals where-        vals = VS.create $ do-            vm <- VSM.replicate (rows * cols) (I.cast (0 `asTypeOf` (head list)))-            forM_ (zip [(col * rows + row) | row <- [0..pred rows], col <- [0..pred cols]] list) $ \(idx, val) ->-                VSM.write vm idx (I.cast val)-            return vm---- | Convert matrix to a list by concatenating rows-toFlatList :: I.Elem a b => Matrix a b -> [a]-toFlatList m@(Matrix rows cols vals)-    | not (valid m) = error "matrix is not valid"-    | otherwise = [I.cast $ vals `VS.unsafeIndex` (col * rows + row) | row <- [0..pred rows], col <- [0..pred cols]]---- | [generate rows cols (λ row col -> val)]------ Create matrix using generator function @λ row col -> val@----generate :: I.Elem a b => Int -> Int -> (Int -> Int -> a) -> Matrix a b-generate rows cols f = Matrix rows cols $ VS.create $ do-    vals <- VSM.new (rows * cols)-    forM_ [0..pred rows] $ \row ->-        forM_ [0..pred cols] $ \col ->-            VSM.write vals (col * rows + row) (I.cast $ f row col)-    return vals---- | The sum of all coefficients of the matrix-sum :: I.Elem a b => Matrix a b -> a-sum = _prop I.sum---- | The product of all coefficients of the matrix-prod :: I.Elem a b => Matrix a b -> a-prod = _prop I.prod---- | The mean of all coefficients of the matrix-mean :: I.Elem a b => Matrix a b -> a-mean = _prop I.mean---- | The trace of a matrix is the sum of the diagonal coefficients and can also be computed as sum (diagonal m)-trace :: I.Elem a b => Matrix a b -> a-trace = _prop I.trace---- | Applied to a predicate and a matrix, all determines if all elements of the matrix satisfies the predicate-all :: I.Elem a b => (a -> Bool) -> Matrix a b -> Bool-all f = VS.all (f . I.cast) . _vals---- | Applied to a predicate and a matrix, any determines if any element of the matrix satisfies the predicate-any :: I.Elem a b => (a -> Bool) -> Matrix a b -> Bool-any f = VS.any (f . I.cast) . _vals---- | Returns the number of coefficients in a given matrix that evaluate to true-count :: I.Elem a b => (a -> Bool) -> Matrix a b -> Int-count f = VS.foldl' (\n x -> if f (I.cast x) then succ n else n) 0 . _vals--{-| For vectors, the l2 norm, and for matrices the Frobenius norm.-    In both cases, it consists in the square root of the sum of the square of all the matrix entries.-    For vectors, this is also equals to the square root of the dot product of this with itself.--}-norm :: I.Elem a b => Matrix a b -> a-norm = _prop I.norm---- | For vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.-squaredNorm :: I.Elem a b => Matrix a b -> a-squaredNorm = _prop I.squaredNorm---- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.-blueNorm :: I.Elem a b => Matrix a b -> a-blueNorm = _prop I.blueNorm---- | The l2 norm of the matrix avoiding undeflow and overflow. This version use a concatenation of hypot calls, and it is very slow.-hypotNorm :: I.Elem a b => Matrix a b -> a-hypotNorm = _prop I.hypotNorm---- | The determinant of the matrix-determinant :: I.Elem a b => Matrix a b -> a-determinant m -- = _prop I.determinant (unrefine m)-    | square m = _prop I.determinant m-    | otherwise = error "Matrix.determinant: non-square matrix"---- | Adding two matrices by adding the corresponding entries together. You can use @(+)@ function as well.-add :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b-add m1 m2-    | dims m1 == dims m2 = _binop const I.add m1 m2-    | otherwise = error "Matrix.add: matrices should have the same size"---- | Subtracting two matrices by subtracting the corresponding entries together. You can use @(-)@ function as well.-sub :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b-sub m1 m2-    | dims m1 == dims m2 = _binop const I.sub m1 m2-    | otherwise = error "Matrix.add: matrices should have the same size"---- | Matrix multiplication. You can use @(*)@ function as well.-mul :: I.Elem a b => Matrix a b -> Matrix a b -> Matrix a b-mul m1 m2-    | cols m1 == rows m2 = _binop (\(rows, _) (_, cols) -> (rows, cols)) I.mul m1 m2-    | otherwise = error "Matrix.mul: number of columns for lhs matrix should be the same as number of rows for rhs matrix"--{- | Apply a given function to each element of the matrix.--Here is an example how to implement scalar matrix multiplication:-->>> let a = fromList [[1,2],[3,4]] :: MatrixXf-->>> a-Matrix 2x2-1.0 2.0-3.0 4.0-->>> map (*10) a-Matrix 2x2-10.0    20.0-30.0    40.0---}-map :: I.Elem a b => (a -> a) -> Matrix a b -> Matrix a b-map f (Matrix rows cols vals) = Matrix rows cols (VS.map (I.cast . f . I.cast) vals)--{- | Apply a given function to each element of the matrix.--Here is an example how upper triangular matrix can be implemented:-->>> let a = fromList [[1,2,3],[4,5,6],[7,8,9]] :: MatrixXf-->>> a-Matrix 3x3-1.0 2.0 3.0-4.0 5.0 6.0-7.0 8.0 9.0-->>> imap (\row col val -> if row <= col then val else 0) a-Matrix 3x3-1.0 2.0 3.0-0.0 5.0 6.0-0.0 0.0 9.0---}--imap :: I.Elem a b => (Int -> Int -> a -> a) -> Matrix a b -> Matrix a b-imap f (Matrix rows cols vals) = Matrix rows cols (VS.imap (\n -> let (c, r) = divMod n rows in I.cast . f r c . I.cast) vals)--data TriangularMode-    -- | View matrix as a lower triangular matrix.-    = Lower-    -- | View matrix as an upper triangular matrix.-    | Upper-    -- | View matrix as a lower triangular matrix with zeros on the diagonal.-    | StrictlyLower-    -- | View matrix as an upper triangular matrix with zeros on the diagonal.-    | StrictlyUpper-    -- | View matrix as a lower triangular matrix with ones on the diagonal.-    | UnitLower-    -- | View matrix as an upper triangular matrix with ones on the diagonal.-    | UnitUpper deriving (Eq, Enum, Show, Read)---- | Triangular view extracted from the current matrix-triangularView :: I.Elem a b => TriangularMode -> Matrix a b -> Matrix a b-triangularView Lower         = imap $ \row col val -> case compare row col of { LT -> 0; _ -> val }-triangularView Upper         = imap $ \row col val -> case compare row col of { GT -> 0; _ -> val }-triangularView StrictlyLower = imap $ \row col val -> case compare row col of { GT -> val; _ -> 0 }-triangularView StrictlyUpper = imap $ \row col val -> case compare row col of { LT -> val; _ -> 0 }-triangularView UnitLower     = imap $ \row col val -> case compare row col of { GT -> val; LT -> 0; EQ -> 1 }-triangularView UnitUpper     = imap $ \row col val -> case compare row col of { LT -> val; GT -> 0; EQ -> 1 }---- | Lower trinagle of the matrix. Shortcut for @triangularView Lower@-lowerTriangle :: I.Elem a b => Matrix a b -> Matrix a b-lowerTriangle = triangularView Lower---- | Upper trinagle of the matrix. Shortcut for @triangularView Upper@-upperTriangle :: I.Elem a b => Matrix a b -> Matrix a b-upperTriangle = triangularView Upper---- | Filter elements in the matrix. Filtered elements will be replaced by 0-filter :: I.Elem a b => (a -> Bool) -> Matrix a b -> Matrix a b-filter f = map (\x -> if f x then x else 0)---- | Filter elements in the matrix. Filtered elements will be replaced by 0-ifilter :: I.Elem a b => (Int -> Int -> a -> Bool) -> Matrix a b -> Matrix a b-ifilter f = imap (\r c x -> if f r c x then x else 0)---- | Reduce matrix using user provided function applied to each element.-fold :: I.Elem a b => (c -> a -> c) -> c -> Matrix a b -> c-fold f a (Matrix _ _ vals) = VS.foldl (\a x -> f a (I.cast x)) a vals---- | Reduce matrix using user provided function applied to each element. This is strict version of 'fold'-fold' :: I.Elem a b => (c -> a -> c) -> c -> Matrix a b -> c-fold' f a (Matrix _ _ vals) = VS.foldl' (\a x -> f a (I.cast x)) a vals---- | Reduce matrix using user provided function applied to each element and it's index-ifold :: I.Elem a b => (Int -> Int -> c -> a -> c) -> c -> Matrix a b -> c-ifold f a (Matrix rows _ vals) = VS.ifoldl (\a n x -> let (c,r) = divMod n rows in f r c a (I.cast x)) a vals---- | Reduce matrix using user provided function applied to each element and it's index. This is strict version of 'ifold'-ifold' :: I.Elem a b => (Int -> Int -> c -> a -> c) -> c -> Matrix a b -> c-ifold' f a (Matrix rows _ vals) = VS.ifoldl' (\a n x -> let (c,r) = divMod n rows in f r c a (I.cast x)) a vals---- | Reduce matrix using user provided function applied to each element.-fold1 :: I.Elem a b => (a -> a -> a) -> Matrix a b -> a-fold1 f = foldl1 f . P.map I.cast . VS.toList . _vals---- | Reduce matrix using user provided function applied to each element. This is strict version of 'fold'-fold1' :: I.Elem a b => (a -> a -> a) -> Matrix a b -> a-fold1' f = L.foldl1' f . P.map I.cast . VS.toList . _vals---- | Diagonal of the matrix-diagonal :: I.Elem a b => Matrix a b -> Matrix a b-diagonal = _unop (\(rows, cols) -> (min rows cols, 1)) I.diagonal--{- | Inverse of the matrix--For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses PartialPivLU decomposition--}-inverse :: I.Elem a b => Matrix a b -> Matrix a b-inverse m -- = _unop id I.inverse (unrefine m)-    | square m = _unop id I.inverse m-    | otherwise = error "Matrix.inverse: non-square matrix"---- | Adjoint of the matrix-adjoint :: I.Elem a b => Matrix a b -> Matrix a b-adjoint = _unop swap I.adjoint---- | Transpose of the matrix-transpose :: I.Elem a b => Matrix a b -> Matrix a b-transpose = _unop swap I.transpose---- | Conjugate of the matrix-conjugate :: I.Elem a b => Matrix a b -> Matrix a b-conjugate = _unop id I.conjugate---- | Nomalize the matrix by deviding it on its 'norm'-normalize :: I.Elem a b => Matrix a b -> Matrix a b-normalize (Matrix rows cols vals) = I.performIO $ do-    vals <- VS.thaw vals-    VSM.unsafeWith vals $ \p ->-        I.call $ I.normalize p (I.cast rows) (I.cast cols)-    Matrix rows cols <$> VS.unsafeFreeze vals---- | Apply a destructive operation to a matrix. The operation will be performed in place if it is safe to do so and will modify a copy of the matrix otherwise.-modify :: I.Elem a b => (forall s. M.MMatrix a b s -> ST s ()) -> Matrix a b -> Matrix a b-modify f (Matrix rows cols vals) = Matrix rows cols (VS.modify (f . M.MMatrix rows cols) vals)---- | Convert matrix to different type using user provided element converter-convert :: (I.Elem a b, I.Elem c d) => (a -> c) -> Matrix a b -> Matrix c d-convert f (Matrix rows cols vals) = Matrix rows cols $ VS.map (I.cast . f . I.cast) vals---- | Yield an immutable copy of the mutable matrix-freeze :: I.Elem a b => PrimMonad m => M.MMatrix a b (PrimState m) -> m (Matrix a b)-freeze (M.MMatrix mrows mcols mvals) = VS.freeze mvals >>= return . Matrix mrows mcols---- | Yield a mutable copy of the immutable matrix-thaw :: I.Elem a b => PrimMonad m => Matrix a b -> m (M.MMatrix a b (PrimState m))-thaw (Matrix rows cols vals) = VS.thaw vals >>= return . M.MMatrix rows cols---- | Unsafe convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.-unsafeFreeze :: I.Elem a b => PrimMonad m => M.MMatrix a b (PrimState m) -> m (Matrix a b)-unsafeFreeze (M.MMatrix mrows mcols mvals) = VS.unsafeFreeze mvals >>= return . Matrix mrows mcols---- | Unsafely convert an immutable matrix to a mutable one without copying. The immutable matrix may not be used after this operation.-unsafeThaw :: I.Elem a b => PrimMonad m => Matrix a b -> m (M.MMatrix a b (PrimState m))-unsafeThaw (Matrix rows cols vals) = VS.unsafeThaw vals >>= return . M.MMatrix rows cols---- | Pass a pointer to the matrix's data to the IO action. The data may not be modified through the pointer.-unsafeWith :: I.Elem a b => Matrix a b -> (Ptr b -> CInt -> CInt -> IO c) -> IO c-unsafeWith m@(Matrix rows cols vals) f-    | not (valid m) = fail "Matrix.unsafeWith: matrix layout is invalid"-    | otherwise = VS.unsafeWith vals $ \p -> f p (I.cast rows) (I.cast cols)--{-# INLINE _prop #-}-_prop :: I.Elem a b => (Ptr b -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> a-_prop f m = I.cast $ I.performIO $ alloca $ \p -> do-    I.call $ unsafeWith m (f p)-    peek p--{-# INLINE _binop #-}-_binop :: I.Elem a b => ((Int, Int) -> (Int, Int) -> (Int, Int)) -> (Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> Matrix a b -> Matrix a b-_binop f g m1 m2 = I.performIO $ do-    m0 <- uncurry M.new $ f (dims m1) (dims m2)-    M.unsafeWith m0 $ \vals0 rows0 cols0 ->-        unsafeWith m1 $ \vals1 rows1 cols1 ->-            unsafeWith m2 $ \vals2 rows2 cols2 ->-                I.call $ g-                    vals0 rows0 cols0-                    vals1 rows1 cols1-                    vals2 rows2 cols2-    unsafeFreeze m0--{-# INLINE _unop #-}-_unop :: I.Elem a b => ((Int,Int) -> (Int,Int)) -> (Ptr b -> CInt -> CInt -> Ptr b -> CInt -> CInt -> IO CString) -> Matrix a b -> Matrix a b-_unop f g m1 = I.performIO $ do-    m0 <- uncurry M.new $ f (dims m1)-    M.unsafeWith m0 $ \vals0 rows0 cols0 ->-        unsafeWith m1 $ \vals1 rows1 cols1 ->-            I.call $ g-                vals0 rows0 cols0-                vals1 rows1 cols1-    unsafeFreeze m0--{-# INLINE _vals #-}-_vals :: I.Elem a b => Matrix a b -> VS.Vector b-_vals (Matrix _ _ vals) = vals

@@ -1,116 +0,0 @@-{-# LANGUAGE RecordWildCards, ScopedTypeVariables #-}-module Data.Eigen.Matrix.Mutable (-    MMatrix(..),-    MMatrixXf,-    MMatrixXd,-    MMatrixXcf,-    MMatrixXcd,-    IOMatrix,-    STMatrix,-    -- * Construction-    new,-    replicate,-    -- * Consistency check-    valid,-    -- * Accessing individual elements-    read,-    write,-    unsafeRead,-    unsafeWrite,-    -- * Modifying matrices-    set,-    copy,-    unsafeCopy,-    -- * Raw pointers-    unsafeWith-) where--import Prelude hiding (read, replicate)-import Control.Monad.Primitive-import Foreign.Ptr-import Foreign.C.Types-import Data.Complex-import Text.Printf-import qualified Data.Vector.Storable.Mutable as VSM-import qualified Data.Eigen.Internal as I---- | Mutable matrix. You can modify elements-data MMatrix a b s = MMatrix {-    mm_rows :: Int,-    mm_cols :: Int,-    mm_vals :: VSM.MVector s b-}---- | Alias for single precision mutable matrix-type MMatrixXf = MMatrix Float CFloat--- | Alias for double precision mutable matrix-type MMatrixXd = MMatrix Double CDouble--- | Alias for single previsiom mutable matrix of complex numbers-type MMatrixXcf = MMatrix (Complex Float) (I.CComplex CFloat)--- | Alias for double prevision mutable matrix of complex numbers-type MMatrixXcd = MMatrix (Complex Double) (I.CComplex CDouble)--type IOMatrix a b = MMatrix a b RealWorld-type STMatrix a b s = MMatrix a b s---- | Verify matrix dimensions and memory layout-valid :: I.Elem a b => MMatrix a b s -> Bool-valid MMatrix{..} = mm_rows >= 0 && mm_cols >= 0 && VSM.length mm_vals == mm_rows * mm_cols---- | Create a mutable matrix of the given size and fill it with 0 as an initial value.-new :: (PrimMonad m, I.Elem a b) => Int -> Int -> m (MMatrix a b (PrimState m))-new rows cols = replicate rows cols 0---- | Create a mutable matrix of the given size and fill it with as an initial value.-replicate :: (PrimMonad m, I.Elem a b) => Int -> Int -> a -> m (MMatrix a b (PrimState m))-replicate rows cols val = do-    vals <- VSM.replicate (rows * cols) (I.cast val)-    return $ MMatrix rows cols vals---- | Set all elements of the matrix to the given value-set :: (PrimMonad m, I.Elem a b) => (MMatrix a b (PrimState m)) -> a -> m ()-set MMatrix{..} val = VSM.set mm_vals (I.cast val)---- | Copy a matrix. The two matrices must have the same size and may not overlap.-copy :: (PrimMonad m, I.Elem a b) => (MMatrix a b (PrimState m)) -> (MMatrix a b (PrimState m)) -> m ()-copy m1 m2-    | not (valid m1) = fail "MMatrix.copy: lhs matrix layout is invalid"-    | not (valid m2) = fail "MMatrix.copy: rhs matrix layout is invalid"-    | mm_rows m1 /= mm_rows m2 = fail "MMatrix.copy: matrices have different number of cols"-    | mm_cols m1 /= mm_cols m2 = fail "MMatrix.copy: matrices have different number of rows"-    | otherwise = VSM.copy (mm_vals m1) (mm_vals m2)---- | Yield the element at the given position.-read :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> m a-read mm@MMatrix{..} row col-    | not (valid mm) = fail "MMatrix.read: matrix layout is invalid"-    | row < 0 || row >= mm_rows = fail $ printf "MMatrix.read: row %d is out of bounds [0..%d)" row mm_rows-    | col < 0 || col >= mm_cols = fail $ printf "MMatrix.read: col %d is out of bounds [0..%d)" col mm_cols-    | otherwise = unsafeRead mm row col---- | Replace the element at the given position.-write :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> a -> m ()-write mm@MMatrix{..} row col val-    | not (valid mm) = fail "MMatrix.write: matrix layout is invalid"-    | row < 0 || row >= mm_rows = fail $ printf "MMatrix.write: row %d is out of bounds [0..%d)" row mm_rows-    | col < 0 || col >= mm_cols = fail $ printf "MMatrix.write: col %d is out of bounds [0..%d)" col mm_cols-    | otherwise = unsafeWrite mm row col val---- | Copy a matrix. The two matrices must have the same size and may not overlap however no bounds check performaned to it may SEGFAULT for incorrect input.-unsafeCopy :: (PrimMonad m, I.Elem a b) => (MMatrix a b (PrimState m)) -> (MMatrix a b (PrimState m)) -> m ()-unsafeCopy m1 m2 = VSM.unsafeCopy (mm_vals m1) (mm_vals m2)---- | Yield the element at the given position. No bounds checks are performed.-unsafeRead :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> m a-unsafeRead MMatrix{..} row col = VSM.unsafeRead mm_vals (col * mm_rows + row) >>= \val -> return (I.cast val)---- | Replace the element at the given position. No bounds checks are performed.-unsafeWrite :: (PrimMonad m, I.Elem a b) => MMatrix a b (PrimState m) -> Int -> Int -> a -> m ()-unsafeWrite MMatrix{..} row col val = VSM.unsafeWrite mm_vals (col * mm_rows + row) (I.cast val)---- | Pass a pointer to the matrix's data to the IO action. Modifying data through the pointer is unsafe if the matrix could have been frozen before the modification.-unsafeWith :: I.Elem a b => IOMatrix a b -> (Ptr b -> CInt -> CInt -> IO c) -> IO c-unsafeWith mm@MMatrix{..} f-    | not (valid mm) = fail "mutable matrix layout is invalid"-    | otherwise = VSM.unsafeWith mm_vals $ \p -> f p (I.cast mm_rows) (I.cast mm_cols)-

@@ -1,25 +0,0 @@-{- |-Some Eigen's algorithms can exploit the multiple cores present in your hardware.-To this end, it is enough to enable OpenMP on your compiler, for instance: GCC: -fopenmp.-You can control the number of thread that will be used using either the OpenMP API or Eiegn's API using the following priority:--1. OMP_NUM_THREADS=n ./my_program-2. setNbThreads n--Unless setNbThreads has been called, Eigen uses the number of threads specified by OpenMP.-You can restore this behaviour by calling @setNbThreads n@--Currently, the following algorithms can make use of multi-threading: general matrix - matrix products PartialPivLU.--}--module Data.Eigen.Parallel where--import Data.Eigen.Internal---- | Sets the max number of threads reserved for Eigen-setNbThreads :: Int -> IO ()-setNbThreads = c_setNbThreads . cast---- | Gets the max number of threads reserved for Eigen-getNbThreads :: IO Int-getNbThreads = fmap cast $ c_getNbThreads

@@ -1,439 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE FlexibleInstances #-}-{-# LANGUAGE ForeignFunctionInterface #-}-{-# LANGUAGE FunctionalDependencies #-}-{-# LANGUAGE MultiParamTypeClasses #-}-{-# LANGUAGE RecordWildCards #-}-{-# LANGUAGE GeneralizedNewtypeDeriving #-}--{- |--This documentation is based on original Eigen page <http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html Solving Sparse Linear Systems>--Eigen currently provides a limited set of built-in MPL2 compatible solvers.-They are summarized in the following table:--@-Sparse solver       Solver kind             Matrix kind         Notes--ConjugateGradient   Classic iterative CG    SPD                 Recommended for large symmetric-                                                                problems (e.g., 3D Poisson eq.)-BiCGSTAB            Iterative stabilized    Square-                    bi-conjugate gradient-SparseLU            LU factorization        Square              Optimized for small and large problems-                                                                with irregular patterns-SparseQR            QR factorization        Any, rectangular    Recommended for least-square problems,-                                                                has a basic rank-revealing feature-@--All these solvers follow the same general concept. Here is a typical and general example:--@-let-    a :: SparseMatrixXd-    a = ... -- fill a--    b :: SparseMatrixXd-    b = ... -- fill b--    validate msg = info >>= (`when` fail msg) . (/= Success)--// solve Ax = b-runSolverT solver $ do-    compute a-    validate "decomposition failed"--    x <- solve b-    validate "solving failed"--    // solve for another right hand side-    x1 <- solve b1-@--In the case where multiple problems with the same sparsity pattern have to be solved, then the "compute" step can be decomposed as follow:--@-runSolverT solver $ do-    analyzePattern a1-    factorize a1-    x1 <- solve b1-    x2 <- solve b2--    factorize a2-    x1 <- solve b1-    x2 <- solve b2-@--Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on.---}--module Data.Eigen.SparseLA (-    -- * Sparse Solvers-    Solver,-    DirectSolver,-    IterativeSolver,-    OrderingMethod(..),-    Preconditioner(..),-    ConjugateGradient(ConjugateGradient),-    BiCGSTAB(BiCGSTAB),-    SparseLU(SparseLU),-    SparseQR(SparseQR),-    ComputationInfo(..),-    SolverT,-    runSolverT,-    -- * The Compute step-    {- |-        In the `compute` function, the matrix is generally factorized: LLT for self-adjoint matrices, LDLT for general hermitian matrices,-        LU for non hermitian matrices and QR for rectangular matrices. These are the results of using direct solvers.-        For this class of solvers precisely, the compute step is further subdivided into `analyzePattern` and `factorize`.--        The goal of `analyzePattern` is to reorder the nonzero elements of the matrix, such that the factorization step creates less fill-in.-        This step exploits only the structure of the matrix. Hence, the results of this step can be used for other linear systems where the-        matrix has the same structure.--        In `factorize`, the factors of the coefficient matrix are computed. This step should be called each time the values of the matrix change.-        However, the structural pattern of the matrix should not change between multiple calls.--        For iterative solvers, the `compute` step is used to eventually setup a preconditioner.-        Remember that, basically, the goal of the preconditioner is to speedup the convergence of an iterative method by solving a modified linear-        system where the coefficient matrix has more clustered eigenvalues.-        For real problems, an iterative solver should always be used with a preconditioner.-    -}-    analyzePattern,-    factorize,-    compute,-    -- * The Solve step-    {- |-    The `solve` function computes the solution of the linear systems with one or many right hand sides.--    @-    x <- solve b-    @--    Here, @b@ can be a vector or a matrix where the columns form the different right hand sides.-    The `solve` function can be called several times as well, for instance when all the right hand sides are not available at once.--    @-    x1 <- solve b1-    -- Get the second right hand side b2-    x2 <- solve b2-    --  ...-    @-    -}-    solve,-    --solveWithGuess,-    info,-    -- * Iterative Solvers-    tolerance,-    setTolerance,-    maxIterations,-    setMaxIterations,-    Data.Eigen.SparseLA.error,-    iterations,-    -- * SparseQR Solver-    matrixR,-    matrixQ,-    rank,-    setPivotThreshold,-    -- * SparseLU Solver-    setSymmetric,-    matrixL,-    matrixU,-    determinant,-    absDeterminant,-    signDeterminant,-    logAbsDeterminant,-) where--import Prelude as P-import Foreign.Ptr-import Foreign.ForeignPtr-import Foreign.Storable-import Foreign.C.String-import Foreign.Marshal.Alloc-import Control.Monad.IO.Class-import Control.Monad.Trans.Reader-import qualified Foreign.Concurrent as FC-#if __GLASGOW_HASKELL__ >= 710-#else-import Control.Applicative-#endif-import qualified Data.Eigen.Internal as I-import qualified Data.Eigen.SparseMatrix as SM--{- | Ordering methods for sparse matrices. They are typically used to reduce the number of elements during the sparse matrix-    decomposition (@LLT@, @LU@, @QR@). Precisely, in a preprocessing step, a permutation matrix @P@ is computed using those ordering methods-    and applied to the columns of the matrix. Using for instance the sparse Cholesky decomposition, it is expected that the nonzeros-    elements in @LLT(A*P)@ will be much smaller than that in @LLT(A)@.--}-data OrderingMethod-    -- | The column approximate minimum degree ordering The matrix should be in column-major and compressed format-    = COLAMDOrdering-    -- | The natural ordering (identity)-    | NaturalOrdering deriving (Show, Read)--data Preconditioner-    {- | A preconditioner based on the digonal entries--        It allows to approximately solve for A.x = b problems assuming A is a diagonal matrix.-        In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for:-        @-        A.diagonal().asDiagonal() . x = b-        @-        This preconditioner is suitable for both selfadjoint and general problems.-        The diagonal entries are pre-inverted and stored into a dense vector.--        A variant that has yet to be implemented would attempt to preserve the norm of each column.-    -}-    = DiagonalPreconditioner-    -- | A naive preconditioner which approximates any matrix as the identity matrix-    | IdentityPreconditioner deriving (Show, Read)---class I.Code s => Solver s where--- | For direct methods, the solution is computed at the machine precision.-class Solver s => DirectSolver s where--- | Sometimes, the solution need not be too accurate.--- In this case, the iterative methods are more suitable and the desired accuracy can be set before the solve step using `setTolerance`.-class Solver s => IterativeSolver s where--{- | A conjugate gradient solver for sparse self-adjoint problems.--    This class allows to solve for @A.x = b@ sparse linear problems using a conjugate gradient algorithm. The sparse matrix @A@ must be selfadjoint.--    The maximal number of iterations and tolerance value can be controlled via the `setMaxIterations` and `setTolerance` methods.-    The defaults are the size of the problem for the maximal number of iterations and @epsilon@ for the tolerance--}-data ConjugateGradient = ConjugateGradient Preconditioner deriving (Show, Read)-instance Solver ConjugateGradient-instance IterativeSolver ConjugateGradient-instance I.Code ConjugateGradient where-    code (ConjugateGradient DiagonalPreconditioner) = 0-    code (ConjugateGradient IdentityPreconditioner) = 1--{- | A bi conjugate gradient stabilized solver for sparse square problems.--    This class allows to solve for @A.x = b@ sparse linear problems using a bi conjugate gradient stabilized algorithm.-    The vectors @x@ and @b@ can be either dense or sparse.--    The maximal number of iterations and tolerance value can be controlled via the `setMaxIterations` and `setTolerance` methods.-    The defaults are the size of the problem for the maximal number of iterations and @epsilon@ for the tolerance--}-data BiCGSTAB = BiCGSTAB Preconditioner deriving (Show, Read)-instance Solver BiCGSTAB-instance IterativeSolver BiCGSTAB-instance I.Code BiCGSTAB where-    code (BiCGSTAB DiagonalPreconditioner) = 2-    code (BiCGSTAB IdentityPreconditioner) = 3--{- | Sparse supernodal LU factorization for general matrices.--    This class implements the supernodal LU factorization for general matrices. It uses the main techniques from the sequential-    <http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ SuperLU package>. It handles transparently real and complex arithmetics with-    single and double precision, depending on the scalar type of your input matrix. The code has been optimized to provide BLAS-3-    operations during supernode-panel updates. It benefits directly from the built-in high-performant Eigen BLAS routines.-    Moreover, when the size of a supernode is very small, the BLAS calls are avoided to enable a better optimization from the compiler.-    For best performance, you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors.--    An important parameter of this class is the ordering method. It is used to reorder the columns-    (and eventually the rows) of the matrix to reduce the number of new elements that are created during-    numerical factorization. The cheapest method available is COLAMD.-    See <http://eigen.tuxfamily.org/dox/group__OrderingMethods__Module.html OrderingMethods module> for the list of-    built-in and external ordering methods.--}-data SparseLU = SparseLU OrderingMethod deriving (Show, Read)-instance Solver SparseLU-instance DirectSolver SparseLU-instance I.Code SparseLU where-    code (SparseLU NaturalOrdering) = 4-    code (SparseLU COLAMDOrdering) = 5--{- | Sparse left-looking rank-revealing QR factorization.--    This class implements a left-looking rank-revealing QR decomposition of sparse matrices. When a column has a norm less than a given-    tolerance it is implicitly permuted to the end. The QR factorization thus obtained is given by @A*P = Q*R@ where @R@ is upper triangular or trapezoidal.--    @P@ is the column permutation which is the product of the fill-reducing and the rank-revealing permutations.--    @Q@ is the orthogonal matrix represented as products of Householder reflectors.--    @R@ is the sparse triangular or trapezoidal matrix. The later occurs when @A@ is rank-deficient.--}-data SparseQR = SparseQR OrderingMethod deriving (Show, Read)-instance Solver SparseQR-instance DirectSolver SparseQR-instance I.Code SparseQR where-    code (SparseQR NaturalOrdering) = 6-    code (SparseQR COLAMDOrdering) = 7---data ComputationInfo-    -- | Computation was successful.-    = Success-    -- | The provided data did not satisfy the prerequisites.-    | NumericalIssue-    -- | Iterative procedure did not converge.-    | NoConvergence-    -- | The inputs are invalid, or the algorithm has been improperly called. When assertions are enabled, such errors trigger an error.-    | InvalidInput-    deriving (Eq, Enum, Show, Read)--type SolverT s a b m = ReaderT (s, ForeignPtr (I.CSolver a b)) m--runSolverT :: (Solver s, MonadIO m, I.Elem a b) => s -> SolverT s a b m c -> m c-runSolverT i f = do-    fs <- liftIO $ alloca $ \ps -> do-        I.call $ I.sparse_la_newSolver i ps-        s <- peek ps-        FC.newForeignPtr s (I.call $ I.sparse_la_freeSolver i s)-    runReaderT f (i,fs)---- | Initializes the iterative solver for the sparsity pattern of the matrix @A@ for further solving @Ax=b@ problems.-analyzePattern :: (Solver s, MonadIO m, I.Elem a b) => SM.SparseMatrix a b -> SolverT s a b m ()-analyzePattern (SM.SparseMatrix fa) = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s ->-    withForeignPtr fa $ \a ->-        I.call $ I.sparse_la_analyzePattern i s a---- | Initializes the iterative solver with the numerical values of the matrix @A@ for further solving @Ax=b@ problems.-factorize :: (Solver s, MonadIO m, I.Elem a b) => SM.SparseMatrix a b -> SolverT s a b m ()-factorize (SM.SparseMatrix fa) = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s ->-    withForeignPtr fa $ \a ->-        I.call $ I.sparse_la_factorize i s a---- | Initializes the iterative solver with the matrix @A@ for further solving @Ax=b@ problems.------ The `compute` method is equivalent to calling both `analyzePattern` and `factorize`.-compute :: (Solver s, MonadIO m, I.Elem a b) => SM.SparseMatrix a b -> SolverT s a b m ()-compute (SM.SparseMatrix fa) = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s ->-    withForeignPtr fa $ \a ->-        I.call $ I.sparse_la_compute i s a---- | An expression of the solution @x@ of @Ax=b@ using the current decomposition of @A@.-solve :: (Solver s, MonadIO m, I.Elem a b) => SM.SparseMatrix a b -> SolverT s a b m (SM.SparseMatrix a b)-solve (SM.SparseMatrix fb) = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s ->-    withForeignPtr fb $ \b ->-    alloca $ \px -> do-        I.call $ I.sparse_la_solve i s b px-        x <- peek px-        SM.SparseMatrix <$> FC.newForeignPtr x (I.call $ I.sparse_free x)--{---- | The solution @x@ of @Ax=b@ using the current decomposition of @A@ and @x0@ as an initial solution.-solveWithGuess :: (MonadIO m, I.Elem a b) => SM.SparseMatrix a b -> SM.SparseMatrix a b -> SolverT s a b m (SM.SparseMatrix a b)-solveWithGuess (SM.SparseMatrix fb) (SM.SparseMatrix fx0) = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s ->-    withForeignPtr fb $ \b ->-    withForeignPtr fx0 $ \x0 ->-    alloca $ \px -> do-        I.call $ I.sparse_la_solveWithGuess i s b x0 px-        x <- peek px-        SM.SparseMatrix <$> FC.newForeignPtr x (I.call $ I.sparse_free x)--}---- |--- * `Success` if the iterations converged or computation was succesful--- * `NumericalIssue` if the factorization reports a numerical problem--- * `NoConvergence` if the iterations are not converged--- * `InvalidInput` if the input matrix is invalid-info :: (Solver s, MonadIO m, I.Elem a b) => SolverT s a b m ComputationInfo-info = _get_prop I.sparse_la_info >>= \x -> return (toEnum x)---- | The tolerance threshold used by the stopping criteria.-tolerance :: (IterativeSolver s, MonadIO m, I.Elem a b) => SolverT s a b m Double-tolerance = _get_prop I.sparse_la_tolerance---- | Sets the tolerance threshold used by the stopping criteria.------   This value is used as an upper bound to the relative residual error: @|Ax-b|/|b|@. The default value is the machine precision given by @epsilon@-setTolerance :: (IterativeSolver s, MonadIO m, I.Elem a b) => Double -> SolverT s a b m ()-setTolerance = _set_prop I.sparse_la_setTolerance---- | The max number of iterations. It is either the value setted by setMaxIterations or, by default, twice the number of columns of the matrix.-maxIterations :: (IterativeSolver s, MonadIO m, I.Elem a b) => SolverT s a b m Int-maxIterations = _get_prop I.sparse_la_maxIterations---- | Sets the max number of iterations. Default is twice the number of columns of the matrix.-setMaxIterations :: (IterativeSolver s, MonadIO m, I.Elem a b) => Int -> SolverT s a b m ()-setMaxIterations = _set_prop I.sparse_la_setMaxIterations---- | The tolerance error reached during the last solve. It is a close approximation of the true relative residual error @|Ax-b|/|b|@.-error :: (IterativeSolver s, MonadIO m, I.Elem a b) => SolverT s a b m Double-error = _get_prop I.sparse_la_error---- | The number of iterations performed during the last solve-iterations :: (IterativeSolver s, MonadIO m, I.Elem a b) => SolverT s a b m Int-iterations = _get_prop I.sparse_la_iterations---- | Returns the @b@ sparse upper triangular matrix @R@ of the QR factorization.-matrixR :: (MonadIO m, I.Elem a b) => SolverT SparseQR a b m (SM.SparseMatrix a b)-matrixR = _get_matrix I.sparse_la_matrixR---- | Returns the matrix @Q@ as products of sparse Householder reflectors.-matrixQ :: (MonadIO m, I.Elem a b) => SolverT SparseQR a b m (SM.SparseMatrix a b)-matrixQ = _get_matrix I.sparse_la_matrixQ---- | Sets the threshold that is used to determine linearly dependent columns during the factorization.------ In practice, if during the factorization the norm of the column that has to be eliminated is below--- this threshold, then the entire column is treated as zero, and it is moved at the end.-setPivotThreshold :: (MonadIO m, I.Elem a b) => Double -> SolverT SparseQR a b m ()-setPivotThreshold = _set_prop I.sparse_la_setPivotThreshold---- | Returns the number of non linearly dependent columns as determined by the pivoting threshold.-rank :: (MonadIO m, I.Elem a b) => SolverT SparseQR a b m Int-rank = _get_prop I.sparse_la_rank---- | Indicate that the pattern of the input matrix is symmetric-setSymmetric :: (MonadIO m, I.Elem a b) => Bool -> SolverT SparseLU a b m ()-setSymmetric = _set_prop I.sparse_la_setSymmetric . fromEnum---- | Returns the matrix @L@-matrixL :: (MonadIO m, I.Elem a b) => SolverT SparseLU a b m (SM.SparseMatrix a b)-matrixL = _get_matrix I.sparse_la_matrixL---- | Returns the matrix @U@-matrixU :: (MonadIO m, I.Elem a b) => SolverT SparseLU a b m (SM.SparseMatrix a b)-matrixU = _get_matrix I.sparse_la_matrixU---- | The determinant of the matrix.-determinant :: (MonadIO m, I.Elem a b) => SolverT SparseLU a b m a-determinant = _get_prop I.sparse_la_determinant---- | The natural log of the absolute value of the determinant of the matrix of which this is the QR decomposition------ This method is useful to work around the risk of overflow/underflow that's inherent to the determinant computation.-logAbsDeterminant :: (MonadIO m, I.Elem a b) => SolverT SparseLU a b m a-logAbsDeterminant = _get_prop I.sparse_la_logAbsDeterminant---- | The absolute value of the determinant of the matrix of which *this is the QR decomposition.------ A determinant can be very big or small, so for matrices of large enough dimension, there is a risk of overflow/underflow.--- One way to work around that is to use `logAbsDeterminant` instead.-absDeterminant :: (MonadIO m, I.Elem a b) => SolverT SparseLU a b m a-absDeterminant = _get_prop I.sparse_la_absDeterminant---- | A number representing the sign of the determinant-signDeterminant :: (MonadIO m, I.Elem a b) => SolverT SparseLU a b m a-signDeterminant = _get_prop I.sparse_la_signDeterminant--_get_prop :: (I.Cast c d, Solver s, MonadIO m, Storable c) => (s -> I.CSolverPtr a b -> Ptr c -> IO CString) -> SolverT s a b m d-_get_prop f = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s -> alloca $ \px -> do-        I.call $ f i s px-        I.cast <$> peek px--_get_matrix :: (Solver s, MonadIO m, I.Elem a b) => (s -> I.CSolverPtr a b -> Ptr (I.CSparseMatrixPtr a b) -> IO CString) -> SolverT s a b m (SM.SparseMatrix a b)-_get_matrix f = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s -> alloca $ \px -> do-        I.call $ f i s px-        x <- peek px-        SM.SparseMatrix <$> FC.newForeignPtr x (I.call $ I.sparse_free x)--_set_prop :: (I.Cast c d, Solver s, MonadIO m, Storable c) => (s -> I.CSolverPtr a b -> d -> IO CString) -> c -> SolverT s a b m ()-_set_prop f x = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s -> I.call $ f i s (I.cast x)

@@ -1,538 +0,0 @@-{-# LANGUAGE CPP #-}-{-# LANGUAGE FlexibleContexts #-}-{-# LANGUAGE GADTs #-}-{-# LANGUAGE ScopedTypeVariables #-}--module Data.Eigen.SparseMatrix (-    -- * SparseMatrix type-    SparseMatrix(..),-    SparseMatrixXf,-    SparseMatrixXd,-    SparseMatrixXcf,-    SparseMatrixXcd,-    -- * Matrix internal data-    values,-    innerIndices,-    outerStarts,-    innerNNZs,-    -- * Accessing matrix data-    cols,-    rows,-    coeff,-    (!),-    getRow,-    getCol,-    getRows,-    getCols,-    squareSubset,-    -- * Matrix conversions-    fromList,-    toList,-    fromVector,-    toVector,-    fromDenseList,-    toDenseList,-    fromMatrix,-    toMatrix,-    ones,-    ident,-    diagCol,-    diagRow,-    fromRows,-    fromCols,-    -- * Matrix properties-    norm,-    squaredNorm,-    blueNorm,-    block,-    nonZeros,-    innerSize,-    outerSize,-    getRowSums,-    getColSums,-    getSum,-    -- * Basic matrix algebra-    add,-    sub,-    mul,-    -- * Matrix transformations-    pruned,-    scale,-    transpose,-    adjoint,-    _map,-    _imap,-    -- * Matrix representation-    compress,-    uncompress,-    compressed,-    -- * Matrix serialization-    encode,-    decode,-    -- * Mutable matricies-    thaw,-    freeze,-    unsafeThaw,-    unsafeFreeze,-) where--import qualified Prelude as P-import Prelude hiding (map)-import qualified Data.List as L-import Data.Complex-import Data.Binary hiding (encode, decode)-import qualified Data.Binary as B-import Foreign.C.Types-import Foreign.C.String-import Foreign.Storable-import Foreign.Ptr-import Foreign.ForeignPtr-import Foreign.Marshal.Alloc-import Control.Monad-#if __GLASGOW_HASKELL__ >= 710-#else-import Control.Applicative-#endif-import qualified Data.Eigen.Matrix as M-import qualified Data.Eigen.Matrix.Mutable as MM-import qualified Data.Eigen.SparseMatrix.Mutable as SMM-import qualified Foreign.Concurrent as FC-import qualified Data.Eigen.Internal as I-import qualified Data.Vector.Storable as VS-import qualified Data.Vector.Storable.Mutable as VSM-import qualified Data.ByteString.Lazy as BSL--{-| A versatible sparse matrix representation.--SparseMatrix is the main sparse matrix representation of Eigen's sparse module.-It offers high performance and low memory usage.-It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme.--It consists of four compact arrays:--* `values`: stores the coefficient values of the non-zeros.-* `innerIndices`: stores the row (resp. column) indices of the non-zeros.-* `outerStarts`: stores for each column (resp. row) the index of the first non-zero in the previous two arrays.-* `innerNNZs`: stores the number of non-zeros of each column (resp. row). The word inner refers to an inner vector that is a column for a column-major matrix, or a row for a row-major matrix. The word outer refers to the other direction.--This storage scheme is better explained on an example. The following matrix--@-0   3   0   0   0-22  0   0   0   17-7   5   0   1   0-0   0   0   0   0-0   0   14  0   8-@--and one of its possible sparse, __column major__ representation:--@-values:         22  7   _   3   5   14  _   _   1   _   17  8-innerIndices:   1   2   _   0   2   4   _   _   2   _   1   4-outerStarts:    0   3   5   8   10  12-innerNNZs:      2   2   1   1   2-@--Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices.-The "\_" indicates available free space to quickly insert new elements. Assuming no reallocation is needed,-the insertion of a random element is therefore in @O(nnz_j)@ where @nnz_j@ is the number of nonzeros of the-respective inner vector. On the other hand, inserting elements with increasing inner indices in a given inner-vector is much more efficient since this only requires to increase the respective `innerNNZs` entry that is a @O(1)@ operation.--The case where no empty space is available is a special case, and is refered as the compressed mode.-It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS).-Any `SparseMatrix` can be turned to this form by calling the `compress` function.-In this case, one can remark that the `innerNNZs` array is redundant with `outerStarts` because we the equality:-@InnerNNZs[j] = OuterStarts[j+1]-OuterStarts[j]@. Therefore, in practice a call to `compress` frees this buffer.--The results of Eigen's operations always produces compressed sparse matrices.-On the other hand, the insertion of a new element into a `SparseMatrix` converts this later to the uncompressed mode.--For more infomration please see Eigen <http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html documentation page>.--}--data SparseMatrix a b where-    SparseMatrix :: I.Elem a b => !(ForeignPtr (I.CSparseMatrix a b)) -> SparseMatrix a b---- | Alias for single precision sparse matrix-type SparseMatrixXf = SparseMatrix Float CFloat--- | Alias for double precision sparse matrix-type SparseMatrixXd = SparseMatrix Double CDouble--- | Alias for single previsiom sparse matrix of complex numbers-type SparseMatrixXcf = SparseMatrix (Complex Float) (I.CComplex CFloat)--- | Alias for double prevision sparse matrix of complex numbers-type SparseMatrixXcd = SparseMatrix (Complex Double) (I.CComplex CDouble)---- | Pretty prints the sparse matrix-instance (I.Elem a b, Show a) => Show (SparseMatrix a b) where-    show m = concat [-        "SparseMatrix ", show (rows m), "x", show (cols m),-        "\n", L.intercalate "\n" $ P.map (L.intercalate "\t" . P.map show) $ toDenseList m, "\n"]---- | Basic sparse matrix math exposed through Num instance: @(*)@, @(+)@, @(-)@, `fromInteger`, `signum`, `abs`, `negate`-instance I.Elem a b => Num (SparseMatrix a b) where-    (*) = mul-    (+) = add-    (-) = sub-    fromInteger x = fromList 1 1 [(0,0,fromInteger x)]-    signum = _map signum-    abs = _map abs-    negate = _map negate--instance I.Elem a b => Binary (SparseMatrix a b) where-    put m = do-        put $ I.magicCode (undefined :: b)-        put $ rows m-        put $ cols m-        put $ toVector m--    get = do-        get >>= (`when` fail "wrong matrix type") . (/= I.magicCode (undefined :: b))-        fromVector <$> get <*> get <*> get---- | Encode the sparse matrix as a lazy byte string-encode :: I.Elem a b => SparseMatrix a b -> BSL.ByteString-encode = B.encode----- | Decode sparse matrix from the lazy byte string-decode :: I.Elem a b => BSL.ByteString -> SparseMatrix a b-decode = B.decode---- | Stores the coefficient values of the non-zeros.-values :: I.Elem a b => SparseMatrix a b -> VS.Vector b-values = _getvec I.sparse_values---- | Stores the row (resp. column) indices of the non-zeros.-innerIndices :: I.Elem a b => SparseMatrix a b -> VS.Vector CInt-innerIndices = _getvec I.sparse_innerIndices---- | Stores for each column (resp. row) the index of the first non-zero in the previous two arrays.-outerStarts :: I.Elem a b => SparseMatrix a b -> VS.Vector CInt-outerStarts = _getvec I.sparse_outerStarts---- | Stores the number of non-zeros of each column (resp. row).--- The word inner refers to an inner vector that is a column for a column-major matrix, or a row for a row-major matrix.--- The word outer refers to the other direction-innerNNZs :: I.Elem a b => SparseMatrix a b -> Maybe (VS.Vector CInt)-innerNNZs m-    | compressed m = Nothing-    | otherwise = Just $ _getvec I.sparse_innerNNZs m---- | Number of rows for the sparse matrix-rows :: I.Elem a b => SparseMatrix a b -> Int-rows = _unop I.sparse_rows (return . I.cast)---- | Number of columns for the sparse matrix-cols :: I.Elem a b => SparseMatrix a b -> Int-cols = _unop I.sparse_cols (return . I.cast)---- | Matrix coefficient at given row and col-coeff :: I.Elem a b => Int -> Int -> SparseMatrix a b -> a-coeff row col (SparseMatrix fp) = I.performIO $ withForeignPtr fp $ \p -> alloca $ \pq -> do-    I.call $ I.sparse_coeff p (I.cast row) (I.cast col) pq-    I.cast <$> peek pq---- | Matrix coefficient at given row and col-(!) :: I.Elem a b => SparseMatrix a b -> (Int, Int) -> a-(!) m (row, col) = coeff row col m--{-| For vectors, the l2 norm, and for matrices the Frobenius norm.-    In both cases, it consists in the square root of the sum of the square of all the matrix entries.-    For vectors, this is also equals to the square root of the dot product of this with itself.--}-norm :: I.Elem a b => SparseMatrix a b -> a-norm = _unop I.sparse_norm (return . I.cast)---- | For vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.-squaredNorm :: I.Elem a b => SparseMatrix a b -> a-squaredNorm = _unop I.sparse_squaredNorm (return . I.cast)---- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.-blueNorm :: I.Elem a b => SparseMatrix a b -> a-blueNorm = _unop I.sparse_blueNorm (return . I.cast)---- | Extract rectangular block from sparse matrix defined by startRow startCol blockRows blockCols-block :: I.Elem a b => Int -> Int -> Int -> Int -> SparseMatrix a b -> SparseMatrix a b-block row col rows cols = _unop (\p pq -> I.sparse_block p (I.cast row) (I.cast col) (I.cast rows) (I.cast cols) pq) _mk---- | Number of non-zeros elements in the sparse matrix-nonZeros :: I.Elem a b => SparseMatrix a b -> Int-nonZeros = _unop I.sparse_nonZeros (return . I.cast)---- | The matrix in the compressed format-compress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-compress = _unop I.sparse_makeCompressed _mk---- | The matrix in the uncompressed mode-uncompress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-uncompress = _unop I.sparse_uncompress _mk---- | Is this in compressed form?-compressed :: I.Elem a b => SparseMatrix a b -> Bool-compressed = _unop I.sparse_isCompressed (return . (/=0))---- | Minor dimension with respect to the storage order-innerSize :: I.Elem a b => SparseMatrix a b -> Int-innerSize = _unop I.sparse_innerSize (return . I.cast)---- | Major dimension with respect to the storage order-outerSize :: I.Elem a b => SparseMatrix a b -> Int-outerSize = _unop I.sparse_outerSize (return . I.cast)---- | Suppresses all nonzeros which are much smaller than reference under the tolerence @epsilon@-pruned :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b-pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (I.cast r) >> I.sparse_prunedRef p pr pq) _mk---- | Multiply matrix on a given scalar-scale :: I.Elem a b => a -> SparseMatrix a b -> SparseMatrix a b-scale x = _unop (\p pq -> alloca $ \px -> poke px (I.cast x) >> I.sparse_scale p px pq) _mk---- | Transpose of the sparse matrix-transpose :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-transpose = _unop I.sparse_transpose _mk---- | Adjoint of the sparse matrix-adjoint :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b-adjoint = _unop I.sparse_adjoint _mk---- | Adding two sparse matrices by adding the corresponding entries together. You can use @(+)@ function as well.-add :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b-add = _binop I.sparse_add _mk---- | Subtracting two sparse matrices by subtracting the corresponding entries together. You can use @(-)@ function as well.-sub :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b-sub = _binop I.sparse_sub _mk---- | Matrix multiplication. You can use @(*)@ function as well.-mul :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b -> SparseMatrix a b-mul = _binop I.sparse_mul _mk---- | Construct sparse matrix of given size from the list of triplets (row, col, val)-fromList :: I.Elem a b => Int -> Int -> [(Int, Int, a)] -> SparseMatrix a b-fromList rows cols = fromVector rows cols . VS.fromList . P.map I.cast---- | Construct sparse matrix of given size from the storable vector of triplets (row, col, val)-fromVector :: I.Elem a b => Int -> Int -> VS.Vector (I.CTriplet b) -> SparseMatrix a b-fromVector rows cols tris = I.performIO $ VS.unsafeWith tris $ \p -> alloca $ \pq -> do-    I.call $ I.sparse_fromList (I.cast rows) (I.cast cols) p (I.cast $ VS.length tris) pq-    peek pq >>= _mk---- | Convert sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included-toList :: I.Elem a b => SparseMatrix a b -> [(Int, Int, a)]-toList = P.map I.cast . VS.toList . toVector---- | Convert sparse matrix to the storable vector of triplets (row, col, val). Compressed elements will not be included-toVector :: I.Elem a b => SparseMatrix a b -> VS.Vector (I.CTriplet b)-toVector m@(SparseMatrix fp) = I.performIO $ do-    let size = nonZeros m-    tris <- VSM.new size-    withForeignPtr fp $ \p ->-        VSM.unsafeWith tris $ \q ->-            I.call $ I.sparse_toList p q (I.cast size)-    VS.unsafeFreeze tris---- | Construct sparse matrix of two-dimensional list of values. Matrix dimensions will be detected automatically. Zero values will be compressed.-fromDenseList :: (I.Elem a b, Eq a) => [[a]] -> SparseMatrix a b-fromDenseList list = fromList rows cols $ do-    (row, vals) <- zip [0..] list-    (col, val) <- zip [0..] vals-    guard $ val /= 0-    return (row, col, val)-    where-        rows = length list-        cols = L.foldl' max 0 $ P.map length list---- | Convert sparse matrix to (rows X cols) dense list of values-toDenseList :: I.Elem a b => SparseMatrix a b -> [[a]]-toDenseList m = [[coeff row col m | col <- [0 .. cols m - 1]] | row <- [0 .. rows m - 1]]---- | Construct sparse matrix from dense matrix. Zero elements will be compressed-fromMatrix :: I.Elem a b => M.Matrix a b -> SparseMatrix a b-fromMatrix m1 = I.performIO $ alloca $ \pm0 ->-    M.unsafeWith m1 $ \vals rows cols -> do-        I.call $ I.sparse_fromMatrix vals rows cols pm0-        peek pm0 >>= _mk---- | Construct dense matrix from sparse matrix-toMatrix :: I.Elem a b => SparseMatrix a b -> M.Matrix a b-toMatrix m1@(SparseMatrix fp) = I.performIO $ do-    m0 <- MM.new (rows m1) (cols m1)-    MM.unsafeWith m0 $ \vals rows cols ->-        withForeignPtr fp $ \pm1 ->-            I.call $ I.sparse_toMatrix pm1 vals rows cols-    M.unsafeFreeze m0---- | Yield an immutable copy of the mutable matrix-freeze :: I.Elem a b => SMM.IOSparseMatrix a b -> IO (SparseMatrix a b)-freeze (SMM.IOSparseMatrix fp) = SparseMatrix <$> _clone fp---- | Yield a mutable copy of the immutable matrix-thaw :: I.Elem a b => SparseMatrix a b -> IO (SMM.IOSparseMatrix a b)-thaw (SparseMatrix fp) = SMM.IOSparseMatrix <$> _clone fp---- | Unsafe convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.-unsafeFreeze :: I.Elem a b => SMM.IOSparseMatrix a b -> IO (SparseMatrix a b)-unsafeFreeze (SMM.IOSparseMatrix fp) = return $! SparseMatrix fp---- | Unsafely convert an immutable matrix to a mutable one without copying. The immutable matrix may not be used after this operation.-unsafeThaw :: I.Elem a b => SparseMatrix a b -> IO (SMM.IOSparseMatrix a b)-unsafeThaw (SparseMatrix fp) = return $! SMM.IOSparseMatrix fp--_unop :: Storable c => (I.CSparseMatrixPtr a b -> Ptr c -> IO CString) -> (c -> IO d) -> SparseMatrix a b -> d-_unop f g (SparseMatrix fp) = I.performIO $-    withForeignPtr fp $ \p ->-        alloca $ \pq -> do-            I.call (f p pq)-            peek pq >>= g--_binop :: Storable c => (I.CSparseMatrixPtr a b -> I.CSparseMatrixPtr a b -> Ptr c -> IO CString) -> (c -> IO d) -> SparseMatrix a b -> SparseMatrix a b -> d-_binop f g (SparseMatrix fp1) (SparseMatrix fp2) = I.performIO $-    withForeignPtr fp1 $ \p1 ->-        withForeignPtr fp2 $ \p2 ->-            alloca $ \pq -> do-                I.call (f p1 p2 pq)-                peek pq >>= g--_getvec :: (I.Elem a b, Storable c) => (Ptr (I.CSparseMatrix a b) -> Ptr CInt -> Ptr (Ptr c) -> IO CString) -> SparseMatrix a b -> VS.Vector c-_getvec f (SparseMatrix fm) = I.performIO $-    withForeignPtr fm $ \m ->-    alloca $ \ps ->-    alloca $ \pq -> do-        I.call $ f m ps pq-        s <- fromIntegral <$> peek ps-        q <- peek pq-        fr <- FC.newForeignPtr q $ touchForeignPtr fm-        return $! VS.unsafeFromForeignPtr0 fr s--_clone :: I.Elem a b => ForeignPtr (I.CSparseMatrix a b) -> IO (ForeignPtr (I.CSparseMatrix a b))-_clone fp = withForeignPtr fp $ \p -> alloca $ \pq -> do-    I.call $ I.sparse_clone p pq-    q <- peek pq-    FC.newForeignPtr q $ I.call $ I.sparse_free q---- | Map over values of a sparse matrix.-_map :: I.Elem a b => (a -> a) -> SparseMatrix a b -> SparseMatrix a b-_map f m = fromVector (rows m) (cols m) . VS.map g . toVector $ m where-    g (I.CTriplet r c v) = I.CTriplet r c $ I.cast $ f $ I.cast v---- | Map over values of a sparse matrix with indices.-_imap-    :: I.Elem a b-    => (Int -> Int -> a -> a) -> SparseMatrix a b -> SparseMatrix a b-_imap f m = fromVector (rows m) (cols m) . VS.map g . toVector $ m-  where-    g (I.CTriplet r c v) =-        I.CTriplet r c $ I.cast $ f (fromIntegral r) (fromIntegral c) $ I.cast v--_mk :: I.Elem a b => Ptr (I.CSparseMatrix a b) -> IO (SparseMatrix a b)-_mk p = SparseMatrix <$> FC.newForeignPtr p (I.call $ I.sparse_free p)---- | Get a row of a sparse matrix.-getRow :: I.Elem a b => Int -> SparseMatrix a b -> SparseMatrix a b-getRow row mat = block row 0 1 m mat-  where-    m = cols mat---- | Get a column of a sparse matrix.-getCol :: I.Elem a b => Int -> SparseMatrix a b -> SparseMatrix a b-getCol col mat = block 0 col n 1 mat-  where-    n = rows mat---- | Get all rows of a sparse matrix.-getRows :: I.Elem a b => SparseMatrix a b -> [SparseMatrix a b]-getRows mat = fmap (flip getRow mat) [0 .. n - 1]-  where-    n = rows mat---- | Get all columns of a sparse matrix.-getCols :: I.Elem a b => SparseMatrix a b -> [SparseMatrix a b]-getCols mat = fmap (flip getCol mat) [0 .. m - 1]-  where-    m = cols mat---- | Get all row sums.-getRowSums :: SparseMatrixXd -> SparseMatrixXd-getRowSums = fromDenseList-           . fmap ((:[]) . sum . fmap (\(_, _, x) -> x) . toList)-           . getRows---- | Get all column sums.-getColSums :: SparseMatrixXd -> SparseMatrixXd-getColSums = fromDenseList-           . (:[])-           . fmap (sum . fmap (\(_, _, x) -> x) . toList)-           . getCols---- | Get sum of matrix.-getSum :: SparseMatrixXd -> Double-getSum = sum . fmap (\(_, _, x) -> x) . toList---- | Get the ones vector.-ones :: Int -> SparseMatrixXd-ones n = transpose . fromDenseList . (:[]) . replicate n $ 1---- | Get the identity matrix.-ident :: Int -> SparseMatrixXd-ident n = fromList n n . zip3 [0..] [0..] . replicate n $ 1---- | Transform a column into a diagonal matrix.-diagCol :: I.Elem a b => Int -> SparseMatrix a b -> SparseMatrix a b-diagCol col mat =-    fromList n m . fmap (\(i, _, v) -> (i, i, v)) . toList . getCol col $ mat-  where-    n = rows mat-    m = rows mat---- | Transform a row into a diagonal matrix.-diagRow :: I.Elem a b => Int -> SparseMatrix a b -> SparseMatrix a b-diagRow row mat =-    fromList n m . fmap (\(_, j, v) -> (j, j, v)) . toList . getRow row $ mat-  where-    n = cols mat-    m = cols mat---- | Update the row indices of a matrix.-updateRowIdx :: I.Elem a b => Int -> SparseMatrix a b -> SparseMatrix a b-updateRowIdx row mat =-    fromList n m . fmap (\(_, j, v) -> (row, j, v)) . toList $ mat-  where-    n = row + 1-    m = cols mat---- | Update the column indices of a matrix.-updateColIdx :: I.Elem a b => Int -> SparseMatrix a b -> SparseMatrix a b-updateColIdx col mat =-    fromList n m . fmap (\(i, _, v) -> (i, col, v)) . toList $ mat-  where-    n = rows mat-    m = col + 1---- | Get a matrix from a list of rows.-fromRows :: I.Elem a b => [SparseMatrix a b] -> SparseMatrix a b-fromRows []   = fromList 0 0 []-fromRows rows =-    fromList n m . concatMap toList . zipWith updateRowIdx [0..] $ rows-  where-    n = length rows-    m = maximum . fmap cols $ rows---- | Get a matrix from a list of cols.-fromCols :: I.Elem a b => [SparseMatrix a b] -> SparseMatrix a b-fromCols []   = fromList 0 0 []-fromCols cols =-    fromList n m . concatMap toList . zipWith updateColIdx [0..] $ cols-  where-    n = maximum . fmap rows $ cols-    m = length cols---- | Get a subset of a square matrix.-squareSubset :: I.Elem a b => [Int] -> SparseMatrix a b -> SparseMatrix a b-squareSubset [] _ = fromList 0 0 []-squareSubset idxs mat = fromCols-                      . (\m -> fmap (flip getCol m) idxs)-                      . fromRows-                      . fmap (flip getRow mat)-                      $ idxs

@@ -1,144 +0,0 @@-{-# LANGUAGE GADTs, RecordWildCards, ScopedTypeVariables #-}-module Data.Eigen.SparseMatrix.Mutable (-    -- * Mutable SparseMatrix-    IOSparseMatrix(..),-    IOSparseMatrixXf,-    IOSparseMatrixXd,-    IOSparseMatrixXcf,-    IOSparseMatrixXcd,-    new,-    reserve,-    -- * Matrix properties-    rows,-    cols,-    innerSize,-    outerSize,-    nonZeros,-    -- * Matrix compression-    compressed,-    compress,-    uncompress,-    -- * Accessing matrix data-    read,-    write,-    setZero,-    setIdentity,-    -- * Changing matrix shape-    resize,-    conservativeResize-) where--import Prelude hiding (read)-import Data.Complex-import Foreign.C.String-import Foreign.C.Types-import Foreign.ForeignPtr-import Foreign.Marshal.Alloc-import Foreign.Ptr-import Foreign.Storable-import qualified Foreign.Concurrent as FC-import qualified Data.Eigen.Internal as I---- | Mutable version of sparse matrix. See `Data.Eigen.SparseMatrix.SparseMatrix` for details about matrix layout.-data IOSparseMatrix a b where-    IOSparseMatrix :: I.Elem a b => !(ForeignPtr (I.CSparseMatrix a b)) -> IOSparseMatrix a b---- | Alias for single precision mutable matrix-type IOSparseMatrixXf = IOSparseMatrix Float CFloat--- | Alias for double precision mutable matrix-type IOSparseMatrixXd = IOSparseMatrix Double CDouble--- | Alias for single previsiom mutable matrix of complex numbers-type IOSparseMatrixXcf = IOSparseMatrix (Complex Float) (I.CComplex CFloat)--- | Alias for double prevision mutable matrix of complex numbers-type IOSparseMatrixXcd = IOSparseMatrix (Complex Double) (I.CComplex CDouble)----- | Creates new matrix with the given size @rows x cols@-new :: I.Elem a b => Int -> Int -> IO (IOSparseMatrix a b)-new rows cols = alloca $ \pm -> do-    I.call $ I.sparse_new (I.cast rows) (I.cast cols) pm-    m <- peek pm-    fm <- FC.newForeignPtr m $ I.call $ I.sparse_free m-    return $! IOSparseMatrix fm---- | Returns the number of rows of the matrix-rows :: I.Elem a b => IOSparseMatrix a b -> IO Int-rows = _prop I.sparse_rows (return . I.cast)---- | Returns the number of columns of the matrix-cols :: I.Elem a b => IOSparseMatrix a b  -> IO Int-cols = _prop I.sparse_cols (return . I.cast)---- | Returns the number of rows (resp. columns) of the matrix if the storage order column major (resp. row major)-innerSize :: I.Elem a b => IOSparseMatrix a b  -> IO Int-innerSize = _prop I.sparse_innerSize (return . I.cast)---- | Returns the number of columns (resp. rows) of the matrix if the storage order column major (resp. row major)-outerSize :: I.Elem a b => IOSparseMatrix a b  -> IO Int-outerSize = _prop I.sparse_outerSize (return . I.cast)---- | Returns whether this matrix is in compressed form.-compressed :: I.Elem a b => IOSparseMatrix a b -> IO Bool-compressed = _prop I.sparse_isCompressed (return . (==1))---- | Turns the matrix into the compressed format.-compress :: I.Elem a b => IOSparseMatrix a b -> IO ()-compress = _inplace I.sparse_compressInplace---- | Turns the matrix into the uncompressed mode.-uncompress :: I.Elem a b => IOSparseMatrix a b -> IO ()-uncompress = _inplace I.sparse_uncompressInplace---- | Reads the value of the matrix at position @i@, @j@.--- This function returns @Scalar(0)@ if the element is an explicit zero.-read :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> IO a-read (IOSparseMatrix fm) row col = withForeignPtr fm $ \m -> alloca $ \px -> do-    I.call $ I.sparse_coeff m (I.cast row) (I.cast col) px-    I.cast <$> peek px--{- | Writes the value of the matrix at position @i@, @j@.-    This function turns the matrix into a non compressed form if that was not the case.--    This is a @O(log(nnz_j))@ operation (binary search) plus the cost of element insertion if the element does not already exist.-        -    Cost of element insertion is sorted insertion in O(1) if the elements of each inner vector are inserted in increasing inner index order, and in @O(nnz_j)@ for a random insertion.--}-write :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> a -> IO ()-write (IOSparseMatrix fm) row col x = withForeignPtr fm $ \m -> alloca $ \px -> do-    I.call $ I.sparse_coeffRef m (I.cast row) (I.cast col) px-    peek px >>= (`poke` I.cast x)---- | Sets the matrix to the identity matrix-setIdentity :: I.Elem a b => IOSparseMatrix a b -> IO ()-setIdentity = _inplace I.sparse_setIdentity---- | Removes all non zeros but keep allocated memory-setZero :: I.Elem a b => IOSparseMatrix a b -> IO ()-setZero = _inplace I.sparse_setZero---- | The number of non zero coefficients-nonZeros :: I.Elem a b => IOSparseMatrix a b -> IO Int-nonZeros = _prop I.sparse_nonZeros (return . I.cast)---- | Preallocates space for non zeros. The matrix must be in compressed mode.-reserve :: I.Elem a b => IOSparseMatrix a b -> Int -> IO ()-reserve m s = _inplace (\p -> I.sparse_reserve p (I.cast s)) m---- | Resizes the matrix to a rows x cols matrix and initializes it to zero.-resize :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> IO ()-resize m rows cols = _inplace (\p -> I.sparse_resize p (I.cast rows) (I.cast cols)) m---- | Resizes the matrix to a rows x cols matrix leaving old values untouched.-conservativeResize :: I.Elem a b => IOSparseMatrix a b -> Int -> Int -> IO ()-conservativeResize m rows cols = _inplace (\p -> I.sparse_conservativeResize p (I.cast rows) (I.cast cols)) m--_inplace :: I.Elem a b => (Ptr (I.CSparseMatrix a b) -> IO CString) -> IOSparseMatrix a b -> IO ()-_inplace f (IOSparseMatrix fm) = withForeignPtr fm $ \m -> I.call $ f m--_prop :: Storable c => (I.CSparseMatrixPtr a b -> Ptr c -> IO CString) -> (c -> IO d) -> IOSparseMatrix a b -> IO d-_prop f g (IOSparseMatrix fp) =-    withForeignPtr fp $ \p ->-        alloca $ \pq -> do-            I.call (f p pq)-            peek pq >>= g-

@@ -190,8 +190,8 @@  // template <class T, class M, class S> // RET sparse_la_solveWithGuess(void* p, void* b, void* x0, void** x) {-//  *x = new M(((S*)p)->solveWithGuess(Matrix<T,Dynamic,Dynamic>(*(M*)b), Matrix<T,Dynamic,Dynamic>(*(M*)x0)));-//  return 0;+//     *x = new M(((S*)p)->solveWithGuess(Matrix<T,Dynamic,Dynamic>(*(M*)b), Matrix<T,Dynamic,Dynamic>(*(M*)x0)));+//     return 0; // } // API_ALL(sparse_la_solveWithGuess, (int code, int s, void* p, void* b, void* x0, void** x), (p,b,x0,x)); 

@@ -1,6 +1,12 @@ #include "eigen-sparse.h" #include <Eigen/LU> +#if __cplusplus > 199711L+  #define SMART_PTR std::unique_ptr+#else+  #define SMART_PTR std::auto_ptr+#endif+ template <class T> RET sparse_new(int rows, int cols, void** pr) {     typedef SparseMatrix<T> M;@@ -21,7 +27,7 @@ RET sparse_fromList(int rows, int cols, void* data, int size, void** pr) {     typedef SparseMatrix<T> M;     typedef Triplet<T> E;-    std::unique_ptr<M> a(new M(rows, cols));+    SMART_PTR<M> a(new M(rows, cols));     a->setFromTriplets((E*)data, (E*)data + size);     *(M**)pr = a.release();     return 0;

@@ -1,10 +1,14 @@-cabal-version:  2.2+cabal-version:  >= 1.10 name:           eigen-version:        3.3.4.1+version:        3.3.4.2 homepage:       https://github.com/chessai/eigen synopsis:       Eigen C++ library (linear algebra: matrices, sparse matrices, vectors, numerical solvers).-description:    This module provides Haskell binding for <http://eigen.tuxfamily.org/ Eigen C++ library>.+description:    __NOTE__: This library does not follow PVP. Instead, it follows Eigen's versioning. Version+                X.Y.Z.W means that the Eigen version is X.Y.Z, and the haskell release is W, where each W+                release is for bug fixes/API improvements.                 .+                This module provides Haskell binding for <http://eigen.tuxfamily.org/ Eigen C++ library>.+                .                 Eigen is versatile.                 .                 * It supports all matrix sizes, from small fixed-size matrices to arbitrarily large dense matrices, and even sparse matrices.@@ -48,14 +52,14 @@                 .                 Documentation at the most extent replicates original <http://eigen.tuxfamily.org/dox/ Eigen documentation>. category:       Data, Math, Algebra, Statistics, Algorithms, Numeric-license:        BSD-3-Clause+license:        BSD3 license-file:   LICENSE copyright:      (c) 2013-2015, Oleg Sidorkin,                 (c) 2018, chessai author:         Oleg Sidorkin <oleg.sidorkin@gmail.com> maintainer:     chessai <chessai1996@gmail.com> build-type:     Simple-tested-with:    GHC == 7.8.3, GHC == 8.0.2, GHC == 8.2.2, GHC == 8.4.3+tested-with:    GHC == 8.2.2, GHC == 8.4.3, GHC == 8.6.1 extra-source-files: cbits/eigen-runtime.h     cbits/eigen-dense.h     cbits/eigen-la.h@@ -689,54 +693,64 @@   library-    exposed-modules:    Data.Eigen.LA-                        Data.Eigen.SparseLA-                        Data.Eigen.Matrix-                        Data.Eigen.Matrix.Mutable-                        Data.Eigen.SparseMatrix-                        Data.Eigen.SparseMatrix.Mutable-                        Data.Eigen.Parallel-                        Data.Eigen.Internal--    ghc-options:        -Wall -fno-warn-name-shadowing-    build-depends:      base >= 4.9 && < 5,-                        vector >= 0.5 && < 0.13,-                        primitive >= 0.1 && < 0.7,-                        binary,-                        bytestring,-                        transformers >= 0.3-    default-language:   Haskell2010+    exposed-modules:+      Eigen.Matrix+      Eigen.Matrix.Mutable+      Eigen.SparseMatrix+      Eigen.SparseMatrix.Mutable+      Eigen.Parallel+      Eigen.Internal+      Eigen.Solver.LA+      Eigen.Solver.SparseLA -    include-dirs:       eigen3, cbits-    c-sources:          cbits/eigen-runtime.cpp-                        cbits/eigen-dense.cpp-                        cbits/eigen-sparse.cpp-                        cbits/eigen-la.cpp-                        cbits/eigen-sparse-la.cpp-    extra-libraries:    stdc++-    cxx-options:        --std=c++14+    hs-source-dirs:+      src+    ghc-options:+      -Wall+    build-depends:+        base >= 4.10 && < 5+      , binary >= 0.8.0.0 && < 0.8.6.0+      , bytestring >= 0.10.4.0 && < 0.11.0.0+      , constraints >= 0.10.0 && < 0.11.0+      , ghc-prim+      , primitive >= 0.6.4.0 && < 0.7+      , transformers >= 0.3 && < 0.6+      , vector >= 0.5 && < 0.13+    default-language:+      Haskell2010+    include-dirs:+        eigen3+      , cbits+    c-sources:+      cbits/eigen-runtime.cpp+      cbits/eigen-dense.cpp+      cbits/eigen-sparse.cpp+      cbits/eigen-la.cpp+      cbits/eigen-sparse-la.cpp+    extra-libraries:+      stdc++  Test-Suite test-solve-    type:               exitcode-stdio-1.0-    main-is:            test/solve.hs-    default-language:   Haskell2010-    build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen+  type:               exitcode-stdio-1.0+  main-is:            test/solve.hs+  default-language:   Haskell2010+  build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen, ghc-prim  Test-Suite test-solve-sparse-    type:               exitcode-stdio-1.0-    main-is:            test/solve-sparse.hs-    default-language:   Haskell2010-    build-depends:      base, primitive, vector, bytestring, transformers, mtl, binary, eigen+  type:               exitcode-stdio-1.0+  main-is:            test/solve-sparse.hs+  default-language:   Haskell2010+  build-depends:      base, primitive, vector, bytestring, transformers, mtl, binary, eigen  Test-Suite test-rank-    type:               exitcode-stdio-1.0-    main-is:            test/rank.hs-    default-language:   Haskell2010-    build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen+  type:               exitcode-stdio-1.0+  main-is:            test/rank.hs+  default-language:   Haskell2010+  build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen, ghc-prim  Test-Suite test-regression-    type:               exitcode-stdio-1.0-    main-is:            test/regression.hs-    default-language:   Haskell2010-    build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen+  type:               exitcode-stdio-1.0+  main-is:            test/regression.hs+  default-language:   Haskell2010+  build-depends:      base, primitive, vector, bytestring, transformers, binary, eigen, ghc-prim 

@@ -0,0 +1,369 @@+--------------------------------------------------------------------------------++{-# OPTIONS_GHC -fno-warn-orphans #-}++--------------------------------------------------------------------------------++{-# LANGUAGE AllowAmbiguousTypes       #-}+{-# LANGUAGE BangPatterns              #-}+{-# LANGUAGE CPP                       #-} +{-# LANGUAGE EmptyDataDecls            #-}+{-# LANGUAGE FlexibleContexts          #-}+{-# LANGUAGE FlexibleInstances         #-}+{-# LANGUAGE ForeignFunctionInterface  #-}+{-# LANGUAGE FunctionalDependencies    #-}+{-# LANGUAGE GADTs                     #-}+{-# LANGUAGE KindSignatures            #-}+{-# LANGUAGE LambdaCase                #-}+{-# LANGUAGE MagicHash                 #-}+{-# LANGUAGE MultiParamTypeClasses     #-}+{-# LANGUAGE Rank2Types                #-}+{-# LANGUAGE ScopedTypeVariables       #-}+{-# LANGUAGE StandaloneDeriving        #-}+{-# LANGUAGE TypeApplications          #-}+{-# LANGUAGE TypeFamilyDependencies    #-}+{-# LANGUAGE TypeInType                #-}+{-# LANGUAGE UnboxedTuples             #-}+{-# LANGUAGE UndecidableInstances      #-}++--------------------------------------------------------------------------------++-- | Internal module to Eigen.+--   Here we define all foreign function calls,+--   and some typeclasses integral to the public and private interfaces+--   of the library.+module Eigen.Internal where --   FIXME: Explicit export list++--------------------------------------------------------------------------------++import           Control.Monad            (when)+import           Data.Binary              (Binary(put,get))+import           Data.Binary.Get          (getByteString, getWord32be)+import           Data.Binary.Put          (putByteString, putWord32be)+import           Data.Bits                (xor)+import           Data.Complex             (Complex((:+)))+import           Data.Kind                (Type)+import           Data.Proxy               (Proxy(Proxy))+import           Foreign.C.String         (CString, peekCString)+import           Foreign.C.Types          (CInt(CInt), CFloat(CFloat), CDouble(CDouble), CChar)+import           Foreign.ForeignPtr       (ForeignPtr, castForeignPtr, withForeignPtr)+import           Foreign.Ptr              (Ptr, castPtr, nullPtr, plusPtr)+import           Foreign.Storable         (Storable(sizeOf, alignment, poke, peek, peekByteOff, peekElemOff, pokeByteOff, pokeElemOff))+import           GHC.TypeLits             (natVal, KnownNat, Nat)+import           System.IO.Unsafe         (unsafeDupablePerformIO)+import qualified Data.Vector.Storable     as VS+import qualified Data.ByteString          as BS+import qualified Data.ByteString.Internal as BSI++--------------------------------------------------------------------------------++-- | Like 'Proxy', but specialised to 'Nat'.+data Row (r :: Nat) = Row+-- | Like 'Proxy', but specialised to 'Nat'.+data Col (c :: Nat) = Col++-- | Used internally. Given a 'KnownNat' constraint, turn the type-level 'Nat' into an 'Int'.+natToInt :: forall n. KnownNat n => Int+{-# INLINE natToInt #-}+natToInt = fromIntegral (natVal @n Proxy)++--------------------------------------------------------------------------------++-- | Cast to and from a C-FFI type+--   'Cast' is a closed typeclass with an associated injective type family.+--   It is closed in the sense that we provide only four types+--   with instances for it; and intend for eigen to only be used+--   with those four types. The injectivity of the type family is+--   then useful for avoiding MPTCs. 'Cast' has two functions; 'toC'+--   and 'fromC', where 'toC' goes from a Haskell type to its associated+--   C type for internal use, with the C FFI, and 'fromC' goes from the+--   associated C type to the Haskell type.+class Cast (a :: Type) where+  type family C a = (result :: Type) | result -> a+  toC   :: a -> C a+  fromC :: C a -> a++instance Cast Int where+  type C Int = CInt+  toC = CInt . fromIntegral+  {-# INLINE toC #-}+  fromC (CInt x) = fromIntegral x+  {-# INLINE fromC #-}++instance Cast Float where+  type C Float = CFloat+  toC = CFloat+  {-# INLINE toC #-}+  fromC (CFloat x) = x+  {-# INLINE fromC #-}++instance Cast Double where+  type C Double = CDouble+  toC = CDouble+  {-# INLINE toC #-}+  fromC (CDouble x) = x+  {-# INLINE fromC #-}++instance Cast a => Cast (Complex a) where+  type C (Complex a) = CComplex (C a)+  toC (a :+ b) = CComplex (toC a) (toC b)+  {-# INLINE toC #-}+  fromC (CComplex a b) = (fromC a) :+ (fromC b)+  {-# INLINE fromC #-}++-- | WARNING! 'toC' is lossy for any Int greater than (maxBound :: Int32)!+instance Cast a => Cast (Int, Int, a) where+  type C (Int, Int, a) = CTriplet a+  {-# INLINE toC #-}+  toC (x, y, z) = CTriplet (toC x) (toC y) (toC z)+  {-# INLINE fromC #-}+  fromC (CTriplet x y z) = (fromC x, fromC y, fromC z)++--------------------------------------------------------------------------------++-- | Complex number for FFI with the same memory layout as std::complex\<T\>+data CComplex a = CComplex !a !a deriving (Show)++instance Storable a => Storable (CComplex a) where+    sizeOf _ = sizeOf (undefined :: a) * 2+    alignment _ = alignment (undefined :: a)+    poke p (CComplex x y) = do+        pokeElemOff (castPtr p) 0 x+        pokeElemOff (castPtr p) 1 y+    peek p = CComplex+        <$> peekElemOff (castPtr p) 0+        <*> peekElemOff (castPtr p) 1++--------------------------------------------------------------------------------++-- | FIXME: Doc+data CTriplet a where+  CTriplet :: Cast a => !CInt -> !CInt -> !(C a) -> CTriplet a++deriving instance (Show a, Show (C a)) => Show (CTriplet a)++instance (Storable a, Elem a) => Storable (CTriplet a) where+    sizeOf _ = sizeOf (undefined :: a) + sizeOf (undefined :: CInt) * 2+    alignment _ = alignment (undefined :: CInt)+    poke p (CTriplet row col val) = do+        pokeElemOff (castPtr p) 0 row+        pokeElemOff (castPtr p) 1 col+        pokeByteOff p (sizeOf (undefined :: CInt) * 2) val+    peek p = CTriplet+        <$> peekElemOff (castPtr p) 0+        <*> peekElemOff (castPtr p) 1+        <*> peekByteOff p (sizeOf (undefined :: CInt) * 2)++--------------------------------------------------------------------------------++-- | `Elem` is a closed typeclass that encompasses the properties+--   eigen expects its values to possess, and simplifies the external+--   API quite a bit.+class (Num a, Cast a, Storable a, Storable (C a), Code (C a)) => Elem a++instance Elem Float+instance Elem Double+instance Elem (Complex Float)+instance Elem (Complex Double)++--------------------------------------------------------------------------------++-- | Encode a C Type as a CInt+--+--   Hack used in FFI wrapper functions when constructing FFI calls+class Code a where; code :: a -> CInt+instance Code CFloat             where; code _ = 0+instance Code CDouble            where; code _ = 1+instance Code (CComplex CFloat)  where; code _ = 2+instance Code (CComplex CDouble) where; code _ = 3++-- | Hack used in constructing FFI calls.+newtype MagicCode = MagicCode CInt deriving Eq++instance Binary MagicCode where+    put (MagicCode _code) = putWord32be $ fromIntegral _code+    get = MagicCode . fromIntegral <$> getWord32be++-- | Hack used in constructing FFI calls.+magicCode :: Code a => a -> MagicCode+magicCode x = MagicCode (code x `xor` 0x45696730)++--------------------------------------------------------------------------------++-- | Machine size of a 'CInt'.+intSize :: Int+intSize = sizeOf (undefined :: CInt)++-- | FIXME: Doc+encodeInt :: CInt -> BS.ByteString+encodeInt x = BSI.unsafeCreate (sizeOf x) $ (`poke` x) . castPtr++-- | FIXME: Doc+decodeInt :: BS.ByteString -> CInt+decodeInt (BSI.PS fp fo fs)+    | fs == sizeOf x = x+    | otherwise = error "decodeInt: wrong buffer size"+    where x = performIO $ withForeignPtr fp $ peek . (`plusPtr` fo)++--------------------------------------------------------------------------------++-- | 'Binary' instance for 'Data.Vector.Storable.Mutable.Vector'+instance Storable a => Binary (VS.Vector a) where+    put vs = put (BS.length bs) >> putByteString bs where+        (fp,fs) = VS.unsafeToForeignPtr0 vs+        es = sizeOf (VS.head vs)+        bs = BSI.fromForeignPtr (castForeignPtr fp) 0 (fs * es)+        +    get = get >>= getByteString >>= \bs -> let+        (fp,fo,fs) = BSI.toForeignPtr bs+        es = sizeOf (VS.head vs)+        -- `plusForeignPtr` is used qualified here to just remind a reader+        -- that it is defined internally within eigen+        vs = VS.unsafeFromForeignPtr0 (Eigen.Internal.plusForeignPtr fp fo) (fs `div` es)+        in return vs++--------------------------------------------------------------------------------++-- | FIXME: Doc+data CSparseMatrix a+-- | FIXME: Doc+type CSparseMatrixPtr a = Ptr (CSparseMatrix a)++-- | FIXME: Doc+data CSolver a+-- | FIXME: Doc+type CSolverPtr a = Ptr (CSolver a)++-- {-# INLINE unholyPerformIO #-}+-- unholyPerformIO :: IO a -> a+-- unholyPerformIO (IO m) = case m realWorld# of (# _, r #) -> r++-- | FIXME: replace with unholyPerformIO (?)+performIO :: IO a -> a+performIO = unsafeDupablePerformIO++-- | FIXME: Doc+plusForeignPtr :: ForeignPtr a -> Int -> ForeignPtr b+plusForeignPtr fp fo = castForeignPtr fp1 where+    vs :: VS.Vector CChar+    vs = VS.unsafeFromForeignPtr (castForeignPtr fp) fo 0+    (fp1, _) = VS.unsafeToForeignPtr0 vs++foreign import ccall "eigen-proxy.h free" c_freeString :: CString -> IO ()++call :: IO CString -> IO ()+call func = func >>= \c_str -> when (c_str /= nullPtr) $+    peekCString c_str >>= \str -> c_freeString c_str >> fail str++foreign import ccall "eigen-proxy.h free" free :: Ptr a -> IO ()++foreign import ccall "eigen-proxy.h eigen_setNbThreads" c_setNbThreads :: CInt -> IO ()+foreign import ccall "eigen-proxy.h eigen_getNbThreads" c_getNbThreads :: IO CInt++--------------------------------------------------------------------------------++#let api1 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall a . Code (C a) => %s\n%s = c_%s (code (undefined :: (C a)))", #name, #name, args, #name, args, #name, #name++#api1 random,        "Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 identity,      "Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 add,           "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 sub,           "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 mul,           "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 diagonal,      "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 transpose,     "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 inverse,       "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 adjoint,       "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 conjugate,     "Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 normalize,     "Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 sum,           "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 prod,          "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 mean,          "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 norm,          "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 trace,         "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 squaredNorm,   "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 blueNorm,      "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 hypotNorm,     "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 determinant,   "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 rank,          "CInt -> Ptr CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 image,         "CInt -> Ptr (Ptr (C a)) -> Ptr CInt -> Ptr CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 kernel,        "CInt -> Ptr (Ptr (C a)) -> Ptr CInt -> Ptr CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 solve,         "CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api1 relativeError, "Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString"++--------------------------------------------------------------------------------++#let api2 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> %s\n%s :: forall a . Code (C a) => %s\n%s = c_%s (code (undefined :: (C a)))", #name, #name, args, #name, args, #name, #name++#api2 sparse_new,           "CInt -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_clone,         "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_fromList,      "CInt -> CInt -> Ptr (CTriplet a) -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_toList,        "CSparseMatrixPtr a -> Ptr (CTriplet a) -> CInt -> IO CString"+#api2 sparse_free,          "CSparseMatrixPtr a -> IO CString"+#api2 sparse_makeCompressed,"CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_uncompress,    "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_isCompressed,  "CSparseMatrixPtr a -> Ptr CInt -> IO CString"+#api2 sparse_transpose,     "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_adjoint,       "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_pruned,        "CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_prunedRef,     "CSparseMatrixPtr a -> Ptr (C a) -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_scale,         "CSparseMatrixPtr a -> Ptr (C a) -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_nonZeros,      "CSparseMatrixPtr a -> Ptr CInt -> IO CString"+#api2 sparse_innerSize,     "CSparseMatrixPtr a -> Ptr CInt -> IO CString"+#api2 sparse_outerSize,     "CSparseMatrixPtr a -> Ptr CInt -> IO CString"+#api2 sparse_coeff,         "CSparseMatrixPtr a -> CInt -> CInt -> Ptr (C a) -> IO CString"+#api2 sparse_coeffRef,      "CSparseMatrixPtr a -> CInt -> CInt -> Ptr (Ptr (C a)) -> IO CString"+#api2 sparse_cols,          "CSparseMatrixPtr a -> Ptr CInt -> IO CString"+#api2 sparse_rows,          "CSparseMatrixPtr a -> Ptr CInt -> IO CString"+#api2 sparse_norm,          "CSparseMatrixPtr a -> Ptr (C a) -> IO CString"+#api2 sparse_squaredNorm,   "CSparseMatrixPtr a -> Ptr (C a) -> IO CString"+#api2 sparse_blueNorm,      "CSparseMatrixPtr a -> Ptr (C a) -> IO CString"+#api2 sparse_add,           "CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_sub,           "CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_mul,           "CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_block,         "CSparseMatrixPtr a -> CInt -> CInt -> CInt -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_fromMatrix,    "Ptr (C a) -> CInt -> CInt -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api2 sparse_toMatrix,      "CSparseMatrixPtr a -> Ptr (C a) -> CInt -> CInt -> IO CString"+#api2 sparse_values,        "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr (C a)) -> IO CString"+#api2 sparse_outerStarts,   "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"+#api2 sparse_innerIndices,  "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"+#api2 sparse_innerNNZs,     "CSparseMatrixPtr a -> Ptr CInt -> Ptr (Ptr CInt) -> IO CString"+#api2 sparse_setZero,       "CSparseMatrixPtr a -> IO CString"+#api2 sparse_setIdentity,   "CSparseMatrixPtr a -> IO CString"+#api2 sparse_reserve,       "CSparseMatrixPtr a -> CInt -> IO CString"+#api2 sparse_resize,        "CSparseMatrixPtr a -> CInt -> CInt -> IO CString"++#api2 sparse_conservativeResize,    "CSparseMatrixPtr a -> CInt -> CInt -> IO CString"+#api2 sparse_compressInplace,       "CSparseMatrixPtr a -> IO CString"+#api2 sparse_uncompressInplace,     "CSparseMatrixPtr a -> IO CString"++--------------------------------------------------------------------------------++#let api3 name, args = "foreign import ccall \"eigen_%s\" c_%s :: CInt -> CInt -> %s\n%s :: forall s a . (Code s, Code (C a)) => s -> %s\n%s s = c_%s (code (undefined :: (C a))) (code s)", #name, #name, args, #name, args, #name, #name++#api3 sparse_la_newSolver,          "Ptr (CSolverPtr a) -> IO CString"+#api3 sparse_la_freeSolver,         "CSolverPtr a -> IO CString"+#api3 sparse_la_factorize,          "CSolverPtr a -> CSparseMatrixPtr a -> IO CString"+#api3 sparse_la_analyzePattern,     "CSolverPtr a -> CSparseMatrixPtr a -> IO CString"+#api3 sparse_la_compute,            "CSolverPtr a -> CSparseMatrixPtr a -> IO CString"+#api3 sparse_la_tolerance,          "CSolverPtr a -> Ptr CDouble -> IO CString"+#api3 sparse_la_setTolerance,       "CSolverPtr a -> CDouble -> IO CString"+#api3 sparse_la_maxIterations,      "CSolverPtr a -> Ptr CInt -> IO CString"+#api3 sparse_la_setMaxIterations,   "CSolverPtr a -> CInt -> IO CString"+#api3 sparse_la_info,               "CSolverPtr a -> Ptr CInt -> IO CString"+#api3 sparse_la_error,              "CSolverPtr a -> Ptr CDouble -> IO CString"+#api3 sparse_la_iterations,         "CSolverPtr a -> Ptr CInt -> IO CString"+#api3 sparse_la_solve,              "CSolverPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+-- #api3 sparse_la_solveWithGuess,     "CSolverPtr a -> CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api3 sparse_la_matrixQ,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api3 sparse_la_matrixR,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api3 sparse_la_setPivotThreshold,  "CSolverPtr a -> CDouble -> IO CString"+#api3 sparse_la_rank,               "CSolverPtr a -> Ptr CInt -> IO CString"+#api3 sparse_la_matrixL,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api3 sparse_la_matrixU,            "CSolverPtr a -> Ptr (CSparseMatrixPtr a) -> IO CString"+#api3 sparse_la_setSymmetric,       "CSolverPtr a -> CInt -> IO CString"+#api3 sparse_la_determinant,        "CSolverPtr a -> Ptr (C a) -> IO CString"+#api3 sparse_la_logAbsDeterminant,  "CSolverPtr a -> Ptr (C a) -> IO CString"+#api3 sparse_la_absDeterminant,     "CSolverPtr a -> Ptr (C a) -> IO CString"+#api3 sparse_la_signDeterminant,    "CSolverPtr a -> Ptr (C a) -> IO CString"++--------------------------------------------------------------------------------

@@ -0,0 +1,544 @@+{-# LANGUAGE BangPatterns #-}+{-# LANGUAGE ConstraintKinds #-}+{-# LANGUAGE CPP #-}+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE GADTs #-}+{-# LANGUAGE KindSignatures #-}+{-# LANGUAGE LambdaCase #-}+{-# LANGUAGE Rank2Types #-}+{-# LANGUAGE ScopedTypeVariables #-}+#if __GLASGOW_HASKELL__ >= 805+{-# LANGUAGE ExplicitNamespaces #-}+{-# LANGUAGE NoStarIsType #-}+#endif+{-# LANGUAGE TypeApplications #-}+{-# LANGUAGE TypeFamilies #-}+{-# LANGUAGE TypeOperators #-}++module Eigen.Matrix+  ( +    -- * Types +    Matrix(..)+  , Vec(..)+  , MatrixXf+  , MatrixXd+  , MatrixXcf+  , MatrixXcd++    -- * Common API+  , Elem+  , C+  , natToInt+  , Row(..)+  , Col(..)++    -- * Encode/Decode a Matrix+  , encode+  , decode++    -- * Querying a Matrix+  , null+  , square+  , rows+  , cols+  , dims+    +    -- * Constructing a Matrix+  , empty+  , constant+  , zero+  , ones+  , identity+  , random+  , diagonal++  , (!)+  , coeff+  , generate+  , sum+  , prod+  , mean+  , trace+  , all+  , any+  , count+  , norm+  , squaredNorm+  , blueNorm+  , hypotNorm+  , determinant+  , add+  , sub+  , mul+  , map+  , imap+  , TriangularMode(..)+  , triangularView+  , filter+  , ifilter+  , length+  , foldl+  , foldl'+  , inverse+  , adjoint+  , transpose+  , conjugate+  , normalize+  , modify+  , block+  , unsafeFreeze+  , unsafeWith+  , fromList+  , toList+  ) where++import Control.Monad (when)+import Control.Monad.ST (ST)+import Prelude hiding+  (map, null, filter, length, foldl, any, all, sum)+import Control.Monad (forM_)+import Control.Monad.Primitive (PrimMonad(..))+import Data.Binary (Binary(..))+import qualified Data.Binary as Binary+import qualified Data.ByteString.Lazy as BSL+import Data.Complex (Complex)+import Data.Constraint.Nat+import Eigen.Internal+  ( Elem+  , Cast(..)+  , natToInt+  , Row(..)+  , Col(..)+  )+import qualified Eigen.Internal as Internal+import qualified Eigen.Matrix.Mutable as M+import qualified Data.List as List+import Data.Kind (Type)+import GHC.TypeLits (Nat, type (*), type (<=), KnownNat)+import Foreign.C.Types (CInt)+import Foreign.C.String (CString)+import Foreign.Marshal.Alloc (alloca)+import Foreign.Ptr (Ptr)+import Foreign.Storable (peek)++import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VSM++-- | Matrix to be used in pure computations.+--+--   * Uses column majour memory layout.+--+--   * Has a copy-free FFI using the <http://eigen.tuxfamily.org Eigen> library.+--+newtype Matrix :: Nat -> Nat -> Type -> Type where+  Matrix :: Vec (n * m) a -> Matrix n m a++-- | Used internally to track the size and corresponding C type of the matrix.+newtype Vec :: Nat -> Type -> Type where+  Vec :: VS.Vector (C a) -> Vec n a++instance forall n m a. (Elem a, Show a, KnownNat n, KnownNat m) => Show (Matrix n m a) where+  show m = List.concat+    [ "Matrix ", show (rows m), "x", show (cols m)+    , "\n", List.intercalate "\n" $ List.map (List.intercalate "\t" . List.map show) $ toList m, "\n"+    ]++instance forall n m a. (KnownNat n, KnownNat m, Elem a) => Binary (Matrix n m a) where+  put (Matrix (Vec vals)) = do+    put $ Internal.magicCode (undefined :: C a)+    put $ natToInt @n+    put $ natToInt @m+    put vals++  get = do+    get >>= (`when` fail "wrong matrix type") . (/= Internal.magicCode (undefined :: C a))+    Matrix . Vec <$> get++-- | Encode the sparse matrix as a lazy bytestring+encode :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> BSL.ByteString+encode = Binary.encode++-- | Decode the sparse matrix from a lazy bytestring+decode :: (Elem a, KnownNat n, KnownNat m) => BSL.ByteString -> Matrix n m a+decode = Binary.decode++-- | Alias for single precision matrix+type MatrixXf n m = Matrix n m Float+-- | Alias for double precision matrix+type MatrixXd n m = Matrix n m Double+-- | Alias for single precision matrix of complex numbers+type MatrixXcf n m = Matrix n m (Complex Float)+-- | Alias for double precision matrix of complex numbers+type MatrixXcd n m = Matrix n m (Complex Double)++-- | Construct an empty 0x0 matrix+empty :: Elem a => Matrix 0 0 a+{-# INLINE empty #-}+empty = Matrix (Vec (VS.empty))++-- | Is matrix empty?+null :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Bool+{-# INLINE null #-}+null m = cols m == 0 && rows m == 0++-- | Is matrix square?+--+square :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Bool+{-# INLINE square #-}+square _ = natToInt @n == natToInt @m++-- | Matrix where all coeffs are filled with the given value+constant :: forall n m a. (Elem a, KnownNat n, KnownNat m) => a -> Matrix n m a+{-# INLINE constant #-}+constant !val =+  let !cval = toC val+  in withDims $ \rs cs -> VS.replicate (rs * cs) cval++-- | Matrix where all coeffs are filled with 0+zero :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a+{-# INLINE zero #-}+zero = constant 0++-- | Matrix where all coeffs are filled with 1+ones :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a+{-# INLINE ones #-}+ones = constant 1++-- | The identity matrix (not necessarily square)+identity :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a+identity =+  Internal.performIO $ do+     m :: M.IOMatrix n m a <- M.new+     Internal.call $ M.unsafeWith m Internal.identity+     unsafeFreeze m++-- | The random matrix of a given size+random :: forall n m a. (Elem a, KnownNat n, KnownNat m) => IO (Matrix n m a)+random = do+  m :: M.IOMatrix n m a <- M.new+  Internal.call $ M.unsafeWith m Internal.random+  unsafeFreeze m++withDims :: forall n m a. (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> VS.Vector (C a)) -> Matrix n m a+{-# INLINE withDims #-}+withDims f =+  let !r = natToInt @n+      !c = natToInt @m+  in Matrix $ Vec $ f r c++-- | The number of rows in the matrix+rows :: forall n m a. KnownNat n => Matrix n m a -> Int+{-# INLINE rows #-}+rows _ = natToInt @n++-- | The number of colums in the matrix+cols :: forall n m a. KnownNat m => Matrix n m a -> Int+{-# INLINE cols #-}+cols _ = natToInt @m++-- | Return Matrix size as a pair of (rows, cols)+dims :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> (Int, Int)+{-# INLINE dims #-}+dims _ = (natToInt @n, natToInt @m)++-- | Return the value at the given position.+(!) :: forall n m a r c. (Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m) => Row r -> Col c -> Matrix n m a -> a+{-# INLINE (!) #-}+(!) = coeff++-- | Return the value at the given position.+coeff :: forall n m a r c. (Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m) => Row r -> Col c -> Matrix n m a -> a+{-# INLINE coeff #-}+coeff _ _ m@(Matrix (Vec vals)) =+  let !row  = natToInt @r+      !col  = natToInt @c+  in fromC $! VS.unsafeIndex vals $! col * rows m + row++unsafeCoeff :: (Elem a, KnownNat n) => Int -> Int -> Matrix n m a -> a+{-# INLINE unsafeCoeff #-}+unsafeCoeff row col m@(Matrix (Vec vals)) = fromC $! VS.unsafeIndex vals $! col * rows m + row++-- | Given a generation function `f :: Int -> Int -> a`, construct a Matrix of known size+--   using points in the matrix as inputs.+generate :: forall n m a. (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> a) -> Matrix n m a+generate f = withDims $ \rs cs -> VS.create $ do+  vals :: VSM.MVector s (C a) <- VSM.new (rs * cs)+  forM_ [0 .. pred rs] $ \r ->+    forM_ [0 .. pred cs] $ \c ->+      VSM.write vals (c * rs + r) (toC $! f r c)+  pure vals++-- | The sum of all coefficients in the matrix+sum :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a+sum = _prop Internal.sum++-- | The product of all coefficients in the matrix+prod :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a+prod = _prop Internal.prod++-- | The arithmetic mean of all coefficients in the matrix+mean :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a+mean = _prop Internal.mean++-- | The trace of a matrix is the sum of the diagonal coefficients.+--   +--   'trace' m == 'sum' ('diagonal' m)+trace :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a+trace = _prop Internal.trace++-- | Given a predicate p, determine if all values in the Matrix satisfy p.+all :: (Elem a, KnownNat n, KnownNat m) => (a -> Bool) -> Matrix n m a -> Bool+all f (Matrix (Vec vals)) = VS.all (f . fromC) vals++-- | Given a predicate p, determine if any values in the Matrix satisfy p.+any :: (Elem a, KnownNat n, KnownNat m) => (a -> Bool) -> Matrix n m a -> Bool+any f (Matrix (Vec vals)) = VS.any (f . fromC) vals++-- | Given a predicate p, determine how many values in the Matrix satisfy p.+count :: (Elem a, KnownNat n, KnownNat m) => (a -> Bool) -> Matrix n m a -> Int+count f (Matrix (Vec vals)) = VS.foldl' (\n x-> if f (fromC x) then (n + 1) else n) 0 vals++norm, squaredNorm, blueNorm, hypotNorm :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> a++{-| For vectors, the l2 norm, and for matrices the Frobenius norm.+    In both cases, it consists in the square root of the sum of the square of all the matrix entries.+    For vectors, this is also equals to the square root of the dot product of this with itself.+-}+norm = _prop Internal.norm++-- | For vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.+squaredNorm = _prop Internal.squaredNorm++-- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.+blueNorm = _prop Internal.blueNorm++-- | The l2 norm of the matrix avoiding undeflow and overflow. This version use a concatenation of hypot calls, and it is very slow.+hypotNorm = _prop Internal.hypotNorm++-- | The determinant of the matrix+determinant :: forall n a. (Elem a, KnownNat n) => Matrix n n a -> a+determinant m = _prop Internal.determinant m++-- | Add two matrices.+add :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a -> Matrix n m a+add m1 m2 = _binop Internal.add m1 m2++-- | Subtract two matrices.+sub :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a -> Matrix n m a+sub m1 m2 = _binop Internal.sub m1 m2++-- | Multiply two matrices.+mul :: (Elem a, KnownNat p, KnownNat q, KnownNat r) => Matrix p q a -> Matrix q r a -> Matrix p r a+mul m1 m2 = _binop Internal.mul m1 m2++{- | Apply a given function to each element of the matrix.+Here is an example how to implement scalar matrix multiplication:+>>> let a = fromList [[1,2],[3,4]] :: MatrixXf 2 2+>>> a+Matrix 2x2+1.0 2.0+3.0 4.0+>>> map (*10) a+Matrix 2x2+10.0    20.0+30.0    40.0+-}+map :: Elem a => (a -> a) -> Matrix n m a -> Matrix n m a+map f (Matrix (Vec vals)) = Matrix $ Vec $ VS.map (toC . f . fromC) vals++{- | Apply a given function to each element of the matrix.+Here is an example how upper triangular matrix can be implemented:+>>> let a = fromList [[1,2,3],[4,5,6],[7,8,9]] :: MatrixXf+>>> a+Matrix 3x3+1.0 2.0 3.0+4.0 5.0 6.0+7.0 8.0 9.0+>>> imap (\row col val -> if row <= col then val else 0) a+Matrix 3x3+1.0 2.0 3.0+0.0 5.0 6.0+0.0 0.0 9.0+-}+imap :: (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> a -> a) -> Matrix n m a -> Matrix n m a+imap f (Matrix (Vec vals)) =+  withDims $ \rs _ ->+    VS.imap (\n ->+      let (c,r) = divMod n rs+      in toC . f r c . fromC) vals++-- | Provide a view of the matrix for extraction of a subset.+data TriangularMode+  -- | View matrix as a lower triangular matrix.+  = Lower+  -- | View matrix as an upper triangular matrix.+  | Upper+  -- | View matrix as a lower triangular matrix with zeros on the diagonal.+  | StrictlyLower+  -- | View matrix as an upper triangular matrix with zeros on the diagonal.+  | StrictlyUpper+  -- | View matrix as a lower triangular matrix with ones on the diagonal.+  | UnitLower+  -- | View matrix as an upper triangular matrix with ones on the diagonal.+  | UnitUpper+  deriving (Eq, Enum, Show, Read)++-- | Triangular view extracted from the current matrix+triangularView :: (Elem a, KnownNat n, KnownNat m) => TriangularMode -> Matrix n m a -> Matrix n m a+triangularView = \case+  Lower         -> imap $ \row col val -> case compare row col of { LT -> 0; _ -> val }+  Upper         -> imap $ \row col val -> case compare row col of { GT -> 0; _ -> val }+  StrictlyLower -> imap $ \row col val -> case compare row col of { GT -> val; _ -> 0 }+  StrictlyUpper -> imap $ \row col val -> case compare row col of { LT -> val; _ -> 0 }+  UnitLower     -> imap $ \row col val -> case compare row col of { GT -> val; LT -> 0; EQ -> 1 }+  UnitUpper     -> imap $ \row col val -> case compare row col of { LT -> val; GT -> 0; EQ -> 1 }++-- | Filter elements in the matrix. Filtered elements will be replaced by 0.+filter :: Elem a => (a -> Bool) -> Matrix n m a -> Matrix n m a+filter f = map (\x -> if f x then x else 0)++-- | Filter elements in the matrix with an indexed predicate. Filtered elements will be replaces by 0.+ifilter :: (Elem a, KnownNat n, KnownNat m) => (Int -> Int -> a -> Bool) -> Matrix n m a -> Matrix n m a+ifilter f = imap (\r c x -> if f r c x then x else 0)++-- | The length of the matrix.+length :: forall n m a r. (Elem a, KnownNat n, KnownNat m, r ~ (n * m), KnownNat r) => Matrix n m a -> Int+length _ = natToInt @r++-- | Left fold of a matrix, where accumulation is lazy.+foldl :: (Elem a, KnownNat n, KnownNat m) => (b -> a -> b) -> b -> Matrix n m a -> b+foldl f b (Matrix (Vec vals)) = VS.foldl (\a x -> f a (fromC x)) b vals++-- | Right fold of a matrix, where accumulation is strict.+foldl' :: Elem a => (b -> a -> b) -> b -> Matrix n m a -> b+foldl' f b (Matrix (Vec vals)) = VS.foldl' (\ !a x -> f a (fromC x)) b vals++-- | Return the diagonal of a matrix.+diagonal :: (Elem a, KnownNat n, KnownNat m, r ~ Min n m, KnownNat r) => Matrix n m a -> Matrix r 1 a+diagonal = _unop Internal.diagonal++{- | Inverse of the matrix+For small fixed sizes up to 4x4, this method uses cofactors. In the general case, this method uses PartialPivLU decomposition+-}+inverse :: forall n a. (Elem a, KnownNat n) => Matrix n n a -> Matrix n n a+inverse = _unop Internal.inverse++-- | Adjoint of the matrix+adjoint :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix m n a+adjoint = _unop Internal.adjoint++-- | Transpose of the matrix+transpose :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix m n a+transpose = _unop Internal.transpose++-- | Conjugate of the matrix+conjugate :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a+conjugate = _unop Internal.conjugate++-- | Normalise the matrix by dividing it on its 'norm'+normalize :: forall n m a. (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> Matrix n m a+normalize (Matrix (Vec vals)) = Internal.performIO $ do+  vals' <- VS.thaw vals+  VSM.unsafeWith vals' $ \p ->+    let !rs = natToInt @n+        !cs = natToInt @m+    in Internal.call $ Internal.normalize p (toC rs) (toC cs)+  Matrix . Vec <$> VS.unsafeFreeze vals'++-- | Apply a destructive operation to a matrix. The operation will be performed in-place, if it is safe+--   to do so - otherwise, it will create a copy of the matrix.+modify :: (Elem a, KnownNat n, KnownNat m) => (forall s. M.MMatrix n m s a -> ST s ()) -> Matrix n m a -> Matrix n m a+modify f (Matrix (Vec vals)) = Matrix $ Vec $ VS.modify (f . M.fromVector ) vals++-- | Extract rectangular block from matrix defined by startRow startCol blockRows blockCols+block :: forall sr sc br bc n m a.+     (Elem a, KnownNat sr, KnownNat sc, KnownNat br, KnownNat bc, KnownNat n, KnownNat m)+  => (sr <= n, sc <= m, br <= n, bc <= m)+  => Row sr -- ^ starting row+  -> Col sc -- ^ starting col+  -> Row br -- ^ block of rows+  -> Col bc -- ^ block of cols+  -> Matrix n m a -- ^ extract from this+  -> Matrix br bc a -- ^ extraction+block _ _ _ _ m =+  let !startRow = natToInt @sr+      !startCol = natToInt @sc+  in generate $ \row col -> unsafeCoeff (startRow + row) (startCol + col) m++-- | Turn a mutable matrix into an immutable matrix without copying.+--   The mutable matrix should not be modified after this conversion.+unsafeFreeze :: (Elem a, KnownNat n, KnownNat m, PrimMonad p) => M.MMatrix n m (PrimState p) a -> p (Matrix n m a)+unsafeFreeze m = VS.unsafeFreeze (M.vals m) >>= pure . Matrix . Vec+  +-- | Pass a pointer to the matrix's data to the IO action. The data may not be modified through the pointer.+unsafeWith  :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> (Ptr (C a) -> CInt -> CInt -> IO b) -> IO b+unsafeWith m@(Matrix (Vec (vals))) f =+  VS.unsafeWith vals $ \p ->+    let !rs = toC $! rows m+        !cs = toC $! cols m+    in f p rs cs++_prop :: (Elem a, KnownNat n, KnownNat m) => (Ptr (C a) -> Ptr (C a) -> CInt -> CInt -> IO CString) -> Matrix n m a -> a+{-# INLINE _prop #-}+_prop f m = fromC $ Internal.performIO $ alloca $ \p -> do+   Internal.call $ unsafeWith m (f p)+   peek p++_binop :: forall n m n1 m1 n2 m2 a. (Elem a, KnownNat n, KnownNat m, KnownNat n1, KnownNat m1, KnownNat n2, KnownNat m2)+  => (Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString)+  -> Matrix n m a+  -> Matrix n1 m1 a+  -> Matrix n2 m2 a+{-# INLINE _binop #-}+_binop g m1 m2 = Internal.performIO $ do+  m0 :: M.IOMatrix n2 m2 a <- M.new+  M.unsafeWith m0 $ \vals0 rows0 cols0 ->+      unsafeWith m1 $ \vals1 rows1 cols1 ->+          unsafeWith m2 $ \vals2 rows2 cols2 ->+              Internal.call $ g+                vals0 rows0 cols0+                vals1 rows1 cols1+                vals2 rows2 cols2+  unsafeFreeze m0++_unop :: forall n m n1 m1 a. (Elem a, KnownNat n, KnownNat m, KnownNat n1, KnownNat m1)+  => (Ptr (C a) -> CInt -> CInt -> Ptr (C a) -> CInt -> CInt -> IO CString)+  -> Matrix n m a+  -> Matrix n1 m1 a+{-# INLINE _unop #-}+_unop g m1 = Internal.performIO $ do+  m0 :: M.IOMatrix n1 m1 a <- M.new+  M.unsafeWith m0 $ \vals0 rows0 cols0 ->+      unsafeWith m1 $ \vals1 rows1 cols1 ->+          Internal.call $ g+              vals0 rows0 cols0+              vals1 rows1 cols1+  unsafeFreeze m0++-- | Convert a matrix to a list.+toList :: (Elem a, KnownNat n, KnownNat m) => Matrix n m a -> [[a]]+{-# INLINE toList #-}+toList m@(Matrix (Vec vals))+  | null m = []+  | otherwise = [[fromC $ vals `VS.unsafeIndex` (col * _rows + row) | col <- [0..pred _cols]] | row <- [0..pred _rows]]+  where+    !_rows = rows m+    !_cols = cols m++-- | Convert a list to a matrix. Returns 'Nothing' if the dimensions of the list do not match that+--   of the matrix.+fromList :: forall n m a. (Elem a, KnownNat n, KnownNat m) => [[a]] -> Maybe (Matrix n m a)+fromList list = do+  let myRows = natToInt @n+  let myCols = natToInt @m+  let _rows  = List.length list+  let _cols  = List.foldl' max 0 (List.map List.length list)+  if ((myRows /= _rows) || (myCols /= _cols))+    then Nothing+    else (Just . Matrix . Vec) $ VS.create $ do+      vm <- VSM.replicate (_rows * _cols) (toC (0 :: a))+      forM_ (zip [0..] list) $ \(row,vals) ->+        forM_ (zip [0..] vals) $ \(col, val) ->+          VSM.write vm (col * _rows + row) (toC val)+      pure vm

@@ -0,0 +1,154 @@+{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE CPP                 #-}+{-# LANGUAGE DataKinds           #-}+{-# LANGUAGE GADTs               #-}+{-# LANGUAGE KindSignatures      #-}+{-# LANGUAGE RecordWildCards     #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications    #-}+{-# LANGUAGE TypeInType          #-}+{-# LANGUAGE TypeOperators       #-}++#if __GLASGOW_HASKELL__ >= 805+{-# LANGUAGE NoStarIsType #-}+{-# LANGUAGE ExplicitNamespaces #-}+#endif++module Eigen.Matrix.Mutable+  ( -- * Types+    MMatrix(..)+  , MMatrixXf+  , MMatrixXd+  , MMatrixXcf+  , MMatrixXcd+  , IOMatrix+  , STMatrix++    -- * Construction+  , new+  , replicate++    -- * Indexing+  , read+  , write++    -- * Modification+  , set+  , copy++    -- * Modification with Pointers+  , unsafeWith+  +    -- * Conversion to Vectors+  , vals+  , fromVector+  ) where++import Eigen.Internal+  ( Elem+  , C(..)+  , natToInt+  , Row(..)+  , Col(..)+  )++import Prelude hiding (replicate, read)+import Control.Monad.Primitive (PrimMonad(..))+import Data.Complex (Complex)+import Data.Kind (Type)+import Foreign.C.Types (CInt)+import Foreign.Ptr (Ptr)+import GHC.Exts (RealWorld)+import GHC.TypeLits (Nat, type (*), type (<=), KnownNat)+import qualified Data.Vector.Storable.Mutable as VSM++-- | A mutable matrix. See 'Eigen.Matrix.Matrix' for+--   details about matrix layout.+newtype MMatrix :: Nat -> Nat -> Type -> Type -> Type where+  MMatrix :: Vec (n * m) s a -> MMatrix n m s a++-- | Used internally to track the size and corresponding C type of the matrix.+newtype Vec :: Nat -> Type -> Type -> Type where+  Vec :: VSM.MVector s (C a) -> Vec n s a++-- | Alias for single precision mutable matrix+type MMatrixXf  n m s = MMatrix n m s Float+-- | Alias for double precision mutable matrix+type MMatrixXd  n m s = MMatrix n m s Double+-- | Alias for single precision mutable matrix of complex numbers+type MMatrixXcf n m s = MMatrix n m s (Complex Float)+-- | Alias for double precision mutable matrix of complex numbers+type MMatrixXcd n m s = MMatrix n m s (Complex Double)++-- | A mutable matrix where the state token is specialised to 'RealWorld'.+type IOMatrix n m a   = MMatrix n m RealWorld a+-- | This type does not differ from MSparseMatrix, but might be desirable for readability.+type STMatrix n m s a = MMatrix n m s a++-- | Create a mutable matrix of the given size and fill it with 0 as an initial value.+new :: (PrimMonad p, Elem a, KnownNat n, KnownNat m) => p (MMatrix n m (PrimState p) a)+{-# INLINE new #-}+new = replicate 0++-- | Create a mutable matrix of the given size and fill it with an initial value.+replicate :: forall n m p a. (PrimMonad p, Elem a, KnownNat n, KnownNat m) => a -> p (MMatrix n m (PrimState p) a)+{-# INLINE replicate #-}+replicate !val = do+  let !mm_rows = natToInt @n+      !mm_cols = natToInt @m+      !cval    = toC val+  _vals <- VSM.replicate (mm_rows * mm_cols) cval+  pure (MMatrix $! Vec $! _vals)++-- | Set all elements of the matrix to a given value.+set :: (PrimMonad p, Elem a) => MMatrix n m (PrimState p) a -> a -> p ()+{-# INLINE set #-}+set (MMatrix (Vec !vec)) !val =+  let !cval = toC val+  in VSM.set vec cval++-- | Copy a matrix.+copy :: (PrimMonad p, Elem a) => MMatrix n m (PrimState p) a -> MMatrix n m (PrimState p) a -> p ()+{-# INLINE copy #-}+copy (MMatrix (Vec m1)) (MMatrix (Vec m2)) = VSM.unsafeCopy m1 m2++-- | Yield the element at the given position.+read :: forall n m p a r c. (PrimMonad p, Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m)+  => Row r -> Col c -> MMatrix n m (PrimState p) a -> p a+{-# INLINE read #-}+read _ _ (MMatrix (Vec m)) =+  let !row = natToInt @r+      !col = natToInt @c+      !mm_rows = natToInt @n+  in VSM.unsafeRead m (col * mm_rows + row) >>= \ !val -> let !cval = fromC val in pure cval++-- | Replace the element at the given position.+write :: forall n m p a r c. (PrimMonad p, Elem a, KnownNat n, KnownNat r, KnownNat c, r <= n, c <= m)+  => Row r -> Col c -> MMatrix n m (PrimState p) a -> a -> p ()+{-# INLINE write #-}+write _ _ (MMatrix (Vec m)) !val =+  let !row = natToInt @r+      !col = natToInt @c+      !mm_rows = natToInt @n+      !cval = toC val+  in VSM.unsafeWrite m (col * mm_rows + row) cval++-- | Pass a pointer to the matrix's data to the IO action.+--   Modifying dat through the pointer is unsafe if the matrix+--   could have been frozen before the modification.+unsafeWith :: forall n m a b. (KnownNat n, KnownNat m, Elem a) => IOMatrix n m a -> (Ptr (C a) -> CInt -> CInt -> IO b) -> IO b+{-# INLINE unsafeWith #-}+unsafeWith (MMatrix (Vec m)) f =+  let !cmm_rows = toC $! natToInt @n+      !cmm_cols = toC $! natToInt @m+  in VSM.unsafeWith m $ \p -> f p cmm_rows cmm_cols++-- | Return a mutable storable 'VSM.MVector' of the corresponding C types to one's mutable matrix.+vals :: MMatrix n m s a -> VSM.MVector s (C a)+{-# INLINE vals #-}+vals (MMatrix (Vec x)) = x++-- | Create a mutable matrix from a mutable storable 'VSM.MVector'.+fromVector :: Elem a => VSM.MVector s (C a) -> MMatrix n m s a+{-# INLINE fromVector #-}+fromVector x = MMatrix (Vec x)

@@ -0,0 +1,28 @@+{- |+Some of Eigen's algorithms can exploit the multiple cores present in your hardware.+To this end, it is enough to enable OpenMP on your compiler, for instance: GCC: -fopenmp.+You can control the number of threads that will be used using either by the OpenMP API or by Eigen's API using the following priority:++1. OMP_NUM_THREADS=n ./my_program+2. setNbThreads n++Unless setNbThreads has been called, Eigen uses the number of threads specified by OpenMP.+You can restore this behaviour by calling @setNbThreads n@++Currently, the following algorithms can make use of multi-threading: general matrix - matrix products PartialPivLU.+-}++module Eigen.Parallel+  ( setNbThreads+  , getNbThreads+  ) where++import Eigen.Internal (Cast(..), c_setNbThreads, c_getNbThreads)++-- | Sets the max number of threads reserved for Eigen+setNbThreads :: Int -> IO ()+setNbThreads = c_setNbThreads . toC++-- | Gets the max number of threads reserved for Eigen+getNbThreads :: IO Int+getNbThreads = fromC <$> c_getNbThreads

@@ -0,0 +1,204 @@+{-# LANGUAGE DataKinds             #-}+{-# LANGUAGE GADTs                 #-}+{-# LANGUAGE FlexibleInstances     #-}+{-# LANGUAGE KindSignatures        #-}+{-# LANGUAGE LambdaCase            #-}+{-# LANGUAGE MagicHash             #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE OverloadedStrings     #-}+{-# LANGUAGE ScopedTypeVariables   #-}+{-# LANGUAGE StandaloneDeriving    #-}++module Eigen.Solver.LA+  ( Decomposition(..)+  , solve+  , relativeError+  , rank+  , kernel+  , image+  , linearRegression+  ) where++import Eigen.Internal (Elem, Cast(..))+import Eigen.Matrix+import Foreign.C.Types (CInt)+import Foreign.Marshal.Alloc (alloca)+import Foreign.Storable (Storable(..))+import GHC.TypeLits (KnownNat)+import GHC.Types+import Prelude+import qualified Eigen.Internal as Internal+import qualified Eigen.Matrix.Mutable as MM+import qualified Eigen.Matrix as M+import qualified Data.List as List+import qualified Data.Vector.Storable as VS+import qualified Foreign.Concurrent as FC+--import qualified Prelude as Prelude++{- |+@+Decomposition           Requirements on the matrix          Speed   Accuracy  Rank  Kernel  Image++PartialPivLU            Invertible                          ++      +         -     -       -+FullPivLU               None                                -       +++       +     +       ++HouseholderQR           None                                ++      +         -     -       -+ColPivHouseholderQR     None                                +       ++        +     -       -+FullPivHouseholderQR    None                                -       +++       +     -       -+LLT                     Positive definite                   +++     +         -     -       -+LDLT                    Positive or negative semidefinite   +++     ++        -     -       -+JacobiSVD               None                                -       +++       +     -       -+@+The best way to do least squares solving for square matrices is with a SVD decomposition ('JacobiSVD')+-}+data Decomposition+  -- | LU decomposition of a matrix with partial pivoting.+  = PartialPivLU+  -- | LU decomposition of a matrix with complete pivoting.+  | FullPivLU+  -- | Householder QR decomposition of a matrix.+  | HouseholderQR+  -- | Householder rank-revealing QR decomposition of a matrix with column-pivoting.+  | ColPivHouseholderQR+  -- | Householder rank-revealing QR decomposition of a matrix with full pivoting.+  | FullPivHouseholderQR+  -- | Standard Cholesky decomposition (LL^T) of a matrix.+  | LLT+  -- | Robust Cholesky decomposition of a matrix with pivoting.+  | LDLT+  -- | Two-sided Jacobi SVD decomposition of a rectangular matrix.+  | JacobiSVD+  deriving (Enum, Eq, Show, Read)++con2CTag :: Decomposition -> CInt+{-# INLINE con2CTag #-}+con2CTag = \case+  PartialPivLU         -> 0+  FullPivLU            -> 1+  HouseholderQR        -> 2+  ColPivHouseholderQR  -> 3+  FullPivHouseholderQR -> 4+  LLT                  -> 5+  LDLT                 -> 6+  JacobiSVD            -> 7++-- | [x = solve d a b] finds a solution @x@ of @ax = b@ equation using decomposition @d@+solve :: (KnownNat n, KnownNat m, KnownNat n1, KnownNat m1, Elem a)+  => Decomposition+  -> Matrix n m a+  -> Matrix n1 m1 a+  -> Matrix m 1 a+solve d a b = Internal.performIO $ do+  x :: MM.IOMatrix m 1 a <- MM.new+  MM.unsafeWith x $ \x_vals x_rows x_cols ->+      unsafeWith a $ \a_vals a_rows a_cols ->+          unsafeWith b $ \b_vals b_rows b_cols ->+              Internal.call $ Internal.solve (con2CTag d)+                x_vals x_rows x_cols+                a_vals a_rows a_cols+                b_vals b_rows b_cols+  unsafeFreeze x++-- | [e = relativeError x a b] computes @norm (ax - b) / norm b@ where @norm@ is L2 norm+relativeError :: (KnownNat n, KnownNat m, KnownNat n1, KnownNat m1, KnownNat n2, KnownNat m2, Elem a)+  => Matrix n m a+  -> Matrix n1 m1 a+  -> Matrix n2 m2 a+  -> a+relativeError x a b = Internal.performIO $+  unsafeWith x $ \x_vals x_rows x_cols ->+    unsafeWith a $ \a_vals a_rows a_cols ->+      unsafeWith b $ \b_vals b_rows b_cols ->+        alloca $ \pe -> do+          Internal.call $ Internal.relativeError pe+            x_vals x_rows x_cols+            a_vals a_rows a_cols+            b_vals b_rows b_cols+          fromC <$> peek pe++-- | The rank of the matrix.+rank :: (KnownNat n, KnownNat m, Elem a)+  => Decomposition+  -> Matrix n m a+  -> Int+rank d m = Internal.performIO $ alloca $ \pr -> do+  Internal.call $ unsafeWith m $ Internal.rank (con2CTag d) pr+  fromC <$> peek pr++-- | Return the matrix whose columns form a basis of the null-space of @A@.+kernel :: forall a n m. (Elem a, KnownNat n, KnownNat m)+  => Decomposition+  -> Matrix n m a+  -> Matrix n m a+kernel d m = Internal.performIO $+  alloca $ \pvals ->+  alloca $ \prows ->+  alloca $ \pcols ->+    unsafeWith m $ \vals1 rows1 cols1 -> do+      Internal.call $ Internal.kernel (con2CTag d)+        pvals prows pcols+        vals1 rows1 cols1+      vals <- peek pvals+      rs <- fromC <$> peek prows+      cs <- fromC <$> peek pcols+      fp <- FC.newForeignPtr vals $ Internal.free vals+      pure $ Matrix . Vec $ VS.unsafeFromForeignPtr0 fp (rs * cs)++-- | Return a matrix whose columns form a basis of the column-space of @A@.+image :: forall a n m. (Elem a, KnownNat n, KnownNat m)+  => Decomposition+  -> Matrix n m a+  -> Matrix n m a+image d m = Internal.performIO $+  alloca $ \pvals ->+  alloca $ \prows ->+  alloca $ \pcols ->+    unsafeWith m $ \vals1 rows1 cols1 -> do+      Internal.call $ Internal.image (con2CTag d)+        pvals prows pcols+        vals1 rows1 cols1+      vals <- peek pvals+      rs <- fromC <$> peek prows+      cs <- fromC <$> peek pcols+      fp <- FC.newForeignPtr vals $ Internal.free vals+      pure $ Matrix . Vec $ VS.unsafeFromForeignPtr0 fp (rs * cs)++{- |+[(coeffs, error) = linearRegression points] computes multiple linear regression @y = a1 x1 + a2 x2 + ... + an xn + b@ using 'ColPivHouseholderQR' decomposition++* point format is @[y, x1..xn]@++* coeffs format is @[b, a1..an]@++* error is calculated using 'relativeError'++@+import Data.Eigen.LA+main = print $ linearRegression (Row @5)+  [+    [-4.32, 3.02, 6.89],+    [-3.79, 2.01, 5.39],+    [-4.01, 2.41, 6.01],+    [-3.86, 2.09, 5.55],+    [-4.10, 2.58, 6.32]+  ]+@++ produces the following output++ @+ Just ([-2.3466569233817127,-0.2534897541434826,-0.1749653335680988],1.8905965120153139e-3)+ @++-}+linearRegression :: forall r. (KnownNat r)+  => Internal.Row r+  -- -> Internal.Col c+  -> [[Double]]+  -> Maybe ([Double], Double)+linearRegression _ points = do+  _a :: MatrixXd r 2 <- M.fromList $ List.map ((1:)  . tail) points+  _b :: MatrixXd r 1 <- M.fromList $ List.map ((:[]) . head) points+  let _x = solve ColPivHouseholderQR _a _b+  let e  = relativeError _x _a _b+  let coeffs = List.map head $ M.toList _x+  return (coeffs, e)

@@ -0,0 +1,449 @@+{-# LANGUAGE CPP #-}+{-# LANGUAGE DeriveFunctor #-}+{-# LANGUAGE FlexibleInstances #-}+{-# LANGUAGE ForeignFunctionInterface #-}+{-# LANGUAGE FunctionalDependencies #-}+{-# LANGUAGE MultiParamTypeClasses #-}+{-# LANGUAGE RecordWildCards #-}+{-# LANGUAGE GeneralizedNewtypeDeriving #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeFamilies #-}++{- |++This documentation is based on original Eigen page <http://eigen.tuxfamily.org/dox/group__TopicSparseSystems.html Solving Sparse Linear Systems>++Eigen currently provides a limited set of built-in MPL2 compatible solvers.+They are summarized in the following table:++@+Sparse solver       Solver kind             Matrix kind         Notes++ConjugateGradient   Classic iterative CG    SPD                 Recommended for large symmetric+                                                                problems (e.g., 3D Poisson eq.)+BiCGSTAB            Iterative stabilized    Square+                    bi-conjugate gradient+SparseLU            LU factorization        Square              Optimized for small and large problems+                                                                with irregular patterns+SparseQR            QR factorization        Any, rectangular    Recommended for least-square problems,+                                                                has a basic rank-revealing feature+@++All these solvers follow the same general concept. Here is a typical and general example:++@+let+    a :: SparseMatrixXd+    a = ... -- fill a++    b :: SparseMatrixXd+    b = ... -- fill b++    validate msg = info >>= (`when` fail msg) . (/= Success)++// solve Ax = b+runSolverT solver $ do+    compute a+    validate "decomposition failed"++    x <- solve b+    validate "solving failed"++    // solve for another right hand side+    x1 <- solve b1+@++In the case where multiple problems with the same sparsity pattern have to be solved, then the "compute" step can be decomposed as follow:++@+runSolverT solver $ do+    analyzePattern a1+    factorize a1+    x1 <- solve b1+    x2 <- solve b2++    factorize a2+    x1 <- solve b1+    x2 <- solve b2+@++Finally, each solver provides some specific features, such as determinant, access to the factors, controls of the iterations, and so on.++-}++module Eigen.Solver.SparseLA+  (+    -- * Sparse Solvers+    Solver +  , DirectSolver+  , IterativeSolver+  , OrderingMethod(..)+  , Preconditioner(..)+  , ConjugateGradient(..)+  , BiCGSTAB(..)+  , SparseLU(..)+  , SparseQR(..)+  , ComputationInfo(..)+  , SolverT(..)+  , runSolverT+  -- * The Compute step+  {- |+      In the `compute` function, the matrix is generally factorized: LLT for self-adjoint matrices, LDLT for general hermitian matrices,+      LU for non hermitian matrices and QR for rectangular matrices. These are the results of using direct solvers.+      For this class of solvers precisely, the compute step is further subdivided into `analyzePattern` and `factorize`.++      The goal of `analyzePattern` is to reorder the nonzero elements of the matrix, such that the factorization step creates less fill-in.+      This step exploits only the structure of the matrix. Hence, the results of this step can be used for other linear systems where the+      matrix has the same structure.++      In `factorize`, the factors of the coefficient matrix are computed. This step should be called each time the values of the matrix change.+      However, the structural pattern of the matrix should not change between multiple calls.++      For iterative solvers, the `compute` step is used to eventually setup a preconditioner.+      Remember that, basically, the goal of the preconditioner is to speedup the convergence of an iterative method by solving a modified linear+      system where the coefficient matrix has more clustered eigenvalues.+      For real problems, an iterative solver should always be used with a preconditioner.+  -}+  , analyzePattern+  , factorize+  , compute+  -- * The Solve step+  {- |+  The `solve` function computes the solution of the linear systems with one or many right hand sides.++  @+  x <- solve b+  @++  Here, @b@ can be a vector or a matrix where the columns form the different right hand sides.+  The `solve` function can be called several times as well, for instance when all the right hand sides are not available at once.++  @+  x1 <- solve b1+  -- Get the second right hand side b2+  x2 <- solve b2+  --  ...+  @+  -}+  , solve+  -- , solveWithGuess+  , info+  -- * Iterative Solvers+  , tolerance+  , setTolerance+  , maxIterations+  , setMaxIterations+  , error+  , iterations+  -- * SparseQR Solver+  , matrixR+  , matrixQ+  , rank+  , setPivotThreshold+  -- * SparseLU Solver+  , setSymmetric+  , matrixL+  , matrixU+  , determinant+  , absDeterminant+  , signDeterminant+  , logAbsDeterminant+  ) where++import Prelude hiding (error)+import Foreign.Ptr+import Foreign.ForeignPtr+import Foreign.Storable+import Foreign.C.String+import Foreign.Marshal.Alloc+import Control.Monad.IO.Class+import Control.Monad.Trans.Class+import Control.Monad.Trans.Reader+import qualified Foreign.Concurrent as FC+#if __GLASGOW_HASKELL__ >= 710+#else+import Control.Applicative+#endif+import qualified Eigen.Internal as I+import qualified Eigen.SparseMatrix as SM++{- | Ordering methods for sparse matrices. They are typically used to reduce the number of elements during the sparse matrix+    decomposition (@LLT@, @LU@, @QR@). Precisely, in a preprocessing step, a permutation matrix @P@ is computed using those ordering methods+    and applied to the columns of the matrix. Using for instance the sparse Cholesky decomposition, it is expected that the nonzeros+    elements in @LLT(A*P)@ will be much smaller than that in @LLT(A)@.+-}+data OrderingMethod+    -- | The column approximate minimum degree ordering The matrix should be in column-major and compressed format+    = COLAMDOrdering+    -- | The natural ordering (identity)+    | NaturalOrdering deriving (Show, Read)++data Preconditioner+    {- | A preconditioner based on the digonal entries++        It allows to approximately solve for A.x = b problems assuming A is a diagonal matrix.+        In other words, this preconditioner neglects all off diagonal entries and, in Eigen's language, solves for:+        @+        A.diagonal().asDiagonal() . x = b+        @+        This preconditioner is suitable for both selfadjoint and general problems.+        The diagonal entries are pre-inverted and stored into a dense vector.++        A variant that has yet to be implemented would attempt to preserve the norm of each column.+    -}+    = DiagonalPreconditioner+    -- | A naive preconditioner which approximates any matrix as the identity matrix+    | IdentityPreconditioner deriving (Show, Read)+++class I.Code s => Solver s where+-- | For direct methods, the solution is computed at the machine precision.+class Solver s => DirectSolver s where+-- | Sometimes, the solution need not be too accurate.+-- In this case, the iterative methods are more suitable and the desired accuracy can be set before the solve step using `setTolerance`.+class Solver s => IterativeSolver s where++{- | A conjugate gradient solver for sparse self-adjoint problems.++    This class allows to solve for @A.x = b@ sparse linear problems using a conjugate gradient algorithm. The sparse matrix @A@ must be selfadjoint.++    The maximal number of iterations and tolerance value can be controlled via the `setMaxIterations` and `setTolerance` methods.+    The defaults are the size of the problem for the maximal number of iterations and @epsilon@ for the tolerance+-}+newtype ConjugateGradient = ConjugateGradient Preconditioner deriving (Show, Read)+instance Solver ConjugateGradient+instance IterativeSolver ConjugateGradient+instance I.Code ConjugateGradient where+    code (ConjugateGradient DiagonalPreconditioner) = 0+    code (ConjugateGradient IdentityPreconditioner) = 1++{- | A bi conjugate gradient stabilized solver for sparse square problems.++    This class allows to solve for @A.x = b@ sparse linear problems using a bi conjugate gradient stabilized algorithm.+    The vectors @x@ and @b@ can be either dense or sparse.++    The maximal number of iterations and tolerance value can be controlled via the `setMaxIterations` and `setTolerance` methods.+    The defaults are the size of the problem for the maximal number of iterations and @epsilon@ for the tolerance+-}+newtype BiCGSTAB = BiCGSTAB Preconditioner deriving (Show, Read)+instance Solver BiCGSTAB+instance IterativeSolver BiCGSTAB+instance I.Code BiCGSTAB where+    code (BiCGSTAB DiagonalPreconditioner) = 2+    code (BiCGSTAB IdentityPreconditioner) = 3++{- | Sparse supernodal LU factorization for general matrices.++    This class implements the supernodal LU factorization for general matrices. It uses the main techniques from the sequential+    <http://crd-legacy.lbl.gov/~xiaoye/SuperLU/ SuperLU package>. It handles transparently real and complex arithmetics with+    single and double precision, depending on the scalar type of your input matrix. The code has been optimized to provide BLAS-3+    operations during supernode-panel updates. It benefits directly from the built-in high-performant Eigen BLAS routines.+    Moreover, when the size of a supernode is very small, the BLAS calls are avoided to enable a better optimization from the compiler.+    For best performance, you should compile it with NDEBUG flag to avoid the numerous bounds checking on vectors.++    An important parameter of this class is the ordering method. It is used to reorder the columns+    (and eventually the rows) of the matrix to reduce the number of new elements that are created during+    numerical factorization. The cheapest method available is COLAMD.+    See <http://eigen.tuxfamily.org/dox/group__OrderingMethods__Module.html OrderingMethods module> for the list of+    built-in and external ordering methods.+-}+newtype SparseLU = SparseLU OrderingMethod deriving (Show, Read)+instance Solver SparseLU+instance DirectSolver SparseLU+instance I.Code SparseLU where+    code (SparseLU NaturalOrdering) = 4+    code (SparseLU COLAMDOrdering) = 5++{- | Sparse left-looking rank-revealing QR factorization.++    This class implements a left-looking rank-revealing QR decomposition of sparse matrices. When a column has a norm less than a given+    tolerance it is implicitly permuted to the end. The QR factorization thus obtained is given by @A*P = Q*R@ where @R@ is upper triangular or trapezoidal.++    @P@ is the column permutation which is the product of the fill-reducing and the rank-revealing permutations.++    @Q@ is the orthogonal matrix represented as products of Householder reflectors.++    @R@ is the sparse triangular or trapezoidal matrix. The later occurs when @A@ is rank-deficient.+-}+newtype SparseQR = SparseQR OrderingMethod deriving (Show, Read)+instance Solver SparseQR+instance DirectSolver SparseQR+instance I.Code SparseQR where+    code (SparseQR NaturalOrdering) = 6+    code (SparseQR COLAMDOrdering) = 7+++data ComputationInfo+    -- | Computation was successful.+    = Success+    -- | The provided data did not satisfy the prerequisites.+    | NumericalIssue+    -- | Iterative procedure did not converge.+    | NoConvergence+    -- | The inputs are invalid, or the algorithm has been improperly called. When assertions are enabled, such errors trigger an error.+    | InvalidInput+    deriving (Eq, Enum, Show, Read)++newtype SolverT s a p c = SolverT (ReaderT (s, ForeignPtr (I.CSolver a)) p c)+  deriving (Functor, Applicative, Monad, MonadTrans)++runSolverT :: (Solver s, MonadIO p, I.Elem a) => s -> SolverT s a p c -> p c+runSolverT i (SolverT f) = do+  fs <- liftIO $ alloca $ \ps -> do+    I.call $ I.sparse_la_newSolver i ps+    s <- peek ps+    FC.newForeignPtr s (I.call $ I.sparse_la_freeSolver i s)+  runReaderT f (i,fs)++forSolverT :: (Solver s, MonadIO p, I.Elem a) => (s -> Ptr (I.CSolver a) -> Ptr (I.CSparseMatrix a) -> IO CString) -> SM.SparseMatrix n m a -> SolverT s a p ()+{-# INLINE forSolverT #-}+forSolverT f (SM.SparseMatrix fa) = SolverT $ ask >>= \(i,fs) -> liftIO $+  withForeignPtr fs $ \s ->+  withForeignPtr fa $ \a ->+    I.call $ f i s a++-- | Initializes the iterative solver for the sparsity pattern of the matrix @A@ for further solving @Ax=b@ problems.+analyzePattern :: (Solver s, MonadIO p, I.Elem a) => SM.SparseMatrix n m a -> SolverT s a p ()+analyzePattern sm = forSolverT I.sparse_la_analyzePattern sm++-- | Initializes the iterative solver with the numerical values of the matrix @A@ for further solving @Ax=b@ problems.+factorize :: (Solver s, MonadIO p, I.Elem a) => SM.SparseMatrix n m a -> SolverT s a p ()+factorize sm = forSolverT I.sparse_la_factorize sm++-- | Initializes the iterative solver with the matrix @A@ for further solving @Ax=b@ problems.+--+-- The `compute` method is equivalent to calling both `analyzePattern` and `factorize`.+compute :: (Solver s, MonadIO p, I.Elem a) => SM.SparseMatrix n m a -> SolverT s a p ()+compute sm = forSolverT I.sparse_la_compute sm++-- | An expression of the solution @x@ of @Ax=b@ using the current decomposition of @A@.+solve :: (Solver s, MonadIO p, I.Elem a) => SM.SparseMatrix n m a -> SolverT s a p (SM.SparseMatrix n m a)+solve (SM.SparseMatrix fb) = SolverT $ ask >>= \(i,fs) -> liftIO $+  withForeignPtr fs $ \s ->+  withForeignPtr fb $ \b ->+    alloca $ \px -> do+      I.call $ I.sparse_la_solve i s b px+      x <- peek px+      SM.SparseMatrix <$> FC.newForeignPtr x (I.call $ I.sparse_free x)++-- | The solution @x@ of @Ax=b@ using the current decomposition of @A@ and @x0@ as an initial solution.+--solveWithGuess :: (Solver s, MonadIO p, I.Elem a) => SM.SparseMatrix n m a -> SM.SparseMatrix n m a -> SolverT s a p (SM.SparseMatrix n m a)+--solveWithGuess (SM.SparseMatrix fb) (SM.SparseMatrix fx0) = SolverT $ ask >>= \(i,fs) -> liftIO $+--  withForeignPtr fs $ \s ->+--  withForeignPtr fb $ \b ->+--  withForeignPtr fx0 $ \x0 ->+--  alloca $ \px -> do+--    I.call $ I.sparse_la_solveWithGuess i s b x0 px+--    x <- peek px+--    SM.SparseMatrix <$> FC.newForeignPtr x (I.call $ I.sparse_free x)++-- |+-- * `Success` if the iterations converged or computation was succesful+-- * `NumericalIssue` if the factorization reports a numerical problem+-- * `NoConvergence` if the iterations are not converged+-- * `InvalidInput` if the input matrix is invalid+info :: (Solver s, MonadIO p, I.Elem a) => SolverT s a p ComputationInfo+info = _get_prop I.sparse_la_info >>= \x -> pure (toEnum x)++-- | The tolerance threshold used by the stopping criteria.+tolerance :: (IterativeSolver s, MonadIO p, I.Elem a) => SolverT s a p Double+tolerance = _get_prop I.sparse_la_tolerance++-- | Sets the tolerance threshold used by the stopping criteria.+--+--   This value is used as an upper bound to the relative residual error: @|Ax-b|/|b|@. The default value is the machine precision given by @epsilon@+setTolerance :: (IterativeSolver s, MonadIO p, I.Elem a) => Double -> SolverT s a p ()+setTolerance = _set_prop I.sparse_la_setTolerance++-- | The max number of iterations. It is either the value setted by setMaxIterations or, by default, twice the number of columns of the matrix.+maxIterations :: (IterativeSolver s, MonadIO p, I.Elem a) => SolverT s a p Int+maxIterations = _get_prop I.sparse_la_maxIterations++-- | Sets the max number of iterations. Default is twice the number of columns of the matrix.+setMaxIterations :: (IterativeSolver s, MonadIO p, I.Elem a) => Int -> SolverT s a p ()+setMaxIterations = _set_prop I.sparse_la_setMaxIterations++-- | The tolerance error reached during the last solve. It is a close approximation of the true relative residual error @|Ax-b|/|b|@.+error :: (IterativeSolver s, MonadIO p, I.Elem a) => SolverT s a p Double+error = _get_prop I.sparse_la_error++-- | The number of iterations performed during the last solve+iterations :: (IterativeSolver s, MonadIO p, I.Elem a) => SolverT s a p Int+iterations = _get_prop I.sparse_la_iterations++-- | Returns the @b@ sparse upper triangular matrix @R@ of the QR factorization.+matrixR :: (MonadIO p, I.Elem a) => SolverT SparseQR a p (SM.SparseMatrix n m a)+matrixR = _get_matrix I.sparse_la_matrixR++-- | Returns the matrix @Q@ as products of sparse Householder reflectors.+matrixQ :: (MonadIO p, I.Elem a) => SolverT SparseQR a p (SM.SparseMatrix n m a)+matrixQ = _get_matrix I.sparse_la_matrixQ++-- | Sets the threshold that is used to determine linearly dependent columns during the factorization.+--+-- In practice, if during the factorization the norm of the column that has to be eliminated is below+-- this threshold, then the entire column is treated as zero, and it is moved at the end.+setPivotThreshold :: (MonadIO p, I.Elem a) => Double -> SolverT SparseQR a p ()+setPivotThreshold = _set_prop I.sparse_la_setPivotThreshold++-- | Returns the number of non linearly dependent columns as determined by the pivoting threshold.+rank :: (MonadIO p, I.Elem a) => SolverT SparseQR a p Int+rank = _get_prop I.sparse_la_rank++-- | Indicate that the pattern of the input matrix is symmetric+setSymmetric :: (MonadIO p, I.Elem a) => Bool -> SolverT SparseLU a p ()+setSymmetric = _set_prop I.sparse_la_setSymmetric . fromEnum++-- | Returns the matrix @L@+matrixL :: (MonadIO p, I.Elem a) => SolverT SparseLU a p (SM.SparseMatrix n m a)+matrixL = _get_matrix I.sparse_la_matrixL++-- | Returns the matrix @U@+matrixU :: (MonadIO p, I.Elem a) => SolverT SparseLU a p (SM.SparseMatrix n m a)+matrixU = _get_matrix I.sparse_la_matrixU++-- | The determinant of the matrix.+determinant :: (MonadIO p, I.Elem a) => SolverT SparseLU a p a+determinant = _get_prop I.sparse_la_determinant++-- | The natural log of the absolute value of the determinant of the matrix of which this is the QR decomposition+--+-- This method is useful to work around the risk of overflow/underflow that's inherent to the determinant computation.+logAbsDeterminant :: (MonadIO p, I.Elem a) => SolverT SparseLU a p a+logAbsDeterminant = _get_prop I.sparse_la_logAbsDeterminant++-- | The absolute value of the determinant of the matrix of which *this is the QR decomposition.+--+-- A determinant can be very big or small, so for matrices of large enough dimension, there is a risk of overflow/underflow.+-- One way to work around that is to use `logAbsDeterminant` instead.+absDeterminant :: (MonadIO p, I.Elem a) => SolverT SparseLU a p a+absDeterminant = _get_prop I.sparse_la_absDeterminant++-- | A number representing the sign of the determinant+signDeterminant :: (MonadIO p, I.Elem a) => SolverT SparseLU a p a+signDeterminant = _get_prop I.sparse_la_signDeterminant++_get_prop :: (I.Cast d, Solver s, MonadIO p, Storable c, c ~ I.C d)+  => (s -> I.CSolverPtr a -> Ptr c -> IO CString)+  -> SolverT s a p d+_get_prop f = SolverT $ ask >>= \(i, fs) -> liftIO $+  withForeignPtr fs $ \s -> alloca $ \px -> do+    I.call (f i s px)+    I.fromC <$> peek px++_get_matrix :: (Solver s, MonadIO p, I.Elem a)+  => (s -> I.CSolverPtr a -> Ptr (I.CSparseMatrixPtr a) -> IO CString)+  -> SolverT s a p (SM.SparseMatrix n m a)+_get_matrix f = SolverT $ ask >>= \(i,fs) -> liftIO $+  withForeignPtr fs $ \s -> alloca $ \px -> do+    I.call (f i s px)+    x <- peek px+    SM.SparseMatrix <$> FC.newForeignPtr x (I.call (I.sparse_free x))++_set_prop :: (I.Cast c, d ~ I.C c, Solver s, MonadIO p, Storable c)+  => (s -> I.CSolverPtr a -> d -> IO CString)+  -> c+  -> SolverT s a p ()+_set_prop f x = SolverT $ ask >>= \(i,fs) -> liftIO $+  withForeignPtr fs $ \s -> I.call (f i s (I.toC x))+

@@ -0,0 +1,504 @@+{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE DataKinds           #-}+{-# LANGUAGE GADTs               #-}+{-# LANGUAGE KindSignatures      #-}+{-# LANGUAGE RecordWildCards     #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications    #-}+{-# LANGUAGE TypeInType          #-}+{-# LANGUAGE TypeOperators       #-}++module Eigen.SparseMatrix+  ( -- * Types+    SparseMatrix(..)+  , SparseMatrixXf+  , SparseMatrixXd+  , SparseMatrixXcf+  , SparseMatrixXcd++    -- * Matrix internal data+  , elems+  , values+  , innerIndices+  , outerStarts+  , innerNNZs++    -- * Accessors+  , cols+  , rows+  , coeff+  , (!)+  , getRow+  , getCol++    -- * Matrix conversions+  , fromList+  , toList+  , fromVector+  , toVector+  , fromDenseList+  , toDenseList+  , fromMatrix+  , toMatrix++    -- * Matrix properties+  , norm+  , squaredNorm+  , blueNorm+  , block+  , nonZeros+  , innerSize+  , outerSize++    -- * Basic matrix algebra+  , add+  , sub+  , mul++    -- * Matrix tranformations+  , pruned+  , scale+  , transpose+  , adjoint+  , map+  , imap++    -- * Matrix representation+  , compress+  , uncompress+  , compressed++    -- * Matrix serialisation+  , encode+  , decode++    -- * Mutable matrices+  , thaw+  , freeze+  , unsafeThaw+  , unsafeFreeze+  ) where++import Control.Monad (when, guard)+import Control.Monad.Primitive (PrimMonad(..), unsafePrimToPrim)+import qualified Prelude+import Prelude hiding (read, map)+import Data.Binary (Binary(..))+import qualified Data.Binary as Binary+import qualified Data.ByteString.Lazy as BSL+import Data.Complex (Complex)+import Data.Kind (Type)+import qualified Data.Vector.Storable as VS+import qualified Data.Vector.Storable.Mutable as VSM+import Foreign.C.String (CString)+import Foreign.C.Types (CInt)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr, touchForeignPtr)+import Foreign.Marshal.Alloc (alloca)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable(..))+import qualified Foreign.Concurrent as FC+import GHC.TypeLits (Nat, KnownNat, type (<=))+import Eigen.Internal+  ( Elem+  , Cast(..)+  , CSparseMatrix+  , CSparseMatrixPtr+  , natToInt+  , Row(..)+  , Col(..)+  , CTriplet(..)+  )+import qualified Data.List as List+import qualified Eigen.Internal as Internal+import qualified Eigen.Matrix as M+import qualified Eigen.Matrix.Mutable as MM+import qualified Eigen.SparseMatrix.Mutable as SMM++{-| A versatible sparse matrix representation.+SparseMatrix is the main sparse matrix representation of Eigen's sparse module.+It offers high performance and low memory usage.+It implements a more versatile variant of the widely-used Compressed Column (or Row) Storage scheme.+It consists of four compact arrays:+* `values`: stores the coefficient values of the non-zeros.+* `innerIndices`: stores the row (resp. column) indices of the non-zeros.+* `outerStarts`: stores for each column (resp. row) the index of the first non-zero in the previous two arrays.+* `innerNNZs`: stores the number of non-zeros of each column (resp. row). The word inner refers to an inner vector that is a column for a column-major matrix, or a row for a row-major matrix. The word outer refers to the other direction.+This storage scheme is better explained on an example. The following matrix+@+0   3   0   0   0+22  0   0   0   17+7   5   0   1   0+0   0   0   0   0+0   0   14  0   8+@+and one of its possible sparse, __column major__ representation:+@+values:         22  7   _   3   5   14  _   _   1   _   17  8+innerIndices:   1   2   _   0   2   4   _   _   2   _   1   4+outerStarts:    0   3   5   8   10  12+innerNNZs:      2   2   1   1   2+@+Currently the elements of a given inner vector are guaranteed to be always sorted by increasing inner indices.+The "\_" indicates available free space to quickly insert new elements. Assuming no reallocation is needed,+the insertion of a random element is therefore in @O(nnz_j)@ where @nnz_j@ is the number of nonzeros of the+respective inner vector. On the other hand, inserting elements with increasing inner indices in a given inner+vector is much more efficient since this only requires to increase the respective `innerNNZs` entry that is a @O(1)@ operation.+The case where no empty space is available is a special case, and is refered as the compressed mode.+It corresponds to the widely used Compressed Column (or Row) Storage schemes (CCS or CRS).+Any `SparseMatrix` can be turned to this form by calling the `compress` function.+In this case, one can remark that the `innerNNZs` array is redundant with `outerStarts` because we the equality:+@InnerNNZs[j] = OuterStarts[j+1]-OuterStarts[j]@. Therefore, in practice a call to `compress` frees this buffer.+The results of Eigen's operations always produces compressed sparse matrices.+On the other hand, the insertion of a new element into a `SparseMatrix` converts this later to the uncompressed mode.+For more infomration please see Eigen <http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html documentation page>.+-}+newtype SparseMatrix :: Nat -> Nat -> Type -> Type where+  SparseMatrix :: ForeignPtr (CSparseMatrix a) -> SparseMatrix n m a++instance forall n m a. (Elem a, Show a, KnownNat n, KnownNat m) => Show (SparseMatrix n m a) where+  show m = concat+    [ "SparseMatrix"+    , show (rows_ m)+    , "x"+    , show (cols_ m)+    , "\n"+    , List.intercalate "\n" $ Prelude.map (List.intercalate "\t" . Prelude.map show) $ toDenseList m+    , "\n"+    ]++instance forall n m a. (Elem a, KnownNat n, KnownNat m) => Binary (SparseMatrix n m a) where+  put mat = do+    put $ Internal.magicCode (undefined :: C a)+    put $ natToInt @n+    put $ natToInt @m+    put $ toVector mat+  +  get = do+    get >>= (`when` fail "wrong matrix type") . (/= Internal.magicCode (undefined :: C a))+    fromVector <$> get++-- | Encode the sparse matrix as a lazy bytestring+encode :: (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> BSL.ByteString+encode = Binary.encode++-- | Decode the sparse matrix from a lazy bytestring+decode :: (Elem a, KnownNat n, KnownNat m) => BSL.ByteString -> SparseMatrix n m a+decode = Binary.decode++-- | Alias for single precision sparse matrix+type SparseMatrixXf  n m = SparseMatrix n m Float+-- | Alias for double precision sparse matrix+type SparseMatrixXd  n m = SparseMatrix n m Double+-- | Alias for single previsiom sparse matrix of complex numbers+type SparseMatrixXcf n m = SparseMatrix n m (Complex Float)+-- | Alias for double prevision sparse matrix of complex numbers+type SparseMatrixXcd n m = SparseMatrix n m (Complex Double)++-- | Get the coefficient values of the non-zeros.+values :: Elem a => SparseMatrix n m a -> VS.Vector a+values = VS.map fromC . _getvec Internal.sparse_values++-- | Get the row indices of the non-zeros.+innerIndices :: Elem a => SparseMatrix n m a -> VS.Vector Int+innerIndices = VS.map fromC . _getvec Internal.sparse_innerIndices++-- | Gets for each column the index of the first non-zero in the previous two arrays.+outerStarts :: Elem a => SparseMatrix n m a -> VS.Vector Int+outerStarts = VS.map fromC . _getvec Internal.sparse_outerStarts++-- | Gets the number of non-zeros of each column.+-- The word inner refers to an inner vector that is a column for a column-major matrix, or a row for a row-major matrix.+-- The word outer refers to the other direction+innerNNZs :: Elem a => SparseMatrix n m a -> Maybe (VS.Vector Int)+innerNNZs m+  | compressed m = Nothing+  | otherwise    = Just $ VS.map fromC $ _getvec Internal.sparse_innerNNZs m++-- | Number of rows in the sparse matrix+rows :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Row n+rows _ = Row @n++-- | Number of colums in the sparse matrix+cols :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Col m+cols _ = Col @m++-- | Number of rows in the sparse matrix+rows_ :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int+rows_ _ = natToInt @n++-- | Number of colums in the sparse matrix+cols_ :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int+cols_ _ = natToInt @m+++-- | Sparse matrix coefficient at the given row and column+coeff :: forall n m r c a. (Elem a, KnownNat n, KnownNat m, KnownNat r, KnownNat c, r <= n, c <= m)+  => Row r -> Col c -> SparseMatrix n m a -> a+coeff _ _ (SparseMatrix fp) =+  let !c_row = toC $! natToInt @r+      !c_col = toC $! natToInt @c+  in Internal.performIO $ withForeignPtr fp $ \p -> alloca $ \pq -> do+       Internal.call $ Internal.sparse_coeff p c_row c_col pq+       fromC <$> peek pq++-- | Matrix coefficient at the given row and column+(!) :: forall n m r c a. (Elem a, KnownNat n, KnownNat m, KnownNat r, KnownNat c, r <= n, c <= m)+  => SparseMatrix n m a -> (Row r, Col c) -> a+(!) m (row,col) = coeff row col m++{-| For vectors, the l2 norm, and for matrices the Frobenius norm.+    In both cases, it consists in the square root of the sum of the square of all the matrix entries.+    For vectors, this is also equals to the square root of the dot product of this with itself.+-}+norm :: Elem a => SparseMatrix n m a -> a+norm = _unop Internal.sparse_norm (pure . fromC)++-- | For vectors, the squared l2 norm, and for matrices the Frobenius norm. In both cases, it consists in the sum of the square of all the matrix entries. For vectors, this is also equals to the dot product of this with itself.+squaredNorm :: Elem a => SparseMatrix n m a -> a+squaredNorm = _unop Internal.sparse_squaredNorm (pure . fromC)++-- | The l2 norm of the matrix using the Blue's algorithm. A Portable Fortran Program to Find the Euclidean Norm of a Vector, ACM TOMS, Vol 4, Issue 1, 1978.+blueNorm :: Elem a => SparseMatrix n m a -> a+blueNorm = _unop Internal.sparse_blueNorm (pure . fromC)++-- | Extract a rectangular block from the sparse matrix, given a startRow, startCol, blockRows, blockCols+block :: forall sr sc br bc n m a.+     (Elem a, KnownNat sr, KnownNat sc, KnownNat br, KnownNat bc, KnownNat n, KnownNat m)+  => (sr <= n, sc <= m, br <= n, bc <= m)+  => Row sr+  -> Col sc+  -> Row br+  -> Col bc+  -> SparseMatrix n m a+  -> SparseMatrix br bc a+block _ _ _ _ =+  let !c_startRow = toC $! natToInt @sr+      !c_startCol = toC $! natToInt @sc+      !c_rows     = toC $! natToInt @br+      !c_cols     = toC $! natToInt @bc+  in _unop (\p pq -> Internal.sparse_block p c_startRow c_startCol c_rows c_cols pq) _mk++-- | Number of non-zeros elements in the sparse matrix+nonZeros :: Elem a => SparseMatrix n m a -> Int+nonZeros = _unop Internal.sparse_nonZeros (pure . fromC)++-- | Number of elements in the sparse matrix, including zeros+elems :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int+elems _ = (natToInt @n) * (natToInt @m)++-- | The matrix in the compressed format+compress :: Elem a => SparseMatrix n m a -> SparseMatrix n m a+compress = _unop Internal.sparse_makeCompressed _mk++-- | The matrix in the uncompressed format+uncompress :: Elem a => SparseMatrix n m a -> SparseMatrix n m a+uncompress = _unop Internal.sparse_uncompress _mk++-- | Is the matrix compressed?+compressed :: Elem a => SparseMatrix n m a -> Bool+compressed = _unop Internal.sparse_isCompressed (pure . (/=0))++-- | Minor dimension with respect to the storage order+innerSize :: Elem a => SparseMatrix n m a -> Int+innerSize = _unop Internal.sparse_innerSize (pure . fromC)++-- | Major dimension with respect to the storage order+outerSize :: Elem a => SparseMatrix n m a -> Int+outerSize = _unop Internal.sparse_outerSize (pure . fromC)++-- | Suppresses all nonzeros which are much smaller than the reference under the tolerance @epsilon@+pruned :: Elem a => a -> SparseMatrix n m a -> SparseMatrix n m a+pruned r = _unop (\p pq -> alloca $ \pr -> poke pr (toC r) >> Internal.sparse_prunedRef p pr pq) _mk++-- | Multiply matrix on a given scalar+scale :: Elem a => a -> SparseMatrix n m a -> SparseMatrix n m a+scale x = _unop (\p pq -> alloca $ \px -> poke px (toC x) >> Internal.sparse_scale p px pq) _mk++-- | Transpose of the sparse matrix+transpose :: Elem a => SparseMatrix n m a -> SparseMatrix m n a+transpose = _unop Internal.sparse_transpose _mk++-- | Adjoint of the sparse matrix+adjoint :: Elem a => SparseMatrix n m a -> SparseMatrix m n a+adjoint = _unop Internal.sparse_adjoint _mk++-- | Add two sparse matrices by adding the corresponding entries together.+add :: Elem a => SparseMatrix n m a -> SparseMatrix n m a -> SparseMatrix n m a+add = _binop Internal.sparse_add _mk++-- | Subtract two sparse matrices by subtracting the corresponding entries together.+sub :: Elem a => SparseMatrix n m a -> SparseMatrix n m a -> SparseMatrix n m a+sub = _binop Internal.sparse_sub _mk++-- | Matrix multiplication.+mul :: Elem a => SparseMatrix p q a -> SparseMatrix q r a -> SparseMatrix p r a+mul = _binop Internal.sparse_mul _mk++-- | Map a function over the 'SparseMatrix'.+map :: (Elem a, Elem b, KnownNat n, KnownNat m) => (a -> b) -> SparseMatrix n m a -> SparseMatrix n m b+map f m = fromVector . VS.map g . toVector $ m where+  g (CTriplet r c v) = CTriplet r c $ (toC . f . fromC) v++-- | Map an indexed function over the 'SparseMatrix'.+imap :: (Elem a, Elem b, KnownNat n, KnownNat m) => (Int -> Int -> a -> b) -> SparseMatrix n m a -> SparseMatrix n m b+imap f m = fromVector . VS.map g . toVector $ m where+  g (CTriplet r c v) =+    let !_r = fromC r+        !_c = fromC c+        !_v = fromC v+    in CTriplet r c $ toC $ f _r _c _v++-- | Construct asparse matrix of the given size from the storable vector of triplets (row, col, val)+fromVector :: forall n m a. (Elem a, KnownNat n, KnownNat m)+  => VS.Vector (CTriplet a)+  -> SparseMatrix n m a+fromVector tris =+  let !c_rs = toC $! natToInt @n+      !c_cs = toC $! natToInt @m+      !len  = toC $! VS.length tris+  in Internal.performIO $ VS.unsafeWith tris $ \p ->+       alloca $ \pq -> do+         Internal.call $ Internal.sparse_fromList c_rs c_cs p len pq+         peek pq >>= _mk++-- | Convert a sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included.+toVector :: Elem a => SparseMatrix n m a -> VS.Vector (CTriplet a)+toVector m@(SparseMatrix fp) = Internal.performIO $ do+  let !size = nonZeros m+  tris <- VSM.new size+  withForeignPtr fp $ \p ->+    VSM.unsafeWith tris $ \q ->+      Internal.call $ Internal.sparse_toList p q (toC size)+  VS.unsafeFreeze tris++-- | Convert a sparse matrix to the list of triplets (row, col, val). Compressed elements will not be included.+toList :: Elem a => SparseMatrix n m a -> [(Int, Int, a)]+toList = Prelude.map fromC . VS.toList . toVector++-- | Construct a sparse matrix from a list of triples (row, val, col)+--+fromList :: (Elem a, KnownNat n, KnownNat m) => [(Int, Int, a)] -> SparseMatrix n m a+fromList = fromVector . VS.fromList . fmap toC++-- | Convert a sparse matrix to a (n X m) dense list of values.+toDenseList :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> [[a]]+toDenseList mat = [[_unsafeCoeff row col mat | col <- [0 .. _unsafeCols mat - 1]] | row <- [0 .. _unsafeRows mat - 1]]++-- | Construct a sparsematrix from a two-dimensional list of values.+--   If the dimensions of the list do not match that of the list, 'Nothing' is returned.+--   Zero values will be compressed.+fromDenseList :: forall n m a. (Elem a, Eq a, KnownNat n, KnownNat m) => [[a]] -> Maybe (SparseMatrix n m a)+fromDenseList list =+  let _rows = List.length list+      _cols = List.foldl' max 0 $ List.map length list+  in if ((_rows /= (natToInt @n)) || (_cols /= (natToInt @m)))+    then Nothing+    else Just $ fromList $ do+      (row, vals) <- zip [0..] list+      (col, val) <- zip [0..] vals+      guard $ val /= 0+      return (row, col, val)+      +  +--  fromList $ do+-- (row, vals) <- zip [0..] list+-- (col, val) <- zip [0..] vals+-- guard $ val /= 0+-- return (row, col, val)+-- where+--   rows = List.length list+--   cols = List.foldl' max 0 $ List.map length list++-- | Construct a dense matrix from a sparse matrix+toMatrix :: (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> M.Matrix n m a+toMatrix (SparseMatrix fp) = Internal.performIO $ do+  m0 :: MM.IOMatrix n m a <- MM.new+  MM.unsafeWith m0 $ \_vals _rows _cols ->+    withForeignPtr fp $ \pm1 ->+      Internal.call $ Internal.sparse_toMatrix pm1 _vals _rows _cols+  M.unsafeFreeze m0++-- | Construct a sparse matrix from a dense matrix. zero-elements will be compressed.+fromMatrix :: (Elem a, KnownNat n, KnownNat m) => M.Matrix n m a -> SparseMatrix n m a+fromMatrix m1 = Internal.performIO $ alloca $ \pm0 ->+  M.unsafeWith m1 $ \_vals _rows _cols -> do+    Internal.call $ Internal.sparse_fromMatrix _vals _rows _cols pm0+    peek pm0 >>= _mk++-- | Yield an immutable copy of the mutable matrix+freeze :: (Elem a, PrimMonad p) => SMM.MSparseMatrix n m (PrimState p) a -> p (SparseMatrix n m a)+freeze (SMM.MSparseMatrix fp) = SparseMatrix <$> _clone fp++-- | Yield a mutable copy of the immutable matrix.+thaw :: (Elem a, PrimMonad p) => SparseMatrix n m a -> p (SMM.MSparseMatrix n m (PrimState p) a)+thaw (SparseMatrix fp) = SMM.MSparseMatrix <$> _clone fp++-- | Unsafely convert a mutable matrix to an immutable one without copying. The mutable matrix may not be used after this operation.+unsafeFreeze :: (Elem a, PrimMonad p) => SMM.MSparseMatrix n m (PrimState p) a -> p (SparseMatrix n m a) +unsafeFreeze (SMM.MSparseMatrix fp) = return $! SparseMatrix fp++-- | Unsafely convert an immutable matrix to a mutable one without copying. The immutable matrix may not be used after this operation.+unsafeThaw :: (Elem a, PrimMonad p) => SparseMatrix n m a -> p (SMM.MSparseMatrix n m (PrimState p) a) +unsafeThaw (SparseMatrix fp) = return $! SMM.MSparseMatrix fp++-- | Return a single row of the sparse matrix.+getRow :: forall n m r a. (Elem a, KnownNat n, KnownNat m, KnownNat r, r <= n, 1 <= n) => Row r -> SparseMatrix n m a -> SparseMatrix 1 m a+getRow row mat = block row (Col @0) (Row @1) (Col @m) mat++-- | Return a single column of the sparse matrix.+getCol :: forall n m c a. (Elem a, KnownNat n, KnownNat m, KnownNat c, c <= m, 1 <= m) => Col c -> SparseMatrix n m a -> SparseMatrix n 1 a+getCol col mat = block (Row @0) col (Row @n) (Col @1) mat++_unop :: Storable b => (CSparseMatrixPtr a -> Ptr b -> IO CString) -> (b -> IO c) -> SparseMatrix n m a -> c+_unop f g (SparseMatrix fp) = Internal.performIO $+  withForeignPtr fp $ \p ->+    alloca $ \pq -> do+      Internal.call (f p pq)+      peek pq >>= g++_binop :: Storable b => (CSparseMatrixPtr a -> CSparseMatrixPtr a -> Ptr b -> IO CString) -> (b -> IO c) -> SparseMatrix n m a -> SparseMatrix n1 m1 a -> c+_binop f g (SparseMatrix fp1) (SparseMatrix fp2) = Internal.performIO $+  withForeignPtr fp1 $ \p1 ->+    withForeignPtr fp2 $ \p2 ->+      alloca $ \pq -> do+        Internal.call (f p1 p2 pq)+        peek pq >>= g++_getvec :: (Elem a, Storable b) => (Ptr (CSparseMatrix a) -> Ptr CInt -> Ptr (Ptr b) -> IO CString) -> SparseMatrix n m a -> VS.Vector b+_getvec f (SparseMatrix fm) = Internal.performIO $+  withForeignPtr fm $ \m ->+    alloca $ \ps ->+      alloca $ \pq -> do+        Internal.call $ f m ps pq+        s <- fromIntegral <$> peek ps+        q <- peek pq+        fr <- FC.newForeignPtr q $ touchForeignPtr fm+        pure $! VS.unsafeFromForeignPtr0 fr s++_clone :: (PrimMonad p, Elem a) => ForeignPtr (CSparseMatrix a) -> p (ForeignPtr (CSparseMatrix a))+_clone fp = unsafePrimToPrim $ withForeignPtr fp $ \p -> alloca $ \pq -> do+  Internal.call $ Internal.sparse_clone p pq+  q <- peek pq+  FC.newForeignPtr q $ Internal.call $ Internal.sparse_free q++_mk :: Elem a => Ptr (CSparseMatrix a) -> IO (SparseMatrix n m a)+_mk p = SparseMatrix <$> FC.newForeignPtr p (Internal.call $ Internal.sparse_free p)++-- | Number of rows in the sparse matrix+_unsafeRows :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int+{-# INLINE _unsafeRows #-}+_unsafeRows _ = natToInt @n++-- | Number of colums in the sparse matrix+_unsafeCols :: forall n m a. (Elem a, KnownNat n, KnownNat m) => SparseMatrix n m a -> Int+{-# INLINE _unsafeCols #-}+_unsafeCols _ = natToInt @m++_unsafeCoeff :: (Elem a, KnownNat n) => Int -> Int -> SparseMatrix n m a -> a+{-# INLINE _unsafeCoeff #-}+_unsafeCoeff !row !col (SparseMatrix fp) =+  let !c_row = toC row+      !c_col = toC col+  in Internal.performIO $ withForeignPtr fp $ \p -> alloca $ \pq -> do+       Internal.call $ Internal.sparse_coeff p c_row c_col pq+       fromC <$> peek pq

@@ -0,0 +1,163 @@+{-# LANGUAGE BangPatterns        #-}+{-# LANGUAGE DataKinds           #-}+{-# LANGUAGE GADTs               #-}+{-# LANGUAGE KindSignatures      #-}+{-# LANGUAGE RecordWildCards     #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications    #-}+{-# LANGUAGE TypeInType          #-}+{-# LANGUAGE TypeOperators       #-}++module Eigen.SparseMatrix.Mutable+  ( +    -- * Mutable SparseMatrix +    MSparseMatrix(..)+  , IOSparseMatrix+  , STSparseMatrix+  , new+  , reserve++    -- * SparseMatrix properties+  , rows+  , cols+  , innerSize+  , outerSize+  , nonZeros++    -- * SparseMatrix compression+  , compressed+  , compress+  , uncompress++    -- * Accessing SparseMatrix data+  , read+  , write+  , setZero+  , setIdentity+  ) where++import Prelude hiding (read)+import Control.Monad.Primitive (PrimMonad(..), unsafePrimToPrim)+import Data.Kind (Type)+import Foreign.C.String (CString)+import Foreign.ForeignPtr (ForeignPtr, withForeignPtr)+import Foreign.Marshal.Alloc (alloca)+import Foreign.Ptr (Ptr)+import Foreign.Storable (Storable(..))+import qualified Foreign.Concurrent as FC+import GHC.Exts (RealWorld)+import GHC.TypeLits (Nat, KnownNat, type (<=))+import Eigen.Internal+  ( Elem+  , Cast(..)+  , CSparseMatrix+  , CSparseMatrixPtr+  , natToInt+  , Row+  , Col+  )+import qualified Eigen.Internal as Internal++-- | Mutable sparse matrix. See 'Eigen.SparseMatrix.SparseMatrix' for+--   details about matrix layout.+newtype MSparseMatrix :: Nat -> Nat -> Type -> Type -> Type where+  MSparseMatrix :: (ForeignPtr (CSparseMatrix a)) -> MSparseMatrix n m s a++-- | A sparse matrix where the state token is specialised to 'ReadWorld'.+type IOSparseMatrix n m   a = MSparseMatrix n m RealWorld a+-- | This type does not differ from 'MSparseMatrix', but might be desirable+--   for readability.+type STSparseMatrix n m s a = MSparseMatrix n m s a++-- | Create a new sparse matrix with the given size @rows x cols@.+new :: forall m n p a. (Elem a, KnownNat n, KnownNat m, PrimMonad p) => p (MSparseMatrix n m (PrimState p) a)+new = unsafePrimToPrim $ alloca $ \pm -> do+  let !c_rs = toC $! natToInt @n+      !c_cs = toC $! natToInt @m+  Internal.call $ Internal.sparse_new c_rs c_cs pm+  m <- peek pm+  fm <- FC.newForeignPtr m $ Internal.call $ Internal.sparse_free m+  pure $! MSparseMatrix fm++-- | Returns the number of rows of the matrix.+rows :: forall n m s a. (Elem a, KnownNat n, KnownNat m) => MSparseMatrix n m s a -> Int+rows _ = natToInt @n++-- | Returns the number of columns of the matrix.+cols :: forall n m s a. (Elem a, KnownNat n, KnownNat m) => MSparseMatrix n m s a -> Int+cols _ = natToInt @m++-- | Returns the number of rows (resp. columns) of the matrix if the storage order is column majour (resp. row majour)+innerSize :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p Int+innerSize = _prop Internal.sparse_innerSize (pure . fromC)++-- | Returns the number of columns (resp. rows) of the matrix if the storage order is column majour (resp. row majour)+outerSize :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p Int+outerSize = _prop Internal.sparse_outerSize (pure . fromC)++-- | Returns whether or not the matrix is in compressed form.+compressed :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p Bool+compressed = _prop Internal.sparse_isCompressed (pure . (== 1))++-- | Turns the matrix into compressed format.+compress :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p ()+compress = _inplace Internal.sparse_compressInplace++-- | Decompresses the matrix.+uncompress :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p ()+uncompress = _inplace Internal.sparse_uncompressInplace++-- | Read the value of the matrix at position i,j. This function returns @Scalar(0)@ if the element is an explicit 0.+read :: forall n m r c p a. (Elem a, PrimMonad p, KnownNat n, KnownNat m, KnownNat r, KnownNat c, r <= n, c <= m)+  => Row r+  -> Col c+  -> MSparseMatrix n m (PrimState p) a+  -> p a+read _ _ (MSparseMatrix fm) =+  let !c_r = toC $! natToInt @r+      !c_c = toC $! natToInt @c+  in unsafePrimToPrim $ withForeignPtr fm $ \m -> alloca $ \px -> do+       Internal.call $ Internal.sparse_coeff m c_r c_c px+       fromC <$> peek px++{- | Writes the value of the matrix at position @i@, @j@.+     This function turns the matrix into a non compressed form if that was not the case.++     This is a @O(log(nnz_j))@ operation (binary search) plus the cost of element insertion if the element does not already exist.++     Cost of element insertion is sorted insertion in O(1) if the elements of each inner vector are inserted in increasing inner index order, and in @O(nnz_j)@ for a random insertion.+-}+write :: forall n m r c p a. (Elem a, PrimMonad p, KnownNat n, KnownNat m, KnownNat r, KnownNat c, r <= n, c <= m)+  => MSparseMatrix n m (PrimState p) a -> Row r -> Col c -> a -> p ()+write (MSparseMatrix fm) _ _ x =+  let !c_r = toC $! natToInt @r+      !c_c = toC $! natToInt @c+  in unsafePrimToPrim $ withForeignPtr fm $ \m -> alloca $ \px -> do+       Internal.call $ Internal.sparse_coeffRef m c_r c_c px+       peek px >>= (`poke` toC x)++-- | Sets the matrix to the identity matrix.+setIdentity :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p ()+setIdentity = _inplace Internal.sparse_setIdentity++-- | Remove all non zeros, but keep allocated memory.+setZero :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p ()+setZero = _inplace Internal.sparse_setZero++-- | Preallocates for non zeros. The matrix must be in compressed mode.+reserve :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> Int -> p ()+reserve m s = _inplace (\p -> Internal.sparse_reserve p (toC s)) m++-- | Returns the number of nonzero coefficients.+nonZeros :: (Elem a, PrimMonad p) => MSparseMatrix n m (PrimState p) a -> p Int+nonZeros = _prop Internal.sparse_nonZeros (pure . fromC)++_inplace :: (Elem a, PrimMonad p) => (Ptr (CSparseMatrix a) -> IO CString) -> MSparseMatrix n m (PrimState p) a -> p ()+_inplace f (MSparseMatrix fm) = unsafePrimToPrim $ withForeignPtr fm $ \m -> Internal.call $ f m++_prop :: (Storable b, PrimMonad p) => (CSparseMatrixPtr a -> Ptr b -> IO CString) -> (b -> IO c) -> MSparseMatrix n m (PrimState p) a -> p c+_prop f g (MSparseMatrix fp) = unsafePrimToPrim $ +  withForeignPtr fp $ \p ->+    alloca $ \pq -> do+      Internal.call (f p pq)+      peek pq >>= g

@@ -1,16 +1,22 @@-import Data.Eigen.Matrix-import Data.Eigen.LA+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-} +import Eigen.Matrix+import Eigen.Solver.LA++a :: Maybe (MatrixXf 3 3)+a = fromList [[1,2,5],[2,1,4],[3,0,3]]++main :: IO () main = do-    let a = fromList [[1,2,5],[2,1,4],[3,0,3]] :: MatrixXf     putStrLn "Here is the matrix A:"     print a      putStrLn "The rank of A is:"-    print $ rank FullPivLU a+    print $ rank FullPivLU <$> a      putStrLn "Here is a matrix whose columns form a basis of the null-space of A:"-    print $ kernel FullPivLU a+    print $ kernel FullPivLU <$> a      putStrLn "Here is a matrix whose columns form a basis of the column-space of A:"-    print $ image FullPivLU a+    print $ image FullPivLU <$> a

@@ -1,51 +1,66 @@-{-# LANGUAGE RecordWildCards #-}-import Data.Eigen.Matrix as M-import Data.Eigen.LA-import Data.List as L+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-}+{-# LANGUAGE TypeApplications #-}++import Control.Applicative (liftA2)+import Eigen.Matrix as M+import Eigen.Solver.LA+import Data.List as List import Control.Monad+import Data.Maybe (fromMaybe) -main = do-    let-        a :: MatrixXd-        a = fromList [-            [1,3.02, 6.89],-            [1,2.01, 5.39],-            [1,2.41, 6.01],-            [1,2.09, 5.55],-            [1,2.58, 6.32]]+a :: Maybe (MatrixXd 5 3)+a = fromList +  [ [1,3.02, 6.89]+  , [1,2.01, 5.39]+  , [1,2.41, 6.01]+  , [1,2.09, 5.55]+  , [1,2.58, 6.32]+  ] -        b = fromList $ Prelude.map return [-4.32,-3.79,-4.01,-3.86,-4.10]+b :: Maybe (MatrixXd 5 1)+b = fromList $ List.map return [-4.32,-3.79,-4.01,-3.86,-4.10] -    print a-    print b-    forM_ [FullPivLU, HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR, JacobiSVD] $ \d -> do-        let x = solve d a b-            e = relativeError x a b-            e' = norm (a*x - b) / norm b-        putStrLn $ replicate 20 '*'-        print d-        print x-        print e-        print e'+main :: IO ()+main = do+  print a+  print b+  +  forM_ [FullPivLU, HouseholderQR, ColPivHouseholderQR, FullPivHouseholderQR, JacobiSVD] $ \d -> do+    let x :: Maybe (MatrixXd 3 1)+        x = liftA2 (solve d) a b+    let e :: Maybe Double+        e = relativeError <$> x <*> a <*> b+        e' :: Maybe Double+        e' = liftA2 (/) (norm <$> (sub <$> (liftA2 mul a x) <*> b)) (norm <$> b)+    putStrLn $ replicate 20 '*'+    --print d+    print x+    print e+    print e'     putStrLn "\n-2.34666 - 0.25349 x1 - 0.174965 x2"     putStrLn "done" -    print $ (identity 4 4 :: MatrixXd)-    print $ M.normalize a-    print $ M.transpose a--    let-        a :: MatrixXd-        a = M.fromList [[0.68,  0.597,  -0.33],[-0.211,  0.823,  0.536],[ 0.566, -0.605, -0.444]]-        b = M.inverse a-    print a-    print b-    print $ a * b-    print $ linearRegression [-            [-4.32, 3.02, 6.89],-            [-3.79, 2.01, 5.39],-            [-4.01, 2.41, 6.01],-            [-3.86, 2.09, 5.55],-            [-4.10, 2.58, 6.32]]+  print $ (identity :: MatrixXd 4 4)+  print $ M.normalize <$> a+  print $ M.transpose <$> a +  let _a :: Maybe (MatrixXd 3 3)+      _a = M.fromList+        [ [ 0.680,  0.597, -0.330]+        , [-0.211,  0.823,  0.536]+        , [ 0.566, -0.605, -0.444]+        ]+      _b = M.inverse <$> _a+  +  print _a+  print _b+  print $ liftA2 M.mul _a _b +  print $ linearRegression (M.Row @5)+    [ [-4.32, 3.02, 6.89]+    , [-3.79, 2.01, 5.39]+    , [-4.01, 2.41, 6.01]+    , [-3.86, 2.09, 5.55]+    , [-4.10, 2.58, 6.32]+    ]

@@ -1,24 +1,29 @@-import qualified Data.Eigen.Matrix as M-import Data.Eigen.SparseMatrix-import Data.Eigen.SparseLA as LA+{-# LANGUAGE DataKinds #-}++import qualified Eigen.Matrix as M+import Eigen.SparseMatrix+import Eigen.Solver.SparseLA as LA import Control.Monad.Trans+import Data.Maybe -main = do-    let-        a :: SparseMatrixXd-        b :: SparseMatrixXd-        a = fromDenseList [[1,2,3], [4,5,6], [7,8,10]]-        b = fromDenseList [[3],[3],[4]]-    putStrLn "Here is the matrix A:"-    print $ a+matrices :: Maybe (SparseMatrixXd 3 3, SparseMatrixXd 3 1)+matrices = do+  a <- fromDenseList [[1,2,3], [4,5,6], [7,8,10]]+  b <- fromDenseList [[3],[3],[4]]+  pure (a, b) -    putStrLn "Here is the vector b:"-    print $ b+main :: IO ()+main = flip (maybe (print "WRONG! WRONG! WRONG!")) matrices $ \(a,b) -> do+  putStrLn "Here is the matrix A:"+  print a -    runSolverT (SparseLU COLAMDOrdering) $ do-        compute a-        x <- solve b-        info >>= lift.print-        determinant >>= lift . print-        lift $ putStrLn "The solution is:"-        lift $ print x+  putStrLn "Here is the vector b:"+  print b++  runSolverT (SparseLU COLAMDOrdering) $ do+    compute a+    x <- solve b+    info >>= lift . print+    determinant >>= lift . print+    lift $ putStrLn "The solution is:"+    lift $ print x

@@ -1,17 +1,15 @@-import Data.Eigen.Matrix-import Data.Eigen.LA+{-# LANGUAGE DataKinds #-}+{-# LANGUAGE ScopedTypeVariables #-} -main = do-    let-        a :: MatrixXd-        a = fromList [[1,2,3], [4,5,6], [7,8,10]]-        b = fromList [[3],[3],[4]]-        x = solve ColPivHouseholderQR a b-    putStrLn "Here is the matrix A:"-    print a+import Eigen.Matrix+import Eigen.Solver.LA+import Data.Maybe (fromMaybe) -    putStrLn "Here is the vector b:"-    print b+solution :: Maybe (MatrixXd 3 1)+solution = do+  _a :: MatrixXd 3 3 <- fromList [[1,2,3], [4,5,6], [7,8,9]]+  _b :: MatrixXd 3 1 <- fromList [[3],[3],[4]]+  pure $ solve ColPivHouseholderQR _a _b -    putStrLn "The solution is:"-    print x+main :: IO ()+main = print solution