3 ms·

If you are interested on the exact same topics (matrix multiplication parallelization) here are other step by step tutorials I used: - UlmBLAS [1], for a HPC c

• gpderetta 7y ago
  Thanks for the pointers, I'll take a look. re fast-math, I was interested on how much additional parallelism the compile can extract from the naive code. Fast-math should,at least in principle, allow the compiler to unroll the loop and add the additional accumulators, although that violates strict IEEE semantics.
  • mratsim 7y ago
    Yes exactly. Usually you can extract 2x to 4x instruction level parallelism on simple add instructions [1] vs [2] For fused-multiply-add even though the latency is higher the instruction is also slower so its the same. BLIS, OpenBLAS and my code extract 2x parallelism (2 accumulators) at the lowest level because we are restricted by the number of registers and the fact that x86 can only issue 2 FMAs per cycle [3]. Details, we divided the work until we have a small micro matrix C or size MR x NR (input matrix of size MxK and KxN so output of MxN), and then here is are the constraints you have to deal with: Registers constraints and micro-kernel tuning - To issue 2xFMAs in parallel we need to use 2x SIMD registers - We want to hold C of size MR * NR completely in SIMD registers as well as each value is reused k times during accumulation C[i, j] += A[i, k] * B[k, j] - We should have enough SIMD registers left to hold the corresponding sections of A and B (at least 4, 2xA and 2xB for FMAs) On x86-64 X SIMD registers that can issue 2xFMAs per cycle: - NbVecs is 2 minimum - RegsPerVec = 2 * NbVecs => 4 minimum (for A and for B) - NR = NbVecs * NbScalarsPerSIMD - C: MR*NR and uses MR*NbVecs SIMD registers - MR*NbVecs + RegsPerVec <= X -> MR*NbVecs + 2 * NbVecs <= X -> (MR+2) * NbVecs <= X Some solutions: - AVX with 16 registers: - MR = 6, NbVecs = 2 FP32: 8xFP32 per SIMD --> NR = 2x8 ukernel = 6x16 FP64, ukernel = 6x8 - MR = 2, NbVecs = 4 FP32: 8xFP32 per SIMD --> NR = 4x8 ukernel = 2x32 FP64, ukernel = 2x16 - AVX512 with 32 registers - MR = 6, NbVecs = 4 FP32 ukernel = 6x64 FP64 ukernel = 6x32 - MR = 2, NbVecs = 8 FP32 ukernel = 2x128 FP64 ukernel = 2x64 - MR = 14, NbVecs = 2 FP32 ukernel = 14x32 FP64 ukernel = 14x16 And in-depth overview of the lowest level details is available in the paper Automating the last mile for High Performance Dense Linear Algebra[5]. In short, the compiler is completely unable to deal with this, and a high performance computing compiler should give an escape hatch to allow hand optimization like Halide[6] or Tiramisu do[7]. [1]: https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1187864323d85527d18c/benchmarks/fp_reduction_latency/reduction_bench.nim#L304-L341 https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1... [2]: https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1187864323d85527d18c/benchmarks/fp_reduction_latency/reduction_bench.nim#L493-L505 https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1... [3]: https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1187864323d85527d18c/laser/primitives/matrix_multiplication/gemm_tiling.nim#L111-L146 https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1... [4]: https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1187864323d85527d18c/laser/primitives/matrix_multiplication/gemm_ukernel_generator.nim#L149-L160 https://github.com/numforge/laser/blob/e660eeeb723426e80a7b1... [5]: https://arxiv.org/pdf/1611.08035.pdf https://arxiv.org/pdf/1611.08035.pdf [6]: https://halide-lang.org/ https://halide-lang.org/ [7]: http://tiramisu-compiler.org/ http://tiramisu-compiler.org/