CRS.f

## CRS.f
! 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
{-# 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