12

I have coded two implementations of an algorithm to calculate all the eigenvalues and eigenvectors of a symmetric matrix. One implementation uses the REPA library

https://github.com/felipeZ/Haskell-abinitio/blob/master/Science/QuantumChemistry/NumericalTools/JacobiMethod.hs

While the other implementation uses mutable Vectors and the ST monad

https://github.com/felipeZ/Haskell-abinitio/blob/master/Science/QuantumChemistry/NumericalTools/JacobiMethodST.hs

The algorithm is the Jacobi method which description can be found at http://en.wikipedia.org/wiki/Jacobi_eigenvalue_algorithm

I have tested the two implementations using 100 matrices of dimension 100 x 100, running the code sequentially and I have found the following times:

                  REPA      Mutable Vectors

Total time(s)    6.7            28.5

GC (s)           0.2             1.2

The Jacobi Algorithm requires an iterative updating of some of the matrix entries, meaning that most of the matrix is left untouched between iterations. Therefore I would have guessed (erroneously) that the cost of copying a new matrix for each iteration in the REPA implementation will be bigger than the cost of mutating the matrix using the monad ST, because as far as I understand REPA doesn't mutate the array but makes a copy of it.

It is REPA fusing all the operation and avoiding to copy of a new array in each iteration? or it is something else?

Can someone please comment on this result?

felipez
  • 1,212
  • 9
  • 21
  • 1
    I don't think the repa developer(s) are very active these days. I'd be curious to see how well yarr (https://hackage.haskell.org/package/yarr) performs but sadly I don't have time at the moment to investigate this myself. – idontgetoutmuch Mar 28 '15 at 18:04
  • 1
    Could you please describe your implementation in more detail? In particular, _repa_ distinguishes unboxed `U` and delayed `D` arrays, and mixing these two properly seems to be important for good performance. I wonder which ones and where you used. -- For having fast `ST` code it's often important to use tight loops and prevent memory allocation inside (which I guess _repa_ takes care of internally). And I guess copying data around a few times isn't that costly on modern machines, provided they're in one block to take advantage of memory caching. – Petr Jan 05 '16 at 17:42

0 Answers0