Improve Performance with Vectorization

This article focuses on the steps to improve software performance with vectorization. Included are examples of full applications along with some simpler cases to illustrate the steps to vectorization. As hardware moves forward adding more cores and wider vector registers, software must modernize or change to match the hardware and deliver higher levels of parallelism and vectorization. This increased performance will solve more complex problems in finer resolution. A previous article, Recognizing and Measuring Vectorization Performance, explained how to measure how effectively the software was vectorized.

In this article the terms SIMD and vectorization are used synonymously. A brief review of vectorization is discussed in the introduction. SIMD stands for Single Instruction Multiple Data; the same instruction is used on multiple data elements simultaneously. Intel® Advanced Vector Extensions 512 (Intel® AVX-512) provides registers that are 512 bits wide. One of these registers may be filled with 8 double precision floating point values, or 16 single precision floating point values, or 16 integers. When the register is fully populated a single instruction is applied that operates on 8 to 16 values (depending on the types of values and instruction used) simultaneously (see Figure 1).

Figure 1: Example of SIMD operations on zmm Intel® Advanced Vector Extensions 512 registers.

When the Intel AVX-512 registers are filled across all available lanes performance can be 8 to 16 times faster. Typically there are many scalar operations to perform, as well as some shuffling and data movement, so a speedup of 8 to 16 is not delivered for a full application (though performance increases significantly). Intel® Advanced Vector Extensions 2 (Intel® AVX2) and Intel® Advanced Vector Extensions (Intel® AVX) instructions are half the width of Intel AVX-512 instructions and Intel® Streaming SIMD Extensions 4 instructions are half the width of Intel AVX and Intel AVX2 instructions. All deliver increasing levels of performance.

Taking software to new performance levels with vectorization

How does a software developer go from scalar code to vector mode? There are several steps to follow:

1. Use vectorized libraries (for example, the Intel® Math Kernel Library)
2. Think SIMD (think what is done repetitively or multiple times)
3. Do a performance analysis (collect data, and then improve)
4. Express SIMD

The first three steps are covered in the article referenced above. This article focuses on expressing SIMD and code modifications and includes the following:

• Use pragmas and directives
• Write SIMD-enabled functions
• Remove loop invariant dependencies
• Expand temporary variables to short vectors
• Pay attention to data alignment
• Improve data layout

Pragmas/Directives and dependence

Compilers do a great job of vectorization. There are still times when a compiler is not able to disambiguate memory references inside a repeated do loop to safely guarantee that there are no dependencies between one loop iteration and another. When the compiler cannot determine data independence it will produce scalar code instead of vectorized code. Other times the analysis for the compiler may become too complex, and the compiler may stop and produce scalar code rather than vectorized code.

The most common method to improve vectorization is the addition of pragmas or directives. A compiler pragma (C/C++) or directive (Fortran*) gives the compiler more information about the code segment than it can discern from the code alone allowing the compiler to perform more optimizations. The most common problem for which pragmas or directives are applied is pointer disambiguation, which is the inability to determine if it is safe to vectorize a for or do loop. If you know that the data entering the do (or for) loop does not overlap, you can add a pragma or a directive instructing the compiler to ignore the potential dependency and vectorize it.

Various pragma/directives were developed by different compiler groups. The need to tell compilers to ignore potential data dependencies is so common that it is easier for software developers to have a common set of pragmas and directives recognized by all compilers rather than proprietary pragmas for each compiler. The OpenMP* organization filled that void by adding a SIMD compiler pragma/directive to the OpenMP 4.0 specifications in 2013. This way you can use a standard directive/pragma and know that it will be recognized by many compilers.

The syntax for this directive/pragma is:

#pragma omp simd
!$omp simd

There are additional options to the SIMD pragma that may also be useful:

safelen(length),linear(list:linear-step),aligned(list[alignment]) , . . .   

These additional parameters may be useful if the code contains something like:

for (i=8;i<n;i++)
{    . . .
  x[i] += x[i-8]* . . . ;  . . .
}

There is a dependence between iterations. However, there is a safe length for this (the dependence only goes back 8 values) so the compiler can safely vectorize for a SIMD length of 8. In this case, use a pragma like:

#pragma omp simd safelen(8)

When you use the SIMD pragma/directive, you’re telling the compiler that you know enough about the program and data access that the assumption regarding no dependencies is safe. If you’re wrong, the code may deliver incorrect results. There is a quick check that some developers use for a sanity check—reverse the order of loops—but this is not a guarantee. If reversing the order of the loop changes the value, the code is not safe to vectorize or thread without changes. Additional specifications such as safelen or an atomic operation or other controls may be required to make the code safe to vectorize.

