Matrix-Matrix multiply is quite slow

While looking with @yingboma into getting a PDE solved just by using DArray we encountered that matrix-matrix multiply is quite slow in DArray.

From discussion with @andreasnoack

1. `Distributed` currently has no `fetch!` e.g. a fetch into a localarray, it is therefore hard to avoid temporaries when working across processes. This causes many copies and requires GC work which creates communication bottlenecks.

2. Our communication layer doesn't support RDMA so there are copies happening in the network-stack, and we use sockets instead of shared memory for commincation on the same node.

3. There are some communications bottlenecks due to how we use the event-loop and it is feasible to to get into a situation where forward progress is hard to make due a machine being busy with computation and not communicating in a timely fashion.

```julia
using Distributed
addprocs(4)

using LinearAlgebra

# Set worker BLAS to one thread onlye
@sync for p in workers()
@async remotecall_wait(LinearAlgebra.BLAS.set_num_threads, p , 1)
end

using BenchmarkTools
using DistributedArrays

const suite = BenchmarkGroup()
suite["Array"] = BenchmarkGroup()
suite["distribute"] = BenchmarkGroup()

function benchmark(T=Array, N=10)
@benchmarkable A * B setup=(A = $T(rand($N, $N)); B = $T(rand($N, $N)))
end

for N in (2^i for i = 5:13)
suite["Array"][N] = benchmark(Array, N)
suite["distribute"][N] = benchmark(distribute, N)
end

tune!(suite)
results = run(suite)
```

I would be interested in gathering numbers from different systems here. My first set of results is from just my local laptop with 2 Cores - 4 Threads and using 4 Julia processes.

There is a lot of overhead for smallish problems, but the results aren't that bad once we get to interesting problem sizes...

```
julia> for (name, trial) in sort(collect(results["Array"]), by=x->time(x[2]))
t = time(trial) / 1e6
println(rpad(name, 25, "."), lpad(string(round(t, digits=2), " ms"), 20, "."))
end
32.....................................0.0 ms
64....................................0.02 ms
128...................................0.06 ms
256...................................0.41 ms
512...................................3.51 ms
1024.................................24.51 ms
2048................................249.21 ms
4096...............................2226.76 ms
8192..............................18990.07 ms

julia> for (name, trial) in sort(collect(results["distribute"]), by=x->time(x[2]))
t = time(trial) / 1e6
println(rpad(name, 25, "."), lpad(string(round(t, digits=2), " ms"), 20, "."))
end
32....................................2.01 ms
64....................................2.06 ms
128...................................2.32 ms
256...................................2.97 ms
512...................................6.63 ms
1024..................................34.2 ms
2048................................261.15 ms
4096...............................2295.89 ms
8192..............................17112.45 ms
```

Start by running the supplied Julia benchmark with Distributed, DistributedArrays, and LinearAlgebra, then trace the distribute path used by matrix multiplication. The issue names communication, temporary allocation, event-loop, socket, and shared-memory concerns but no source files or tests. Done would require a scoped bottleneck fix with benchmark evidence across relevant matrix sizes and systems.

Written by the indexing model from the issue text.