mmMult :: Monad m
=> Array U DIM2 Double
-> Array U DIM2 Double
-> m (Array U DIM2 Double)
mmMult a b = sumP (Repa.zipWith (*) aRepl bRepl)
where
t = transpose2D b
aRepl = extend (Z :.All :.colsB :.All) a
bRepl = extend (Z :.rowsA :.All :.All) t
(Z :.colsA :.rowsA) = extent a
(Z :.colsB :.rowsB) = extent b
But that is misleading because there are two kinds of arrays in Haskell, and the kind (namely, mutable) I used in my literal answer is the less common kind. In fact, it is so uncommon and so contrary to the spirit of Haskell that the language implementations do not take the trouble to make it fast. Last time I played with the Glasgow Haskell Compiler, for example, a related operation (readIORef) was over 100 times slower than its counterpart in C.
Parenthetically, because you need to use "a monad" to perform any sort of mutable operation in Haskell, the translation of
m[i,j] := m[j,i]
into Haskell is not
writeIOArray m (i,j) (readIOArray m (j,i))
like someone unfamiliar with monadic style might think, but rather
readIOArray m (j,i) >>= writeIOArray m (i,j)
or equivalently
do
element <- readIOArray m (j,i)
writeIOArray m (i,j) element
But that is a side issue. The central fact is that even though it is technically possible, it is considered bad style in Haskell to access or update individual array elements.
> The central fact is that even though it is technically possible, it is considered bad style in Haskell to access or update individual array elements.
And another reason why Haskell will not fair well. A lot of numerical mathematics is linear algebra which involves, almost exclusively, updating array elements.
I'm sure Haskell is good at some things but for many of the types of problems that numerical analysts face Haskell's paradigm just doesn't make sense. Especially when we have languages like Fortran/Matlab which were explicitly designed for computational mathematics.
Yeah, I was just thinking that as /u/warfangle/ said that
> with Haskell the code would look more like the actual equations you're working from
so at least in matrix algorithms, and finite difference methods for partial differential equations, the textbooks usually explain the algorithms using a notation that refers to matrix elements. So in these cases, an imperative implementation would be close to the math textbook representation, and the Haskell code would not.
I believe the Haskell way of doing numerical things would be thinking in terms of what you're trying to do rather than "Step x is get the value from [i,j] of the matrice". I could be very wrong, but for general programming I had to stop thinking in terms of accessing specific elements of an array to write idiomatic Haskell.
