`each_cache_line(A,R)` iteration for efficient multidimensional algorithms

Working on AxisAlgorithms led me to wonder if we might benefit from yet more fancy iterators. Let's consider writing a matrix-multiplication algorithm that's conceptually equivalent to (A*B')', i.e., multiplication along rows of B. Obviously, A_mul_Bt does all of this except the final transpose, but take this as a simple prototype for a multidimensional generalization, e.g., multiplying along dimension 4 of a 6-dimensional array. In all of this, A is a matrix.

It seems that the way to write a tiled, high-performance algorithm would be something like this (here C is the output array):

for I2 in R2  # R2 is a CartesianRange covering the "trailing" dimensions of B, e.g., dimensions 5 and 6
    for tilek = 1:8:size(B, dim)   # dim is the multiplication dimension
        for tilei = 1:8:size(C, dim)
            for I1 in each_cache_line(B, R1)   # R1 is a CartesianRange covering the "leading" dimensions 1-3
                for k = tilek:min(tilek+7,size(B,dim)), i = tilei:min(tilei+7,size(C,dim)
                    ...

I haven't thought this through very deeply and the above could be flawed, but the main point is that since you have to reuse the same values repeatedly that are scattered over various non-adjacent spots in memory, it makes sense to iterate over chunks that have the alignment and size of a cache line.

There are numerous issues to tackle, for example how should this generalize to StridedArrays, etc. Tackling this issue would be an ambitious task, but one that might yield sizable performance benefits.

This might be fun for anyone who thinks we should implement BLAS operations in julia. Over in my unpublished KernelTools.jl package, I found that a well-tuned matmatmul, at least for Float32, could get within spitting distance of single-threaded BLAS computations (largely thanks to @simd). So it's not a crazy thing to think about.

CC @ArchRobison, @lindahua