public
Last active

Sparse matrix-vector multiplication in Haskell and Fortran 77

  • Download Gist
CRS.f
FORTRAN
1 2 3 4 5 6 7 8 9 10 11 12
! 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;
CRS.hs
Haskell
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
{-# 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

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.