massiv icon indicating copy to clipboard operation
massiv copied to clipboard

Support for Sparse matrices

Open lehins opened this issue 5 years ago • 1 comments

I believe it is almost a requirement for an array library to have support for sparse matrices, so they need to be implemented in massiv as well.

With a little bit of research that I have done in this area, so far I envision it as another array representation:

data CSR r = CSR r

data instance Array (CSR r) Ix2 e = CSRArray
  { csrSize :: Ix2
  , csrIA :: ByteArray
  , csrJA :: ByteArray
  , csrData :: Array r Ix1 e }

Or something along these lines. Similarly CSC representation could be described.

lehins avatar Nov 02 '18 12:11 lehins

How difficult would it be to have different underlying representations of the sparse matrix? It's often easiest to have your in-memory Haskell-side representation be a match for some external library so that your FFI work is simpler and more performant.

There could be a few built in, tagged by different CS... or there could they be user specified somehow?

For example, the SuiteSparse library (natively) uses compressed-column representation:

(From the Cholmod documentation): "In the basic type (“packed,” which corresponds to how MATLAB stores its sparse matrices), an nrow-by-ncol matrix with nzmax entries is held in three arrays: p of size ncol+1, i of size nzmax, and x of size nzmax. Row indices of nonzero entries in column j are held in i [p[j] ... p[j+1]-1], and their corresponding numerical values are held in x [p[j] ... p[j+1]-1]. The first column starts at location zero (p[0]=0). There may be no duplicate entries. Row indices in each column may be sorted or unsorted (the A->sorted flag must be false if the columns are unsorted). The A->stype determines the storage mode: 0 if the matrix is unsymmetric, 1 if the matrix is symmetric with just the upper triangular part stored, and -1 if the matrix is symmetric with just the lower triangular part stored. In “unpacked” form, an additional array nz of size ncol is used. The end of column j in i and x is given by p[j]+nz[j]. Columns not need be in any particular order (p[0] need not be zero), and there may be gaps between the columns"

adamConnerSax avatar Jul 30 '20 13:07 adamConnerSax