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
```

首先使用 Distributed、DistributedArrays 和 LinearAlgebra 运行提供的 Julia 基准测试，然后跟踪矩阵乘法使用的 distribute path。该 issue 提到了通信、临时分配、event-loop、socket 和共享内存方面的问题，但没有指出源文件或测试。要视为完成，需要针对瓶颈进行范围明确的修复，并提供在相关矩阵大小和系统上的基准测试证据。

由索引模型根据 Issue 内容生成。