Version 36 (modified by gckeller, 5 years ago) (diff)




Strict Arrays, Delayed Array and Shape Data Type




The Shape class


  dim   :: sh -> Int           
  -- ^number of dimensions (>= 0)
  size  :: sh -> Int           
  toIndex :: sh -> sh -> Int     
  -- ^corresponding index into a linear representation of 
  -- the array (first argument is the shape)

  fromIndex:: sh -> Int -> sh   
  -- ^given index into linear row major representation, calculates 
  -- index into array                               

  range      :: sh -> U.Array sh

  inRange      :: sh -> sh -> Bool

  zeroDim      :: sh

  intersectDim :: sh -> sh -> sh

  next:: sh -> sh -> Maybe sh


The following instances are defined

instance Shape () 



  addDim     :: sh -> sh' -> sh    
  modDim     :: sh -> sh' -> sh    
  inject     :: sh -> sh' -> sh


Representation, Order of Elements, and Lifted Values



this means that, for example


A scalar value cis isomorphic to a zero-dimensional array

  DArray () (\_ -> c)

Lifting a scalar value over a shape `dim':


Operations on Arrays and Delayed Arrays

Array Creation and Conversion


data Array dim e where
          , arrayData     :: U.Array e          -- flat parallel array
           }  -> Array dim e
  deriving Show



data DArray dim e where 

Delayed arrays can be converted to and from strict arrays:




Singular (zero-dimensional) arrays are isomorphic to scalar values and can be converted to one using the following function:

toScalar:: U.Elt e  => DArray () e -> DArray () e

Note that in contrast to all the previous operations, toScalar requires the array to be of a particular dimensionality.

Collection Oriented Operations

Elementwise Application of Functions






fold :: (U.Elt e, A.Shape dim) => 


Support for Parallel Execution

Since the implementation of DArrays builds on the DPH library, all the array operations can be executed in parallel. That is, we compiling a DArray program, the compiler generates a sequential as well as a parallel executable. All collective operations, like map, fold, and so on are executed in parallel.

Shape Polymorphic Computations on Arrays







data Index a initialDim projectedDim where
  IndexNil   :: Index a initialDim initialDim
  IndexAll   :: (Shape init, Shape proj) =>      


type SelectIndex = Index Int
type MapIndex    = Index ()

Given this definition, the type of select now becomes







arr' :: Array dim Double
arr' = select arr (IndexFixed 3 IndexNil)




simpleReplicate arr n =
  replicate arr (IndexFixed n IndexNil)


elementwiseReplicate:: (U.Elt e, Shape dim) => 
elementwiseReplicate arr n =
  replicate arr (IndexAll (IndexFixed n IndexNil))




`Nesting' Array Functions





and those functions can be trivially mapped since

 map f = f

The function fold, which we introduced earlier, is an example of a mappable library function. Applied to a matrix, fold (+) 0 will calculate the sum of all rows. If run in parallel, fold itself is run in parallel, and all the rows are processed in parallel.

fold :: (U.Elt e, A.Shape dim) => 


  selectEvenInd (DArray (sh :*: n ) f =

In this case dim could simply be unit, if and selectEven extracts all elements with an even index, or it could be any other shape, and thus

  map selectEvenInd = selectEvenInd

where selectEvenInd on the left and right-hand side of the equation are two different instances of the function. Now, lets try and write function selectEvenElems, which selects all even elements from an array. To determine the shape of the result, it is not sufficient to look at the shape of the argument. Instead, we have to calculate the new size by counting the even elements in the array using fold. The function fold is mappable, and returns an array. If the argument is a one-dimensional array, the result is a singular array, which then can be converted to a scalar using toScalar. The necessary application of toScalar, however, also restricts sh, which can now only be unit, and therefore the whole operation selectEvenElems is restricted to one-dimensional arrays, and not mappable.

  selectEvenElems (DArray  (sh :*: n) f) = DArray (sh :*: newSize) <......>
          newSize = toScalar $ fold (+) 0 $ map (\x -> if (even x) then 1 else 0) arr  

Example 1: Matrix multiplication


  fold (+) 0 arrDP



  fold (+) 0 arrDP
    arr2T = forceDArray $ transpose arr2

transpose (DArray (() :*: m :*: n) f) =
  DArray (() :*: n :*: m) (\(() :*: i :*: j) -> f (() :*: j :*: i))


As it is, mmMult can only take two-dimensional arrays as arguments, and is not mappable. If we look at the implementation closely, we can see that the restriction to two-dimensional arrays is unnecessary. All we have to do to generalise it is to adjust the type signatures and replace () with an arbitrary shape variable:

mmMult1:: Shape dim => 
mmMult1 arr1@(DArray (sh :*: m1 :*: n1) _) arr2@(DArray (sh' :*: m2 :*: n2) _) = 
  fold (+) 0 arrDP
    arr2T = forceDArray $ transpose arr2
    arrDP = DArray (sh:*: m1 :*: n2 :*:n1) 
       (\(sh :*: i :*: j :*: k) -> (index arr1 (sh :*: i :*: k)) * (index arr2T (sh :*: j:*: k)))

transpose:: Shape dim =>
transpose (DArray (sh:*: m :*: n) f) =
  DArray (sh :*: n :*: m) (\(sh :*: i :*: j) -> f (sh :*: j :*: i))

An alternative way to define matrix-matrix multiplication is in terms of the collective library functions provided. First, we expand both arrays and, in case of arr2 transpose it such that the elements which have to be multiplied match up. Then, we calculate the products using zipWith, and then use fold to compute the sums:

mmMult2:: (Array.RepFun dim, Array.InitShape dim, Array.Shape dim) => 
   fold (+) 0 (arr1Ext * arr2Ext)
    arr2T   = forceDArray $ transpose arr2  
    arr2Ext = replicate arr2T
                 (Array.IndexAll (Array.IndexAll (Array.IndexFixed n1 Array.IndexNil)))

In this implementation, transpose is necessary to place the elements at the right position for zipWith, and we call forceDArray for the same reason as in the previous implementation, to improve locality. Also, mmMult2' outperforms mmMult1, as the use of replicate` exposes the structure of the communication, whereas the general index calculations in mmMult1 hide this structure, and thus are less efficient.

Performance of Matrix-Matrix Multiplication

  size                                |       256 |       512 |    1024 | 
  mmMult1                             |       675 |      5323 |   42674 |
  mmMult2                             |       345 |      2683 |   21442 |
  mmMult2 (without forceDArray)       |       974 |      8376 |   73368 |
  C (hand written)                    |        34 |       514 |    7143 |
  C (MacOS Accelerated Vector ops)    |        33 |       510 |    6949 |

Open Questions

Do we need array comprehension on DArrays?

Wave computations