File DenseMatrix.h

↰ Parent directory (ADS/numeric)

The engineered dense matrices, beyond the book: DenseMatrix keeps the m x n entries by rows in one contiguous vector, instead of the vector of rows of File Matrix.h [see File Matrix.h], and the product, the transposition and the solution of linear systems run on it at the speed of the processor.

Definition (ADS/numeric/DenseMatrix.h)

Detailed Description

  • gemm() is :math:` C leftarrow alpha A B + beta C ` on matrices by rows with leading dimensions, organized as the GotoBLAS and BLIS libraries do: the loop on the panels of NC columns of B, the one on the slices of KC of the common dimension, whose panel of B is packed in micro-panels of NR columns (KC x NC in the last level cache), the one on the blocks of MC rows of A, packed in micro-panels of MR rows (MC x KC in the second level), and two loops of micro-tiles MR x NR of C, each the sum of KC products of a column of MR entries of A by a row of NR entries of B, held in registers: MR NR / w vector accumulators, w the entries of a vector, each step loading NR / w vectors of B and broadcasting MR entries of A into fused products and sums. The micro-kernel is written with the intrinsics of AVX2 for double (6 x 16, 24 accumulators of 4) and float, and in portable C++ otherwise, where the compiler vectorizes what it can; ADS_GEMM_AVX512 selects a kernel of AVX-512 (12 x 16, 24 accumulators of 8) for the processors with two units of fused products of 512 bits, since on those with one (e.g., the Xeon Gold 5120) both widths make 16 operations a cycle and the vectors of 256 bits run at a higher clock. KC = 384, MC = 480 and NC = 4096 came from a search on the same processor. The product makes :math:` 2 m n k ` operations in :math:` Theta( m n k ) ` time; a step of the micro-kernel makes :math:` M_R N_R ` fused operations on :math:` M_R + N_R ` entries loaded from the first level ( :math:` 96 / 22 ` for 6 x 16), the packing moves :math:` m k lceil n / N_C rceil + k n ` entries, a fraction :math:` O( 1 / M_C + 1 / N_C ) ` of the work, and the buffers take :math:` K_C ( M_C + N_C ) ` entries, kept by the thread between the calls. The partial tiles at the borders are packed with zeros, and only their entries inside C are written. A small product ( :math:` m n k le 2^{21} `) packs B alone and reads A where it is, the tiles of its last rows and columns by kernels of their shape.

  • transpose() moves blocks of 32 x 32 entries, which stay in the cache of first level both read by columns and written by rows, instead of running through whole columns: :math:` Theta( m n ) ` time and :math:` O( m n / L ) ` misses of lines of L entries when a block of :math:` 32 L ` entries fits the cache (the tall cache), the bound of the cache-oblivious transposition.

  • multiply() is the product A B, a new matrix, by gemm().

  • multiplyStrassen() is the product by the variant of Winograd of the algorithm of Strassen [see strassen() of File Matrix.h], on the views of the quadrants: with :math:` S_1 = A_{21} + A_{22} , :math: S_2 = S_1 - A_{11} , :math: S_3 = A_{11} - A_{21} , :math: S_4 = A_{12} - S_2 , :math: T_1 = B_{12} - B_{11} , :math: T_2 = B_{22} - T_1 , :math: T_3 = B_{22} - B_{12} , :math: T_4 = T_2 - B_{21} ` and the products :math:` P_1 = A_{11} B_{11} , :math: P_2 = A_{12} B_{21} , :math: P_3 = S_4 B_{22} , :math: P_4 = A_{22} T_4 , :math: P_5 = S_1 T_1 , :math: P_6 = S_2 T_2 , :math: P_7 = S_3 T_3 , the quadrants of C are :math: C_{11} = P_1 + P_2 , :math: C_{12} = U + P_5 + P_3 , :math: C_{21} = U + P_7 - P_4 , :math: C_{22} = U + P_7 + P_5 ` with :math:` U = P_1 + P_6 : 7 products and 15 sums of halves, :math: T( n ) = 7 T( n / 2 ) + O( n^2 ) = O( n^{log_2 7} ) = O( n^{2.81} ) , until a dimension falls below the cutoff, where gemm() takes over; an odd dimension leaves its last row or column to gemm() (dynamic peeling). The bound of the error grows as :math: ( n / n_0 )^{log_2 18} n_0^2 u | A | | B | ` for the cutoff :math:` n_0 ` and the unit roundoff u, against the :math:` n u | A | | B | ` of the classical product (Higham), which is the price of the fewer products: with the cutoff 3072, found on the benchmark, one level spares 1/8 of the operations of a product of 4096; the schedule of the sums needs three temporaries of a quadrant.

The contracts [see File Contract.h]: the sizes of the factors agree, the entries read are inside the matrices; gemm() needs leading dimensions at least the columns. The functions are defined in DenseMatrix.tpp.

  1. Goto, R. A. van de Geijn, Anatomy of high-performance matrix multiplication, ACM Transactions on Mathematical Software 34(3), 12:1-25, 2008, doi:10.1145/1356052.1356053

    1. Van Zee, R. A. van de Geijn, BLIS: a framework for rapidly instantiating BLAS functionality, ACM Transactions on Mathematical Software 41(3), 14:1-33, 2015, doi:10.1145/2764454

    1. Low, F. D. Igual, T. M. Smith, E. S. Quintana-Orti, Analytical modeling is enough for high-performance BLIS, ACM Transactions on Mathematical Software 43(2), 12:1-18, 2016, doi:10.1145/2925987

  1. Strassen, Gaussian elimination is not optimal, Numerische Mathematik 13(4), 354-356, 1969, doi:10.1007/BF02165411

  1. Winograd, On multiplication of 2 x 2 matrices, Linear Algebra and its Applications 4(4), 381-388, 1971, doi:10.1016/0024-3795(71)90009-7

    1. Higham, Exploiting fast matrix multiplication within the level 3 BLAS, ACM Transactions on Mathematical Software 16(4), 352-368, 1990, doi:10.1145/98267.98290

  1. Huang, T. M. Smith, G. M. Henry, R. A. van de Geijn, Strassen’s algorithm reloaded, SC16: International Conference for High Performance Computing, Networking, Storage and Analysis, 690-701, 2016, doi:10.1109/SC.2016.58

  1. Frigo, C. E. Leiserson, H. Prokop, S. Ramachandran, Cache-oblivious algorithms, ACM Transactions on Algorithms 8(1), 4:1-22, 2012, doi:10.1145/2071379.2071383

Donato Meoli

Dipartimento di Informatica Universita’ di Pisa

by Donato Meoli

Includes

Included By

Namespaces

Classes

Functions