Skip to content

Instantly share code, notes, and snippets.

@bos
Created April 7, 2011 01:09
Show Gist options
  • Star 2 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save bos/906853 to your computer and use it in GitHub Desktop.
Save bos/906853 to your computer and use it in GitHub Desktop.
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