Sparse matrix-vector multiplication in Haskell and Fortran 77
! From http://netlib.org/linalg/html_templates/node98.html | |
! | |
! val: non-zero values from the matrix | |
! col_ind: indices into val at which columns start in successive rows | |
! row_ptr: indices into val at which successive rows start | |
for i = 1, n | |
y(i) = 0 | |
for j = row_ptr(i), row_ptr(i+1) - 1 | |
y(i) = y(i) + val(j) * x(col_ind(j)) | |
end; | |
end; |
{-# LANGUAGE RecordWildCards #-} | |
import Data.Vector.Unboxed as U | |
-- | A compressed row storage (CRS) sparse matrix. | |
data CRS a = CRS { | |
crsValues :: Vector a | |
, colIndices :: Vector Int | |
, rowIndices :: Vector Int | |
} deriving (Show) | |
multiplyMV :: CRS Double -> Vector Double -> Vector Double | |
multiplyMV CRS{..} x = generate (U.length x) outer | |
where outer i = U.sum . U.map inner $ U.enumFromN start (end-start) | |
where inner j = (crsValues ! j) * (x ! (colIndices ! j)) | |
start = rowIndices ! i | |
end = rowIndices ! (i+1) | |
(!) a b = unsafeIndex a b |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
This comment has been minimized.
I like this better http://darcs.haskell.org/packages/dph/dph-examples/spectral/SMVM/dph/SMVMVectorised.hs