What's wrong with using Identity monad with mmultP when using repa?

Question

I do not understand why this program using repa:

import Data.Array.Repa
import Data.Array.Repa.Algorithms.Matrix
import Data.Functor.Identity

go = runIdentity $ do
  let mat = fromListUnboxed (ix2 2 2) [1..4]
  let ins = fromListUnboxed (ix2 2 1) [1, 1]
  mmultP mat ins

is giving me the following warning:

Data.Array.Repa: Performing nested parallel computation sequentially.
  You've probably called the 'compute' or 'copy' function while another
  instance was already running. This can happen if the second version
  was suspended due to lazy evaluation. Use 'deepSeqArray' to ensure
  that each array is fully evaluated before you 'compute' the next one.

I have no nested computations, I didn't call compute or copy, and everything that I used to do the computation is inside the same monad. Is it something to do with lazy evaluation? If so, how do I make the parallel computation happen while using the Identity monad (to keep the overall computation pure) ?

For reference, replacing runIdentity with runST makes it work, although in either case the specific monad's functionality isn't being used at all.

Answer 1

The reason for having the Monad constraint in computeP and similar parallel operations is to force sequential computation where required. This is described in [Parallel and Concurrent Programming in Haskell], in subsection Monads and computeP.

In your case, the problem seems to be caused by the internal implementation of mmultP:

mmultP  :: Monad m
        => Array U DIM2 Double 
        -> Array U DIM2 Double 
        -> m (Array U DIM2 Double)

mmultP arr brr 
 = [arr, brr] `deepSeqArrays` 
   do   trr      <- transpose2P brr
        let (Z :. h1  :. _)  = extent arr
        let (Z :. _   :. w2) = extent brr
        computeP 
         $ fromFunction (Z :. h1 :. w2)
         $ \ix   -> R.sumAllS 
                  $ R.zipWith (*)
                        (unsafeSlice arr (Any :. (row ix) :. All))
                        (unsafeSlice trr (Any :. (col ix) :. All))

It calls first transpose2P and then computeP, and transpose2P internally calls computeUnboxedP. If you use the Identity monad, there is no sequencing forced, so both these parallel computations can run in parallel, hence the nested parallelism.

If you want to keep things pure and also don't want to use ST, you can replace Identity with Eval, which is a strict version of Identity:

import Control.Parallel.Strategies
...
go = runEval $ do ...