Large-data memory layout for equation systems

Zero-runtime-cost data-layout wrappers with type-based compile-time dispatch

Tags for selecting particle-, or property-major data layout policies:

const bool ParticleMajor = true;
const bool PropertyMajor = false;

Essentially, the implementation is as follows:

template< bool Major >
class Data {

  private:
   // Transform a compile-time bool into a type
   template< bool m >
   struct int2type {
     enum { value = m };
   };

   // Overloads for particle-, and property-major accesses
   inline
   tk::real& access( int particle, int property, int2type<ParticleMajor> ) {
     return *(m_ptr + particle*m_nprop + property);
   }
   inline
   tk::real& access( int particle, int property, int2type<PropertyMajor> ) {
     // This is the same for now, not called, irrelevant in zero-cost-test
     return *(m_ptr + particle*m_nprop + property);
   }

   tk::real* const m_ptr;
   const int m_nprop;

  public:
    // Constructor
    Data( tk::real* const ptr, int nprop ) :
      m_ptr(ptr), m_nprop(nprop) {}

    // Access dispatch
    inline tk::real& operator()( int particle, int property ) {
      return access( particle, property, int2type<Major>() );
    }
};

Test of zero-cost

Test by adding to Dirichlet (derived equation class) constructor (client code):

Data< Layout > d( particles, m_nprop );
Model::aa = d( 34, 3 );
Model::bb = *(m_particles + 34*m_nprop + 3);

Add to Model:

Model {
  ...
  public:
    tk::real aa;
    tk::real bb;
  ...
}

Add to Physics constructor:

std::cout << m_mix->aa << m_mix->bb;

This is so the optimizing compiler cannot entirely optimize the assignments of aa and bb away.

Debug assembly

Generated assembly code with CMAKE_BUILD_TYPE=DEBUG (i.e., without optimization) of the assignments of aa (line 42) and bb (line 43) in Dirichlet's constructor, with clang -g -S -mllvm --x86-asm-syntax=intel, gnu and intel generate very similar code:

     .loc    143 42 20  ; <quinoa>/src/DiffEq/Dirichlet.h:42:20
.Ltmp27038:
     lea     RDI, QWORD PTR [RBP - 56]
     mov     ESI, 34
     mov     EDX, 3
     call    _ZN6quinoa4DataILb1EEclEii
.Ltmp27039:
     mov     QWORD PTR [RBP - 176], RAX # 8-byte Spill
     jmp     .LBB2550_7
.LBB2550_7:
     mov     RAX, QWORD PTR [RBP - 176] # 8-byte Reload
     movsd   XMM0, QWORD PTR [RAX]
     mov     RCX, QWORD PTR [RBP - 64] # 8-byte Reload
     movsd   QWORD PTR [RCX + 8], XMM0
     .loc    143 43 0    ; <quinoa>/src/DiffEq/Dirichlet.h:43:0
     mov     RDX, QWORD PTR [RCX + 32]
     imul    ESI, DWORD PTR [RCX + 48], 34
     movsxd  RDI, ESI
     shl     RDI, 3
     add     RDX, RDI
     movsxd  RDI, DWORD PTR [RCX + 52]
     movsd   XMM0, QWORD PTR [RDX + 8*RDI + 24]
     movsd   QWORD PTR [RCX + 16], XMM0

Line 42 translates to register loads and a function call into tk::Data, while line 43 translates to some integer arithmetic of the address and loads.

Excerpt from the Intel® 64 and IA-32 Architectures Software Developer Manual

The LEA (load effective address) instruction computes the effective address in memory (offset within a segment) of a source operand and places it in a general-purpose register. This instruction can interpret any of the processor’s addressing modes and can perform any indexing or scaling that may be needed. It is especially useful for initializing the ESI or EDI registers before the execution of string instructions or for initializing the EBX register before an XLAT instruction.

The MOVSXD instruction operates on 64-bit data. It sign-extends a 32-bit value to 64 bits. This instruction is not encodable in non-64-bit modes.

A common type of operation on packed integers is the conversion by zero- or sign-extension of packed integers into wider data types. SSE4.1 adds 12 instructions that convert from a smaller packed integer type to a larger integer type (PMOVSXBW, PMOVZXBW, PMOVSXBD, PMOVZXBD, PMOVSXWD, PMOVZXWD, PMOVSXBQ, PMOVZXBQ, PMOVSXWQ, PMOVZXWQ, PMOVSXDQ, PMOVZXDQ). The source operand is from either an XMM register or memory; the destination is an XMM register.

IMUL Signed multiply. The IMUL instruction multiplies two signed integer operands. The result is computed to twice the size of the source operands; however, in some cases the result is truncated to the size of the source operands.

SAL/SHL Shift arithmetic left/Shift logical left. The SAL (shift arithmetic left), SHL (shift logical left), SAR (shift arithmetic right), SHR (shift logical right) instructions perform an arithmetic or logical shift of the bits in a byte, word, or doubleword. The SAL and SHL instructions perform the same operation. They shift the source operand left by from 1 to 31 bit positions. Empty bit positions are cleared. The CF flag is loaded with the last bit shifted out of the operand.

Optimized assembly

Generated assembly code with CMAKE_BUILD_TYPE=RELWITHDEBINFO (i.e., with optimization) of the assignments of aa (line 42) and bb (line 43) in Dirichlet's constructor, with clang -O2 -g DNDEBUG -S -mllvm --x86-asm-syntax=intel, gnu and intel generate very similar optimized code:

.loc    144 42 20  ; <quinoa>/src/DiffEq/Dirichlet.h:42:20
movsd   XMM0, QWORD PTR [R14 + 8*RAX + 24]
movsd   QWORD PTR [R13 + 8], XMM0
.loc    144 43 0    ; <quinoa>/src/DiffEq/Dirichlet.h:43:0
mov     RCX, QWORD PTR [R13 + 32]
movsd   XMM0, QWORD PTR [RCX + 8*RAX + 24]
movsd   QWORD PTR [R13 + 16], XMM0

