packages feed

eigen 2.1.3 → 2.1.4

raw patch · 9 files changed

+438/−63 lines, 9 filesdep ~basedep ~primitivedep ~vector

Dependency ranges changed: base, primitive, vector

Files

Data/Eigen/Internal.hsc view
@@ -154,6 +154,8 @@  #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"@@ -165,11 +167,11 @@ #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_diagonal,      "CSparseMatrixPtr a 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"@@ -181,7 +183,19 @@ #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"@@ -205,7 +219,6 @@ #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_simplicialFactorize,"CSolverPtr a b -> CSparseMatrixPtr a b -> 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"
Data/Eigen/Matrix.hs view
@@ -148,7 +148,7 @@         "\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`+-- | Basic matrix math exposed through Num instance: @(*)@, @(+)@, @(-)@, `fromInteger`, `signum`, `abs`, `negate` instance I.Elem a b => Num (Matrix a b) where     (*) = mul     (+) = add@@ -441,24 +441,24 @@     -- | 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 +    | 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 } +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@ +-- | Upper trinagle of the matrix. Shortcut for @triangularView Upper@ upperTriangle :: I.Elem a b => Matrix a b -> Matrix a b upperTriangle = triangularView Upper
Data/Eigen/SparseLA.hs view
@@ -137,7 +137,6 @@     rank,     setPivotThreshold,     -- * SparseLU Solver-    simplicialFactorize,     setSymmetric,     matrixL,     matrixU,@@ -169,14 +168,14 @@     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 +    -- | 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:         @@@ -184,7 +183,7 @@         @         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@@ -193,7 +192,7 @@   class I.Code s => Solver s where--- | For direct methods, the solution are computed at the machine precision.+-- | 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`.@@ -237,11 +236,11 @@     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  \link OrderingMethods_Module the OrderingMethods module \endlink for the list of -    built-in and external ordering methods. +    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@@ -313,7 +312,7 @@     withForeignPtr fa $ \a ->         I.call $ I.sparse_la_compute i s a --- | An expression of the solution x of @A x = b@ using the current decomposition of @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 ->@@ -324,7 +323,7 @@         SM.SparseMatrix <$> FC.newForeignPtr x (I.call $ I.sparse_free x)  {---- | The solution @x@ of @A x = b@ using the current decomposition of @A@ and @x0@ as an initial solution.+-- | 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 ->@@ -336,7 +335,11 @@         SM.SparseMatrix <$> FC.newForeignPtr x (I.call $ I.sparse_free x) -} --- | Success if the iterations converged, and NoConvergence otherwise.+-- |+-- * `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) @@ -346,7 +349,7 @@  -- | 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+--   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 @@ -375,7 +378,7 @@ 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 ()@@ -385,21 +388,15 @@ rank :: (MonadIO m, I.Elem a b) => SolverT SparseQR a b m Int rank = _get_prop I.sparse_la_rank -simplicialFactorize :: (MonadIO m, I.Elem a b) => (SM.SparseMatrix a b) -> SolverT SparseLU a b m ()-simplicialFactorize (SM.SparseMatrix fa) = ask >>= \(i,fs) -> liftIO $-    withForeignPtr fs $ \s ->-    withForeignPtr fa $ \a ->-        I.call $ I.sparse_la_simplicialFactorize i s a- -- | 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+-- | 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+-- | 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 @@ -408,13 +405,13 @@ 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
Data/Eigen/SparseMatrix.hs view
@@ -5,12 +5,16 @@  module Data.Eigen.SparseMatrix (     -- * SparseMatrix type-    -- | SparseMatrix aliases follows Eigen naming convention     SparseMatrix(..),     SparseMatrixXf,     SparseMatrixXd,     SparseMatrixXcf,     SparseMatrixXcd,+    -- * Matrix internal data+    values,+    innerIndices,+    outerStarts,+    innerNNZs,     -- * Accessing matrix data     cols,     rows,@@ -28,9 +32,6 @@     squaredNorm,     blueNorm,     block,-    compress,-    uncompress,-    compressed,     nonZeros,     innerSize,     outerSize,@@ -43,9 +44,18 @@     scale,     transpose,     adjoint,+    -- * Matrix representation+    compress,+    uncompress,+    compressed,     -- * Matrix serialization     encode,     decode,+    -- * Mutable matricies+    thaw,+    freeze,+    unsafeThaw,+    unsafeFreeze, ) where  import qualified Prelude as P@@ -65,6 +75,7 @@ #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@@ -75,19 +86,52 @@  {-| A versatible sparse matrix representation. -This class implements a more versatile variants of the common compressed row/column storage format.-Each colmun's (resp. row) non zeros are stored as a pair of value with associated row (resp. colmiun) index.-All the non zeros are stored in a single large buffer.-Unlike the compressed format, there might be extra space inbetween the nonzeros of two successive colmuns-(resp. rows) such that insertion of new non-zero can be done with limited memory reallocation and copies.+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. -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.+It consists of four compact arrays: -A call to the function 'compress' turns the matrix into the standard compressed format compatible with many library.+* `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. -Implementation deails of SparseMatrix are intentionally hidden behind ForeignPtr bacause Eigen doesn't provide mapping over plain data for sparse matricies.+This storage scheme is better explained on an example. The following matrix -For more infomration please see Eigen documentation page: <http://eigen.tuxfamily.org/dox/classEigen_1_1SparseMatrix.html>+@+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@@ -108,7 +152,7 @@         "SparseMatrix ", show (rows m), "x", show (cols m),         "\n", L.intercalate "\n" $ P.map (L.intercalate "\t" . P.map show) $ toDenseList m, "\n"] --- | Shortcuts for basic matrix math+-- | Basic sparse matrix math exposed through Num instance: @(*)@, @(+)@, @(-)@, `fromInteger`, `signum`, `abs`, `negate` instance I.Elem a b => Num (SparseMatrix a b) where     (*) = mul     (+) = add@@ -125,6 +169,26 @@ 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) +-- | 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)@@ -166,15 +230,15 @@ nonZeros :: I.Elem a b => SparseMatrix a b -> Int nonZeros = _unop I.sparse_nonZeros (return . I.cast) --- | Turns the matrix into the compressed format+-- | The matrix in the compressed format compress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b compress = _unop I.sparse_makeCompressed mk --- | not exposed currently+-- | The matrix in the uncompressed mode uncompress :: I.Elem a b => SparseMatrix a b -> SparseMatrix a b uncompress = _unop I.sparse_uncompress mk --- | not exposed currently+-- | Is this in compressed form? compressed :: I.Elem a b => SparseMatrix a b -> Bool compressed = _unop I.sparse_isCompressed (return . (/=0)) @@ -315,6 +379,22 @@         I.call $ I.sparse_fromList rows cols p size pq         peek pq >>= mk +-- | 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 ->@@ -329,3 +409,20 @@             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
+ Data/Eigen/SparseMatrix/Mutable.hs view
@@ -0,0 +1,144 @@+{-# 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+
cbits/eigen-sparse-la.cpp view
@@ -292,14 +292,6 @@ API_SPARSE_QR(sparse_la_rank, (int code, int s, void* p, double* x), (p,x));  template <class T, class M, class S>-RET sparse_la_simplicialFactorize(void* p, void* a) {-    ((S*)p)->simplicialfactorize(*(M*)a);-    return 0;-}-API_SPARSE_LU(sparse_la_simplicialFactorize, (int code, int s, void* p, void* a), (p,a));---template <class T, class M, class S> RET sparse_la_setSymmetric(void* p, int x) {     ((S*)p)->isSymmetric(x);     return 0;
cbits/eigen-sparse.cpp view
@@ -3,6 +3,22 @@ #include <Eigen/LeastSquares>  template <class T>+RET sparse_new(int rows, int cols, void** pr) {+    typedef SparseMatrix<T> M;+    *(M**)pr = new M(rows, cols);+    return 0;+}+API(sparse_new, (int code, int rows, int cols, void** pr), (rows,cols,pr));++template <class T>+RET sparse_clone(void* p, void** q) {+    typedef SparseMatrix<T> M;+    *(M**)q = new M(*(M*)p);+    return 0;+}+API(sparse_clone, (int code, void* p, void** q), (p,q));++template <class T> RET sparse_fromList(int rows, int cols, void* data, int size, void** pr) {     typedef SparseMatrix<T> M;     typedef Triplet<T> E;@@ -98,6 +114,13 @@ } API(sparse_coeff, (int code, void* p, int row, int col, void* pr), (p, row, col, pr)); +template <class T>+RET sparse_coeffRef(void* p, int row, int col, void** pr) {+    *(T**)pr = &((SparseMatrix<T>*)p)->coeffRef(row, col);+    return 0;+}+API(sparse_coeffRef, (int code, void* p, int row, int col, void** pr), (p, row, col, pr));+ #define SPARSE_PROP(name,type)\ template <class T>\ RET sparse_##name(void* p, void* pr) {\@@ -153,6 +176,7 @@                 dst->insert(row, col) = val;         }     }+    dst->makeCompressed();     *(M**)pq = dst.release();     return 0; }@@ -165,3 +189,73 @@     return 0; } API(sparse_toMatrix, (int code, void* p, void* q, int rows, int cols), (p,q,rows,cols));++template <class T>+RET sparse_values(void* p, int* s, void** q) {+    SparseMatrix<T>* m = (SparseMatrix<T>*)p;+    *s = m->outerIndexPtr()[m->outerSize()];+    *q = m->valuePtr();+    return 0;+}+API(sparse_values, (int code, void*p, int* s, void** q), (p,s,q));++template <class T>+RET sparse_outerStarts(void* p, int* s, void** q) {+    SparseMatrix<T>* m = (SparseMatrix<T>*)p;+    *s = m->outerSize() + 1;+    *q = m->outerIndexPtr();+    return 0;+}+API(sparse_outerStarts, (int code, void*p, int* s, void** q), (p,s,q));++template <class T>+RET sparse_innerIndices(void* p, int* s, void** q) {+    SparseMatrix<T>* m = (SparseMatrix<T>*)p;+    *s = m->outerIndexPtr()[m->outerSize()];+    *q = m->innerIndexPtr();+    return 0;+}+API(sparse_innerIndices, (int code, void*p, int* s, void** q), (p,s,q));++template <class T>+RET sparse_innerNNZs(void* p, int* s, void** q) {+    SparseMatrix<T>* m = (SparseMatrix<T>*)p;+    *s = m->outerSize();+    *q = m->innerNonZeroPtr();+    return 0;+}+API(sparse_innerNNZs, (int code, void*p, int* s, void** q), (p,s,q));++#define SPARSE_INPLACE(name,op)\+template <class T>\+RET sparse_##name(void* p) {\+    ((SparseMatrix<T>*)p)->op();\+    return 0;\+}\+API(sparse_##name, (int code, void* p), (p));++SPARSE_INPLACE(setIdentity, setIdentity);+SPARSE_INPLACE(setZero, setZero);+SPARSE_INPLACE(compressInplace, makeCompressed);+SPARSE_INPLACE(uncompressInplace, uncompress);++template <class T>+RET sparse_reserve(void* p, int s) {+    ((SparseMatrix<T>*)p)->reserve(s);+    return 0;+}+API(sparse_reserve, (int code, void* p, int s), (p,s));++template <class T>+RET sparse_resize(void* p, int r, int c) {+    ((SparseMatrix<T>*)p)->resize(r,c);+    return 0;+}+API(sparse_resize, (int code, void* p, int r, int c), (p,r,c));++template <class T>+RET sparse_conservativeResize(void* p, int r, int c) {+    ((SparseMatrix<T>*)p)->conservativeResize(r,c);+    return 0;+}+API(sparse_conservativeResize, (int code, void* p, int r, int c), (p,r,c));
eigen.cabal view
@@ -1,5 +1,5 @@ name:           eigen-version:        2.1.3+version:        2.1.4 homepage:       https://github.com/osidorkin/haskell-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>.@@ -7,25 +7,33 @@                 Eigen is versatile.                 .                 * It supports all matrix sizes, from small fixed-size matrices to arbitrarily large dense matrices, and even sparse matrices.+                .                 * It supports all standard numeric types, including std::complex, integers, and is easily extensible to custom numeric types.+                .                 * It supports various <http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html matrix decompositions> and <http://eigen.tuxfamily.org/dox/group__TutorialGeometry.html geometry features>.+                .                 * Its ecosystem of <http://eigen.tuxfamily.org/dox/unsupported/index.html unsupported modules> provides many specialized features such as non-linear optimization, matrix functions, a polynomial solver, FFT, and much more.                 .                 Eigen is fast.                 .                 * Expression templates allow to intelligently remove temporaries and enable <http://eigen.tuxfamily.org/dox/TopicLazyEvaluation.html lazy evaluation>, when that is appropriate.-                * <http://eigen.tuxfamily.org/index.php?title=FAQ#Vectorization Explicit vectorization> is performed for SSE 2/3/4, ARM NEON, and AltiVec instruction sets, with graceful fallback to non-vectorized code.+                .+                * <http://eigen.tuxfamily.org/index.php?title=FAQ#Vectorization Explicit vectorization> is performed for SSE 2\/3\/4, ARM NEON, and AltiVec instruction sets, with graceful fallback to non-vectorized code.+                .                 * Fixed-size matrices are fully optimized: dynamic memory allocation is avoided, and the loops are unrolled when that makes sense.+                .                 * For large matrices, special attention is paid to cache-friendliness.                 .                 Eigen is reliable.                 .                 * Algorithms are carefully selected for reliability. Reliability trade-offs are <http://eigen.tuxfamily.org/dox/group__TopicLinearAlgebraDecompositions.html clearly documented> and <http://eigen.tuxfamily.org/dox/classEigen_1_1JacobiSVD.html extremely> <http://eigen.tuxfamily.org/dox/classEigen_1_1FullPivHouseholderQR.html safe> <http://eigen.tuxfamily.org/dox/classEigen_1_1FullPivLU.html decompositions> are available.+                .                 * Eigen is thoroughly tested through its own <http://eigen.tuxfamily.org/index.php?title=Tests test suite> (over 500 executables), the standard BLAS test suite, and parts of the LAPACK test suite.                 .                 Eigen is elegant.                 .                 * The API is extremely <http://eigen.tuxfamily.org/index.php?title=API_Showcase clean and expressive> while feeling natural to C++ programmers, thanks to expression templates.+                .                 * Implementing an algorithm on top of Eigen feels like just copying pseudocode.                 .                 Eigen has good compiler support as we run our test suite against many compilers to guarantee reliability and work around any compiler bugs.@@ -1292,6 +1300,7 @@                         Data.Eigen.Matrix                         Data.Eigen.Matrix.Mutable                         Data.Eigen.SparseMatrix+                        Data.Eigen.SparseMatrix.Mutable                         Data.Eigen.Parallel                         Data.Eigen.Internal @@ -1313,6 +1322,11 @@ Test-Suite test-solve     type:               exitcode-stdio-1.0     main-is:            test/solve.hs+    build-depends:      base, primitive, vector, eigen++Test-Suite test-solve-sparse+    type:               exitcode-stdio-1.0+    main-is:            test/solve-sparse.hs     build-depends:      base, primitive, vector, eigen  Test-Suite test-rank
+ test/solve-sparse.hs view
@@ -0,0 +1,24 @@+import qualified Data.Eigen.Matrix as M+import Data.Eigen.SparseMatrix+import Data.Eigen.SparseLA as LA+import Control.Monad.Trans++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++    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