proc sparseDiag (ref X: [?Xdom] ?eltype, in p: int, const val: eltype = 0)
                                                        where isSparseArr(X) {
  if Adom.rank != 2 then
    halt("Wrong rank for matrix inverse");
  var start, end = 0;
  if (p >= 0) { // upper or main diagonal
    start = 1;
    end = Xdom.shape(1) - p;
  }
  else { // lower diagonal
    start = 1 - p;
    end = Xdom.shape(1);
  }
  var indices : [start..end] (Xdom.idxType, Xdom.idxType);
  forall ind in {start..end} {
    indices[ind] = (ind, ind+p);
  }
  D.bulkAdd(indices, dataSorted=true, isUnique=true, preserveInds=false);
  coforall l in Locales do on l {
    const localDomain = X.localSubdomain();
    var (low_i, low_j) = localDomain.low;
    var (high_i, high_j) = localDomain.high;
    var localindices = indices;
    forall ij in localindices[low_i..high_i] {
      // Prevent unnecessary cross-locale communication
      if (localDomain.contains(ij)) {
         X(ij) = val;
      }
    }
  }
}