If the results are correct, it might be that for this particular dataset the order of operations did not matter, but it could matter with a different dataset. So failing the check clearly shows dependence, passing the check is not a guarantee of safety. Tools such as Intel® Advisor XE provide a more rigorous analysis and can help identify dependencies and give you hints as to when it is safe to apply SIMD pragmas or directives.

Functions

Frequently a function or subroutine call within a loop prevents vectorization. There are a couple of possibilities that may help in this situation. A short function can be written as a vectorizable function. Vector versions of many transcendental functions exist in math libraries used by the compiler. With the proper syntax, the compiler will generate both vector and scalar versions of a user-written function.

Let’s say there is a function that squares a number and adds a constant to the value. The function might be like this:

double sqadd(double a, double r)
{
   double t ;
   t = a*a + r ;
   return(t) ;
}

To tell the compiler to create a vector version of this function as well as a scalar version write it like this:

#pragma omp declare simd // tell compiler to generate vector versiondouble sqadd(double a, double r)
{
   double t ;
   t = a*a + r ;
   return(t) ;
}

Next, at the places where the function is called, further instruct the compiler to vectorize the loop in the presence of a function call like this:

#pragma omp simd // instruct compiler to generate vector code even though loop invokes function
for (i=0 ; i<n; ++)
{
   . . .
  anarray[i] = sqadd(r, t) ;
  . . .
}

In Fortran the correct directive to use is: 

!$omp declare simd(subroutine name) 

For example the above function in Fortran may be written as shown below to generate vectorizable implementations:

     real*8 function sqadd(r,t) result(qreal)
!$omp declare simd(sqadd)
         real*8 :: r
         real*8 :: t
         qreal = r*r + t
       end function sqadd

In the file where the function is invoked I defined it like this:

      INTERFACE
         real*8 function sqadd(r,t)
!$omp declare simd(sqadd)
         real*8 r
         real*8 t
         end function sqadd
      END INTERFACE

The compiler vectorized the do loop containing the call to sqadd. If modules are used, and the function or subroutine is declared SIMD enabled, any file that uses that module can see the function/subroutines are SIMD enabled and then generate vectorized code. For more information on creating SIMD functions/subroutines in Fortran, see this article on explicit vector programming in Fortran.

Loop invariant dependencies

It is common to see do/for loops with an if/else statement in the middle of them where the condition of the if/else statement does not change within the loop. The presence of the conditional may prevent vectorization. This can be rewritten so that it is easier to vectorize. For example if the code looks something like:

for (int ii = 0 ; ii < n ; ++ii)
{
     . . .
     if (method == 0) 
        ts[ii] = . . . . ;
    else
         ss[ii] = . . . ;     
     . . . 
}

This can be rewritten as:

if (method == 0)
    for (int ii = 0; ii < n ; ++ii)
    {
        . . . 
       ts[ii] = . . . ;
       . . . 
    }
else
   for (int ii = 0 ; ii < n; ++ii)
   {
       . . . 
      ss[ii] = . . . ;
       . . . 
   }

The above two techniques were applied to the MPAS-Oⁱ oceanic code. The MPAS-O code simulates the earth's ocean system and can work on time scales of months to thousands of years. It handles regions below 1 km as well as global circulation. Doug Jacobsenⁱⁱ from LANL participated in a short collaboration with ParaTools where we used the TAU Performance System* to measure vectorization intensity of this routine on the Intel® Xeon Phi™ coprocessor. Based on the data we looked where the report indicated low vector intensity and then at the code and compiler report.