Both lines 42 and 43 translate to very similar SSE loads with pointer arithmetic, i.e., line 42 costs the same as line 43.

Excerpt from the Intel® 64 and IA-32 Architectures Software Developer Manual

The MOVS instruction moves the string element addressed by the ESI register to the location addressed by the EDI register. The assembler recognizes three “short forms” of this instruction, which specify the size of the string to be moved: MOVSB (move byte string), MOVSW (move word string), and MOVSD (move doubleword string).

The MOVSD (move scalar double-precision floating-point) instruction transfers a 64-bit double-precision floating- point operand from memory to the low quadword of an XMM register or vice versa, or between XMM registers. Alignment of the memory address is not required, unless alignment checking is enabled.

Quinoa::Walker input file used for the benchmark

Note that walker is no longer part of Quinoa. Walker's code can be pulled from Quinoa's git history.

We will integrate for the duration of a 100,000 time steps a system of 100 coupled non-linear stochastic differential equations (SDEs) whose statistically stationary solution converge to the Dirichlet distribution and measure the wall-clock time. For more on the Dirichlet SDE, see src/DiffEq/Dirichlet.h.

walker

  nstep 100000  # Max number of time steps
  term  140.0   # Max time
  dt    0.05    # Time step size
  npar  40000   # Number of particles

  ttyi  100     # TTY output interval

  rngs
    mkl_mrg32k3a seed 0 end
  end

  dirichlet
    ncomp 100  # = K = N-1
    b     0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
          0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5 0.1 1.5
    end
    S     0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
          0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4 0.625 0.4
    end
    kappa 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
          0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3 0.0125 0.3
    end
    rng mkl_mrg32k3a
  end

  statistics    # Estimate statistics
    <Y1>        # mean of Y1
    <Y2>
    <y1y1>      # variance of Y1 = <(Y1-<Y1>)^2> = <y1^2>
    <y2y2>
    <y1y2>
  end

end

Ptr - Working with raw pointers

This algorithm gets the starting raw pointer from which the given particle data is (contiguously) accessible in memory and simply adds integers to the address to access and update the 100 components specified above. The algorithm assumes a particular data layout – it only works with the particle-major storage – a logically 3-dimensional array with [ particle ] [ equation ] [ component ] layout.

Layout-dependent algorithm:

//! Advance particles
void advance(int p, int tid, tk::real dt) override {
  // Get access to particle scalars
  tk::real* y = m_particles.ptr() + p*m_nprop;
  
  // Compute Nth scalar
  tk::real yn = 1.0 - y[0];
  for (int i=1; i<m_ncomp; ++i) yn -= y[i];
  
  // Generate Gaussian random numbers with zero mean and unit variance
  tk::real dW[m_ncomp];
  m_rng->gaussian( tid, m_ncomp, dW );
  
  // Advance first m_ncomp (K=N-1) scalars
  for (int i=0; i<m_ncomp; ++i) {
    tk::real d = m_k[i]*y[i]*yn*dt;
    if (d > 0.0) d = sqrt(d); else d = 0.0;
    y[i] += 0.5*m_b[i]*(m_S[i]*yn - (1.0-m_S[i])*y[i])*dt + d*dW[i];
  }
}

Par - Access via particle-major layout policy

This algorithm accesses particle data via the wrapper class, tk::Data, in a data-layout-agnostic fashion. Access itself via this class is demonstrably "zero-cost", i.e., an optimizing compiler completely optimizes the abstraction away: see Data layout implementation and assembly for the assembly generated by 3 compilers. However, writing an SDE-advance algorithm in a data-layout-agnostic manner, requires index calculations at every particle-access compared to working with raw pointers, as described above. Thus the following tests are designed to measure only the additional index calculations that the layout-agnostic access entails.

Layout-independent algorithm:

//! Advance particles
void advance(int p, int tid, tk::real dt) override {
  // Compute Nth scalar
  tk::real yn = 1.0 - m_particles(p, 0);
  for (int i=1; i<m_ncomp; ++i) yn -= m_particles(p, i);

  // Generate Gaussian random numbers with zero mean and unit variance
  tk::real dW[m_ncomp];
  m_rng->gaussian( tid, m_ncomp, dW );

  // Advance first m_ncomp (K=N-1) scalars
  for (int i=0; i<m_ncomp; ++i) {
    tk::real d = m_k[i] * m_particles(p, i) * yn * dt;
    if (d > 0.0) d = sqrt(d); else d = 0.0;
    m_particles(p, i) +=
      0.5*m_b[i]*(m_S[i]*yn - (1.0-m_S[i]) * m_particles(p, i) )*dt
      + d*dW[i];
  }
}

Comparison of the algorithms

DEBUG mode uses -O0 and does not optimize function calls away for all of three compilers tested. RELEASE mode uses -O3 and the abstraction is completely optimized away. However, index calculations still remain compared to a layout-dependent advance algorithm.

Total time measured in micro-seconds, run on a Lenovo laptop with Intel Core i7, 8 compute cores:

Data layout benchmark discussion

• As expected, inlining has a significant effect on performance: going from DEBUG to RELEASE mode yields a significant speedup with all three compilers, see last, DEBUGcode{.cmake}RELEASE, rows.
• As expected, the additional index calculations required by layout-agnostic access do take a performance hit: though only 6% with clang, and 3% with gnu, see last, Par/Ptr, columns.
• Surprisingly, the layout-agnostic access is even a tiny bit faster than the layout-dependent algorithm with the Intel compiler with -O3.