Lots of small matrices

[tera]

I didn’t implement the final micro-optimization (which is unpredictable when there’s two equal pivots in the column?).

Yes, behavior is undefined in that case. I am undecided myself, as for my own codes I certainly insist on reproducible results. For the sz=32 case results seemed to be stable though, and as the code looks so much nicer I decided to present it here.

I found using blockDim.y=5 worked best for the 32x32 case, so I stuck with that. Altering it to blockDim.y = 6 or 7 gave somewhat better results for N=64 (1.58s vs 1.81s), but I didn’t feel inclined to tune for so many cases.

Yeah, as the shared memory constraint allows fewer and fewer blocks per SM, we can have more and more warps per block without (further) limiting the number of blocks.

The general trend seems to be that the GPU has a big advantage for the smaller matrix sizes, but this vanishes at bigger matrix sizes (somewhere around N=50 the two are even). I’m guessing this reflects shared memory resources per multiprocessor getting strained by the larger matrix sizes. Would this be correct?

Yes, I’d think so too. The CPU is in advantage there because of its larger caches.

My understanding at right now is that the GPU excels at enormous matrices and —ironically— really small matrices. For somewhat larger matrices the advantage of one over the other starts becoming unclear.

I think this reflects the memory hierarchy: At huge sizes, the GPU excels through its about 10x larger memory bandwidth. At small sizes, the multiple SMs (and execution units) provide a large bandwidth advantage over the CPU as well, but there is an intermediate regime where the matrices still fit into CPU caches, but not into on-chip memory on the GPU.

__shared__ float A[sz][sz + !(sz%16)];

    /* !(sz%16) is added to avoid bank conflicts during pivoting

       when the matrix size is a multiple of 16 */

Interesting generalization. I had anticipated having a table of optimal padding sizes, but using [font=“Courier New”]sz | 1[/font] seems to be enough to avoid bank conflicts for arbitrary size.

I’ve found that exchanging rows explicitly (while still keeping the piv vector so we can return it) is slightly beneficial. It would also be advantageous for my proposal below.

We see a pattern that one matrix per block starts to loose in performance as the matrices becomes larger.

Indeed, as matrix size increases occupancy sucks more and more.

At this point it might be interesting to change the parallelization scheme a bit. How about going towards a few (1-4, whatever is best ) as the matrices become larger and reach a critical size ?

I don’t think that an intermediate scheme between one matrix per block and one matrix per thread would buy us much, unless one introduces really clever memory management (one might reuse the top rows which have finished processing as soon as they are enough to hold a full new matrix).

I’ve been thinking about using a blocked algorithm, or using mixed memory schemes (use registers or local memory to hold part of the matrix). But all seems to get ugly pretty quickly.

Unless one uses Fermi’s unified pointers to transparently store the lower, more often accessed part of the matrix in shared and the upper part in local memory.