Support runtime GEMM dimensions for dynamic shapes

README.md

@@ -96,7 +96,45 @@ sofieBLAS<alpaka::TagGpuHipRt> blas(queue);
 blas.matmul('N', 'N', size, size, size, 1.0f, dA, dB, 0.0f, dC);
 ```
 
-The GPU backends (`BlasCuda`, `BlasHip`) additionally expose `gemmrelu`/`gemmgelu` (fused bias + activation via cuBLASLt/hipBLASLt epilogues), `gemmStridedBatched`, and `addLayoutConfig` (used to pre-register cuBLASLt/hipBLASLt matrix layouts for a given shape before the first `matmul`/`gemm` call on that shape).
+The GPU backends (`BlasCuda`, `BlasHip`) additionally expose 
+- `gemmrelu`/`gemmgelu` (fused bias + activation via cuBLASLt/hipBLASLt epilogues)
+- `gemmStridedBatched` for batched gemm operations through strides
+- `addOperationConfig` that creates the matrix layouts and resolves the multiply algorithm for a call site's shape ahead of its first call (see below).
+
+## GEMM call instantiation and the algorithm cache
+
+A GEMM call computes `C = alpha * op(A) * op(B) + beta * C`, where A and B are the input matrices, C the output, and `op` an optional transpose. To run one, cuBLASLt and hipBLASLt need three kinds of objects besides the data:
+
+- a **matrix layout** per matrix: a descriptor holding its rows, columns and leading dimension;
+- a **matmul descriptor**: the operation settings (the transposes and the epilogue);
+- an **algorithm**: the concrete GEMM kernel the library selects for the given settings and dimensions, obtained by querying its heuristic (`cublasLtMatmulAlgoGetHeuristic` / `hipblasLtMatmulAlgoGetHeuristic`). The query runs on the host and is not free.
+
+The CUDA backend (`BlasCuda`, over cuBLASLt) and the HIP backend (`BlasHip`, over hipBLASLt) behave identically: all three objects are created the first time a combination appears and cached, keyed by the exact dimensions plus, for descriptors and algorithms, the transposes and the epilogue. One instance therefore serves GEMM calls at sizes that vary at runtime: a size seen for the first time creates and caches its objects, and a repeated size reuses them without another heuristic query.
+
+### addOperationConfig
+
+`addOperationConfig(m, n, k, lda, ldb, ldc, transa, transb, epilogue)` creates all three objects for one operation (the matrix layouts, the matmul descriptor and the algorithm) for the given dimensions, transposes and epilogue, before the corresponding call is made. It is optional: a combination that was never configured is created and cached on its first call. The `epilogue` argument is the `Epilogue` enum from `sofieBLAS/core.hpp` and names which call the site will make, because the fused epilogue is part of the selected kernel:
+
+| `Epilogue` value | call it configures |
+| --- | --- |
+| `Epilogue::Default` | `matmul` (no bias) |
+| `Epilogue::Bias` | `gemm` (adds the bias vector) |
+| `Epilogue::ReluBias` | `gemmrelu` (bias, then ReLU) |
+| `Epilogue::GeluBias` | `gemmgelu` (bias, then GELU) |
+
+### Initializing the cache limit
+
+The algorithm cache is unbounded by default. Passing a limit as the second constructor argument caps the number of cached algorithms; when an insertion would exceed the limit, the least recently used entries are evicted. Choose a limit at least as large as the number of distinct shapes the workload uses regularly, or leave it unbounded. `algoCacheSize()` returns the current number of entries.
+
+```cpp
+sofieBLAS<alpaka::TagGpuCudaRt> blas(queue);       // unbounded algorithm cache (default)
+sofieBLAS<alpaka::TagGpuCudaRt> capped(queue, 32); // at most 32 entries, LRU eviction
+
+blas.addOperationConfig(64, 3, 5, 64, 5, 64, 'N', 'N', Epilogue::Default);
+blas.matmul('N', 'N', 64, 3, 5, 1.0f, dA, dB, 0.0f, dC); // created by addOperationConfig: cache hit
+blas.matmul('N', 'N', 37, 3, 5, 1.0f, dA, dB, 0.0f, dC); // new size: created on first use
+blas.gemmrelu('N', 'N', 64, 3, 5, 1.0f, dA, dB, 0.0f, dBias, dC); // same size, other epilogue: layouts reused, descriptor and algorithm created on first use
+```
 
 
 ## Contributing

benchmark/bench.cc

@@ -139,7 +139,7 @@ static void runCudaBench(const BenchOptions &opt) {
     alpaka::memcpy(queue, dB, hB);
     alpaka::wait(queue);
 
-    blas.addLayoutConfig(M, N, K, M, K, M, 'N', 'N');
+    blas.addOperationConfig(M, N, K, M, K, M, 'N', 'N', Epilogue::Default);
 
     for (int i = 0; i < opt.warmup; ++i)
       blas.matmul('N', 'N', M, N, K, 1.f, dA, dB, 0.f, dC);
@@ -187,7 +187,7 @@ static void runHipBench(const BenchOptions &opt) {
     alpaka::memcpy(queue, dB, hB);
     alpaka::wait(queue);
 
-    blas.addLayoutConfig(M, N, K, M, K, M, 'N', 'N');
+    blas.addOperationConfig(M, N, K, M, K, M, 'N', 'N', Epilogue::Default);
 
     for (int i = 0; i < opt.warmup; ++i)
       blas.matmul('N', 'N', M, N, K, 1.f, dA, dB, 0.f, dC);