diff --git a/Data/Eigen/Internal.hsc b/Data/Eigen/Internal.hsc --- a/Data/Eigen/Internal.hsc +++ b/Data/Eigen/Internal.hsc @@ -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" diff --git a/Data/Eigen/Matrix.hs b/Data/Eigen/Matrix.hs --- a/Data/Eigen/Matrix.hs +++ b/Data/Eigen/Matrix.hs @@ -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 diff --git a/Data/Eigen/SparseLA.hs b/Data/Eigen/SparseLA.hs --- a/Data/Eigen/SparseLA.hs +++ b/Data/Eigen/SparseLA.hs @@ -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 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 diff --git a/Data/Eigen/SparseMatrix.hs b/Data/Eigen/SparseMatrix.hs --- a/Data/Eigen/SparseMatrix.hs +++ b/Data/Eigen/SparseMatrix.hs @@ -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: +@ +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 . -} 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 diff --git a/Data/Eigen/SparseMatrix/Mutable.hs b/Data/Eigen/SparseMatrix/Mutable.hs new file mode 100644 --- /dev/null +++ b/Data/Eigen/SparseMatrix/Mutable.hs @@ -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 + diff --git a/cbits/eigen-sparse-la.cpp b/cbits/eigen-sparse-la.cpp --- a/cbits/eigen-sparse-la.cpp +++ b/cbits/eigen-sparse-la.cpp @@ -292,14 +292,6 @@ API_SPARSE_QR(sparse_la_rank, (int code, int s, void* p, double* x), (p,x)); template -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 RET sparse_la_setSymmetric(void* p, int x) { ((S*)p)->isSymmetric(x); return 0; diff --git a/cbits/eigen-sparse.cpp b/cbits/eigen-sparse.cpp --- a/cbits/eigen-sparse.cpp +++ b/cbits/eigen-sparse.cpp @@ -3,6 +3,22 @@ #include template +RET sparse_new(int rows, int cols, void** pr) { + typedef SparseMatrix M; + *(M**)pr = new M(rows, cols); + return 0; +} +API(sparse_new, (int code, int rows, int cols, void** pr), (rows,cols,pr)); + +template +RET sparse_clone(void* p, void** q) { + typedef SparseMatrix M; + *(M**)q = new M(*(M*)p); + return 0; +} +API(sparse_clone, (int code, void* p, void** q), (p,q)); + +template RET sparse_fromList(int rows, int cols, void* data, int size, void** pr) { typedef SparseMatrix M; typedef Triplet E; @@ -98,6 +114,13 @@ } API(sparse_coeff, (int code, void* p, int row, int col, void* pr), (p, row, col, pr)); +template +RET sparse_coeffRef(void* p, int row, int col, void** pr) { + *(T**)pr = &((SparseMatrix*)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 \ 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 +RET sparse_values(void* p, int* s, void** q) { + SparseMatrix* m = (SparseMatrix*)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 +RET sparse_outerStarts(void* p, int* s, void** q) { + SparseMatrix* m = (SparseMatrix*)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 +RET sparse_innerIndices(void* p, int* s, void** q) { + SparseMatrix* m = (SparseMatrix*)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 +RET sparse_innerNNZs(void* p, int* s, void** q) { + SparseMatrix* m = (SparseMatrix*)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 \ +RET sparse_##name(void* p) {\ + ((SparseMatrix*)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 +RET sparse_reserve(void* p, int s) { + ((SparseMatrix*)p)->reserve(s); + return 0; +} +API(sparse_reserve, (int code, void* p, int s), (p,s)); + +template +RET sparse_resize(void* p, int r, int c) { + ((SparseMatrix*)p)->resize(r,c); + return 0; +} +API(sparse_resize, (int code, void* p, int r, int c), (p,r,c)); + +template +RET sparse_conservativeResize(void* p, int r, int c) { + ((SparseMatrix*)p)->conservativeResize(r,c); + return 0; +} +API(sparse_conservativeResize, (int code, void* p, int r, int c), (p,r,c)); diff --git a/eigen.cabal b/eigen.cabal --- a/eigen.cabal +++ b/eigen.cabal @@ -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 . @@ -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 and . + . * Its ecosystem of 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 , when that is appropriate. - * is performed for SSE 2/3/4, ARM NEON, and AltiVec instruction sets, with graceful fallback to non-vectorized code. + . + * 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 and are available. + . * Eigen is thoroughly tested through its own (over 500 executables), the standard BLAS test suite, and parts of the LAPACK test suite. . Eigen is elegant. . * The API is extremely 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 diff --git a/test/solve-sparse.hs b/test/solve-sparse.hs new file mode 100644 --- /dev/null +++ b/test/solve-sparse.hs @@ -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