Excerpts of the code are shown below:

 do k=3,maxLevelCell(iCell)-1
    if(vert4thOrder) then
       high_order_vert_flux(k, iCell) = &
          mpas_tracer_advection_vflux4( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
          tracer_cur(k,iCell),tracer_cur(k+1,iCell), w(k,iCell))
    else if(vert3rdOrder) then
       high_order_vert_flux(k, iCell) = &       
          mpas_tracer_advection_vflux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
          tracer_cur(k,iCell),tracer_cur(k+1,iCell), w(k,iCell), coef_3rd_order )
    else if (vert2ndOrder) then
       verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) +&
          verticalCellSize(k-1, iCell))
       verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) +&
          verticalCellSize(k-1, iCell))
       high_order_vert_flux(k,iCell) = w(k,iCell) * (verticalWeightK * tracer_cur(k,iCell) +&
          verticalWeightKm1 * tracer_cur(k-1,iCell))
    end if
    tracer_max(k,iCell) = max(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
    tracer_min(k,iCell) = min(tracer_cur(k-1,iCell),tracer_cur(k,iCell),tracer_cur(k+1,iCell))
 end do

We quickly saw that this code contained both a subroutine call and invariant conditional within the do loop. So we made the subroutine vectorizable and swapped the positions of the loop and conditional. The new code looks like this:

 ! Example flipped loop
 if ( vert4thOrder ) then
    do k = 3, maxLevelCell(iCell) - 1
       high_order_vert_flux(k, iCell) = &
          mpas_tracer_advection_vflux4( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
          tracer_cur(k,iCell),tracer_cur(k+1,iCell), w(k,iCell)) 
    end do
 else if ( vert3rdOrder ) then
    do k = 3, maxLevelCell(iCell) - 1
       high_order_vert_flux(k, iCell) = &
          mpas_tracer_advection_vflux3( tracer_cur(k-2,iCell),tracer_cur(k-1,iCell), &
          tracer_cur(k,iCell),tracer_cur(k+1,iCell), w(k,iCell), coef_3rd_order )
    end do
 else if ( vert2ndOrder ) then
    do k = 3, maxLevelCell(iCell) - 1
        verticalWeightK = verticalCellSize(k-1, iCell) / (verticalCellSize(k, iCell) +&
           verticalCellSize(k-1, iCell))
        verticalWeightKm1 = verticalCellSize(k, iCell) / (verticalCellSize(k, iCell) +&
           verticalCellSize(k-1, iCell))
        high_order_vert_flux(k,iCell) = w(k,iCell) * (verticalWeightK * tracer_cur(k,iCell) +&
           verticalWeightKm1 * tracer_cur(k-1,iCell))
    end do
 end if

In this case we combined two of the techniques listed in this article to get vectorized code and improve performance: make subroutines/functions SIMD enabled or vectorizable and move out invariant conditionals. You will typically find that you must apply multiple techniques to tune your code. Each change is a step in the correct direction; sometimes the performance doesn't really jump ahead until all the changes are made; change data layout, apply pragmas or directives, align data and such. Each element is a step to improve performance. Several more techniques follow in the upcoming sections.

Expand temporary scalars into short arrays

Temporary values are commonly calculated in the middle of a for or do loop. Sometimes this is done because an intermediate value is used in several calculations, and the common subset of computations is calculated and then held for a short time for reuse. Sometimes it just makes the code easier to read. If you apply “think SIMD” to this situation it means to make the temporary arrays into short vectors that are as long as the widest SIMD register the code will run on.

This technique was applied to finite element and finite volume methods solving a system of partial differential equations related to calcium dynamics in a heart cell. Calcium ions play a crucial role in driving the heartbeat. Calcium is released into the heart cell at a lattice of discrete positions throughout the cell, known as calcium release units (CRUs). The probability that calcium will be released from a CRU depends on the calcium concentration at that CRU. When a CRU releases calcium (it “fires”), the local concentration of calcium ions increases sharply, and the calcium that diffuses raises the probability for release at neighboring sites. As a sequence of CRUs begins to release calcium throughout the cell, the release can self-organize into a wave of increasing concentration. A wave that is triggered among the normal physiological signaling (for example, the cardiac action potential) can lead to irregular heartbeats and possibly a life-threatening ventricular fibrillation. These dynamics are simulated by a system of three time-dependent partial differential equations developed by Leighton T. Izu.

Consider special-purpose MPI code that simulates these calcium dynamics developed by Matthias K. Gobbert and his collaborators (www.umbc.edu/~gobbert/calcium). This code solves the system of partial differential equations using finite element or finite volume methods and matrix-free linear solvers. Using a profiler such as the TAU Performance System, we find that most of the runtime is spent in the matrix vector multiplication function.

Here is a representative code snippetⁱⁱⁱ from this function (variables or functions changed are in bold italics):

struct car_t {
    void randomize();

    color_t color;
    model_t model;
    std::string regnr;
    double weight;
    double hp;
};

class database_2_vector_t : public single_question_database_t
{
public:
    virtual void add_one(const car_t &c);
    virtual void finalize_adding();
    virtual void ask_one_question(query_t &query) const;

private:
    typedef std::vector<car_t> cars_t;
    cars_t cars;
};

class database_3_hot_cold_vector_t : public single_question_database_t
{
public:
    virtual void add_one(const car_t &c);
    virtual void finalize_adding();
    virtual void ask_one_question(query_t &query) const;

private:
    typedef std::vector<color_t> colors_t;
    typedef std::vector<model_t> models_t;
    typedef std::vector<double> weights_t;

    colors_t colors;
    models_t models;
    weights_t weights;

    typedef std::vector<car_t> cars_t;
    cars_t cars;
};