Batched mapAccumL and mapAccumR

I have an array of band matrices and want to perform a batched Cholesky decomposition. I know statically the band-width and thus I can e.g. represent a symmetric 50x50 matrix with band-width 5
not only as 50x3 array of `Double`, but even more strict as 50-element array of `(Double,Double,Double)`.
The Cholesky decomposition then becomes a `mapAccumL`. Here is the running code implemented for Haskell lists:

```
cholesky3 :: (Fractional a) => [(a,a,a)] -> [(a,(a,a))]
cholesky3 =
snd .
mapAccumL
(\((d1,l31), (d2,(l32,l42))) (a33,a43,a53) ->
let d3 = a33 - l31^2 * d1 - l32^2 * d2
l43 = (a43 - l42*l32*d2) / d3
l53 = a53 / d3
c3 = (d3,(l43,l53))
in (((d2,l42), c3), c3))
((0,0), (0,(0,0)))
```

The input list contains the diagonal and sub-diagonal non-zero elements of each column and the output list contains the diagonal elements of the diagonal matrices and the sub-diagonal elements column by column.

Now I like to port that to Accelerate. Unfortunately Accelerate's variants of `scan` all require associative accumulator functions in order to get parallelization from splitting the scan. In contrast to that, I like to perform the `mapAccumL` sequentially and I want to obtain parallisation by running it on hundreds of matrices in parallel. That is, what I need is a function like this one:

```
mapAccumL ::
(Exp acc -> Exp x -> (Exp acc, Exp y)) ->
Exp acc ->
Acc (Array (ix :. Int) x) ->
Acc (Array (ix :. Int) y)
```

See the according post to the Accelerate mailing list: https://groups.google.com/forum/#!topic/accelerate-haskell/BoSWpm0FeJ4