Parallel Programming on a CPU with AVX-512

This text is the second of a two-part collection that presents two distinctly totally different approaches to parallel programming. Within the two articles, I take advantage of totally different approaches to unravel the identical downside: discovering the best-fitting line (or regression line) for a set of factors.

The 2 totally different approaches to parallel programming offered on this and the previous Insights article (Parallel Programming on an NVIDIA GPU | Physics Boards) use these applied sciences.

• Single-instruction multiple-thread (SIMT) programming is offered on the Nvidia® household of graphics processing items (GPUs). In SIMT programming, a single instruction is executed concurrently on a whole lot of microprocessors on a graphics card.
• Single-instruction a number of knowledge (SIMD) as offered on x64 processors from Intel® and AMD® (this text). In SIMD programming, a single instruction operates on vast registers that may comprise vectors of numbers concurrently.

On this article, I describe a program that makes use of Intel AVX-512 meeting directions and features a comparability of the outcomes from each packages.

Regression line calculations

Given a set of factors (X_i, Y_i), the place i ranges from 1 to N, the slope m and y-intercept b of the regression line will be discovered from the next formulation.

$$m = frac{N left(sum_{i = 1}^N X_iY_i proper)~ -~  left(sum_{i = 1}^N X_i proper) left(sum_{i = 1}^N Y_i proper)}{Nleft(sum_{i = 1}^N X_i^2right) – left(sum_{i = 1}^N X_iright)^2 }$$

$$b = frac{sum_{i = 1}^N Y_i – sum_{i = 1}^N X_i} N$$

As will be seen from these formulation, there are numerous calculations that should be carried out. This system should calculate the sum of the x coordinates, in addition to the sum of the y coordinates. It should additionally calculate the sum of the squares of the x coordinates and the sum of the term-by-term xy merchandise. For the latter two sums, it’s handy to create two new vectors. The primary vector consists of the term-by-term merchandise of the x coordinates with themselves. The second consists of the term-by-term product of the x- and y-coordinates.

Just like the GPU model within the previous article, this system I’m presenting right here makes use of parallel programming methods (by way of the CPU’s AVX-512 meeting directions) to calculate the term-by-term vector merchandise##< X_i Y_i>## and ##< X_i^2>## . In contrast to the GPU model of this system within the previous article, this program additionally makes use of AVX-512 meeting routines to calculate the 4 sums ##sum X_i##, ##sum Y_i##, ##sum X_iY_i##, and ##sum X_i^2##.

Foremost program duties

The primary program, written in C, performs the next duties:

• Declare exterior (i.e., assembly-coded) capabilities and declare the arrays and variables that shall be used.
• Initialize enter vectors.
• Name an meeting routine to calculate the term-by-term product vectors ##<X_iY_i>## and ##<X_i^2>##.
• Name a special meeting routine to calculate the 4 vector sums ##sum X_i##, ##sum Y_i##, ##sum X_iY_i##, and ##sum X_i^2##.
• Calculate the slope m and y-intercept b of the regression line.
• Show the outcomes of the calculations.

Every of those steps is described within the following sections.

Operate prototypes and variable declarations

The 2 prototypes under inform the compiler that the definitions of sumVector and vecProduct seem elsewhere within the undertaking (in a separate file). These two capabilities are written in pure meeting code, with a lot of it in AVX-512 meeting code.

The 4 arrays declared under symbolize the vectors ##<X_i>, <Y_i>, <X_i^2>,## and ##<X_iY_i>##. The primary two vectors are, respectively, the x- and y-coordinates of 262,144 factors. The 4 arrays are aligned on 64-byte boundaries, as required by the related AVX-512 directions.

// Prototypes for meeting routines.
extern "C" double sumVector(double[], int);
extern "C" void vecProduct(double output[], double arr1[], double arr2[], int rely);

const int NUMELEMENTS = 512 * 512;
// Allocate the enter vectors X, Y, XY, XX
__declspec(align(64)) double X[NUMELEMENTS];
__declspec(align(64)) double Y[NUMELEMENTS];
__declspec(align(64)) double XX[NUMELEMENTS];
__declspec(align(64)) double XY[NUMELEMENTS];

Enter vector initialization

To make issues easy, this system generates factors that lie on a line. Doing issues this manner makes it easy to verify that the calculated slope and y-intercept of the regression line are appropriate. Within the code under that initializes the 2 coordinate vectors, I’ve “unrolled” the loop by calculating 4 units of factors per loop iteration, fairly than doing only one pair of coordinates at a time.

// Initialize the enter vectors.
// All factors lie on the road y = 1*x + 0.5.
const double slope = 1.0;
const double y_int = 0.5;

for (int i = 0; i < NUMELEMENTS; i += 4)
{
    X[i] = slope * i;
    Y[i] = X[i] + y_int;
    X[i+1] = slope * (i+1);
    Y[i+1] = X[i+1] + y_int;
    X[i+2] = slope * (i+2);
    Y[i+2] = X[i+2] + y_int;
    X[i+3] = slope * (i+3);
    Y[i+3] = X[i+3] + y_int; 
}

Calls to vecProduct

At this level the 2 enter vectors have been initialized with 262,144 (= 512 X 512) factors. After vecProduct returns within the first line, XY will comprise the term-by-term merchandise of the corresponding components of the X and Y vectors.

The second name to vecProduct will set the vector XX equally.

// Calculate vector product X*Y. 
vecProduct(XY, X, Y, NUMELEMENTS);

// Calculate vector product X*X.
vecProduct(XX, X, X, NUMELEMENTS);

Calls to sumVector

The calculation of the regression line parameters additionally requires 4 vector sums: ##sum X_i##, ##sum Y_i##, ##sum X_iY_i##, and ##sum X_i^2##.

double sumX = 0;
sumX = sumVector(X, NUMELEMENTS);

The opposite three sums are calculated similarly.

Slope and intercept calculations

In any case 4 vector sums are calculated, it’s a easy matter of plugging them into the formulation for the slope (m) and y-intercept (b) of the regression line.

m = (double)(NUMELEMENTS * sumXY - sumX * sumY) / (NUMELEMENTS * sumXX - sumX * sumX);
b = (double)(sumY - sumX) / NUMELEMENTS;

Show the outcomes

Reasonably than exhibiting the code that shows the outcomes, I’ll simply present this system output. The final 4 strains are a sanity examine: a comparability of the computed values for the slope and y-intercept in opposition to the know values of those line parameters.

[Linear regression of 262144 points]
Sum of x: 34359607296.000000
Sum of y: 34359738368.000000
Sum of xy: 6004782323269632.000000
Sum of x^2: 6004765143465984.000000
Processed 262144 factors
Predicted worth of m: 1.000000
Computed worth of m: 1.0000000000
Predicted worth of b: 0.500000
Computed worth of b: 0.5000000000

Meeting code

A word on AVX-512 registers

An AVX-512 register resembling ZMM0 is 64 bytes or 512 bits vast and might comprise a vector of 8 double values, 16 float or int values, 32 brief values, or 64 char values. Its decrease half of ZMM0, which incorporates 32 bytes, has an alias of YMM0. In distinction, the higher half of any ZMM register doesn’t have an alias. Equally, the decrease half of YMM0 has an alias of XMM0, however the higher half of any YMM register doesn’t have an alias. Determine 1 reveals the relationships between ZMM, YMM, and XMM registers.

Fig. 1 ZMM, YMM, and XMM registers

Information part

The info part defines storage for STRIDE. That is used to increment the supply and vacation spot index registers RSI and RDI by the variety of bytes in 8 values of sort double (64 bytes).

.knowledge

STRIDE dq 64

vecProduct routine

The vecProduct routine takes two vectors as enter and shops the term-by-term merchandise in a 3rd vector.

vecProduct PROC C
    ;-----------------------------------------------------------------------------------------------
    ; C prototype: void vecProduct(double * outArr, double * inArr1, double * inArr2, int array_len)
    ; The vecProduct proc calculates the term-by-term merchandise of the weather
    ; of inArr1 and inArr2, and shops the merchandise in outArr. 
    ; 
    ; Enter args: 
    ;     outArr - tackle of output array, handed in RCX.
    ;     inArr1 - tackle of first enter array, handed in RDX.
    ;     inArr2 - tackle of second enter array, handed in R8.
    ;     arr_len - variety of components of every array, handed in R9.
    ; Return worth: None.
    ; On return, RAX holds the tackle of the vector product.
    ; In every iteration, 8 doubles are learn from reminiscence and saved in ZMM1, 
    ;   and the following 8 doubles are learn from reminiscence and saved in ZMM2.
    ;-----------------------------------------------------------------------------------------------

; Prologue
    push rdi         ; RDI is a non-volatile register, so put it aside.
    sub rsp, 20h 
    mov r10, rsp
    push rsi          ; Save RSI as properly.

; Initialization 
    vzeroall          ; Zero out all ZMM registers.
    mov rax, rcx      ; Copy output array tackle to RAX.
    shr r9, 3         ; r9 <- Variety of octets of doubles.
    mov rcx, r9       ; Copy no. of 16-float chunks to RCX
    xor rsi, rsi      ; Zero RSI. 
    xor rdi, rdi      ; Zero RDI. 
    vorpd xmm3, xmm3, xmmword ptr [HIMASK] ; Set bits in decrease half of XMM3 for later use.

ProcessChunks:
    ;-----------------------------------------------------------------------------------------------
    ; Iterate via each arrays, doing term-by-term multiplication.
    vmovapd zmm1, zmmword ptr [rdx + rsi]    ; Retailer 8 doubles from first array to ZMM1.
    vmovapd zmm2, zmmword ptr [r8 + rsi]     ; Retailer 8 doubles from second array to ZMM2.

    ; Do term-by-term multiplication of ZMM1 and ZMM2, storing end in ZMM1. 
    vmulpd zmm0, zmm1, zmm2                  ; ZMM0 = ZMM1 * ZMM2

    ; Retailer merchandise from earlier loop iteration into output array reminiscence.
    vmovapd zmmword ptr[rax + rdi], zmm0 
    add rsi, STRIDE           ; Advance 8 doubles (64 bytes) larger in first array.
    add rdi, STRIDE           ; Advance 8 doubles (64 bytes) larger in second array.
    loop ProcessChunks 
    ;-----------------------------------------------------------------------------------------------

; Epilogue
    pop rsi              ; Restore RSI.
    mov rsp, r10
    add rsp, 20h         ; Regulate the stack again to unique state.
    pop rdi              ; Restore RDI. 
    ret 
vecProduct ENDP

How vecProduct works

Of the 2 meeting routines, vecProduct is rather a lot less complicated to elucidate. The primary and final sections, labelled Prologue and Epilogue, are boilerplate sections of code which might be required in 64-bit meeting routines. Suffice it to say that they simply put aside some stack reminiscence after which later reclaim it.

Initialization part

The Initialization part zeroes the entire 512-bit ZMM registers, and copies the output and enter vector addresses from the registers used to move them into totally different registers. It additionally units up the RCX register to the variety of iterations of a loop that shall be used to course of the entire values within the enter vectors. The loop will multiply 8 values of sort double from every of the 2 enter vectors in every iteration, and retailer the outcome within the output vector. Because of this, RCX is about to the variety of double octets (64 bytes or 512 bits from every enter vector) that should be processed.

ProcessChunks part

Within the ProcessChunks part, the code copies (strikes) 8 double values to ZMM1, after which copies 8 extra to ZMM2. These directions use a pair of AVX-512 directions, vmovapd. (Be aware that the primary objective of this instruction is to move values from one place to a different.) Lots of the AVX-512 directions have a v prefix, which signifies that they work on vectors of numbers, versus single numbers. The apd suffix is an abbreviation for aligned (reminiscence), packed double.

The 8 double values in every of ZMM1 and ZMM2 are multiplied pairwise utilizing the vmulpd instruction (a multiply operation), with the outcomes being positioned in ZMM0. Once more, the pd suffix is brief for packed double. Determine 2 reveals the outcomes of the vmulpd instruction, multiplying every corresponding double within the two supply registers, with the outcomes being saved in ZMM0.

Fig. 2 Vector multiplication

The leads to ZMM0 are then copied to the reminiscence for the output vector utilizing vmovapd once more, the place the output vector’s base vacation spot tackle is in RAX.

After this, the index registers RSI and RDI are incremented by STRIDE (64 bytes, the dimensions of 8 doubles), and the loop instruction decrements RCX by 1. Management flows again to the primary instruction after the ProcessChunks label, offered that RCX continues to be nonzero. When RCX lastly will get to zero, management drops via to the following line after the loop instruction.

When vecProduct completes execution, RAX holds the tackle of the vector of element-by-element merchandise of the 2 enter vectors.

sumVector routine

The sumVector routine sums the entire components of its enter vector and returns this sum as a double worth.

sumVector PROC C
    ;----------------------------------------------------------------------------------------------- 
    ; C prototype: float sumVector(double * inputArray, int array_len)
    ; The sumVector proc provides the weather of a handed array. 
    ; Parameters
    ;     inputArray - tackle of an array of doubles, handed in RCX
    ;     array_len - variety of components of inputArray, handed in RDX
    ; Return worth
    ; Returns the sum of the weather within the handed array in XMM0.
    ; In every loop iteration, zmm1 is about with 8 doubles. 
    ;-----------------------------------------------------------------------------------------------

; Prologue
    push rdi ; RDI is a non-volatile register, so put it aside.
    sub rsp, 20h 
    mov rdi, rsp
    push rsi ; RSI additionally must be saved.

; Initialization
    vzeroall ; Zero out all ZMM registers
    mov rax, rcx ; Copy array tackle to RAX.
    mov rcx, rdx ; Copy factor rely to RCX.
    shr rcx, 4 ; RCX <- Variety of 2 X 8-double chunks.
    xor rsi, rsi ; Zero RSI. 
    
PartialSumsLoop:
    ;-----------------------------------------------------------------------------------------------
    ; In every iteration, accumulate 8 partial sums of doubles.
    ; When loop terminates, ZMM0 will maintain 8 doubles that, 
    ; when added, would be the general sum of the given array.
    vmovapd zmm1, zmmword ptr [rax + rsi] ; Get 8 doubles and retailer in ZMM1.
    vaddpd zmm0, zmm0, zmm1 ; ZMM0 = ZMM0 + ZMM1. 
    add rsi, STRIDE
    vmovapd zmm4, zmmword ptr [rax + rsi] ; Get 8 extra doubles and retailer in ZMM4. 
   
    ; Accumulate sum from earlier loop iteration.
    vaddpd zmm0, zmm0, zmm4 ; ZMM0 = ZMM0 + ZMM4.
    add rsi, STRIDE ; Increment RSI to advance to subsequent block of 8 doubles within the array.
    loop PartialSumsLoop 
    ;-----------------------------------------------------------------------------------------------

CombineSums: 
    ; Use vextractf64x4 to extract 4 doubles (256 bits) from the higher half of ZMM0 to YMM1.
    ; YMM0 already incorporates the decrease 256 bits (4 doubles) of ZMM0.
    vextractf64x4 ymm1, zmm0, 1 
    vaddpd ymm1, ymm1, ymm0 ; YMM0 <-- YMM1 + YMM0.
    ; Assuming YMM0 incorporates y3, y2, y1, y0, and YMM1 incorporates x3, x2, x1, x0, 
    ; YMM1 now incorporates as qwords y3+x3, y2+x2, y1+x1, y0+x0

    ; Use horizontal addition. I am not conscious of any model of vhaddpd for AVX-512.
    ; So I am restricted to working with YMM registers, however not ZMM registers. 
    vhaddpd ymm1, ymm1, ymm1
    ; After vhaddpd, YMM1 now incorporates as qwords y3+x3+y2+x2 (two instances), y1+x1+y0+x0 (two instances).

    vmovapd ymm2, ymm1 ; Copy qwords in YMM1 to YMM2

    ; Permute the bits in ymm2, primarily shifting the quadword at index 2 to index 0, 
    ; with all different indexes left unchanged. After permutation, 
    ; YMM2 will comprise y3+x3+y2+x2 at index 0, and YMM1 nonetheless has y1+x1+y0+x0 at index 0.
    ; Indexes 1, 2, and three are do not care.
    vpermpd ymm2, ymm2, 2 

    ; Add YMM1 and YMM2, storing end in YMM0. 
    vaddpd ymm0, ymm1, ymm2    

    ; At this level, YMM0's decrease half, i.e., XMM0, is able to be returned. Its decrease half, a double,
    ;    which is what the caller of this routine is anticipating, incorporates the sum of all components of the enter vector.
; Epilogue
    pop rsi
    add rsp, 20h      ; Regulate the stack again to unique state
    pop rdi           ; Restore RDI 
    ret 
    sumVector ENDP

How sumVector works

The sumVector routine is kind of a bit tougher to explain as a result of including the weather of a ZMM register will not be so simple as pairwise arithmetic on two such registers.

Prologue and Epilogue sections

As was the case with the beforehand described vecProduct routine, the Prologue and Epilogue sections are boilerplate. The descriptions are much like these for vecProduct.

Initialization part

The part with the label Initialization takes care of essential chores, resembling transferring values from the enter registers to different registers and initializing the AVX-512 registers.

PartialSumsLoop part

The code within the part marked PartialSumsLoop cycles via the vector being summed, working with 16 double values in every loop iteration, utilizing two ZMM registers. The code accumulates 16 double values into register ZMM0. As a substitute of working with solely 8 double values per loop iteration, I’m “unrolling” the loop barely within the curiosity of conserving the instruction pipeline as full as doable. Determine 3 under offers a visible clarification of how the 16 double values are added to the register that’s appearing because the accumulator. The loop is completed when the RCX register lastly reaches zero.

Fig. 3 Vector addition

After this system exhausts the entire components of the vector to be summed, a 512-bit ZMM register incorporates 8 partial sums. The following step is so as to add these 8 double values to get the grand complete. Oddly sufficient, this isn’t a simple matter, principally as a result of lack of an instruction that may add all 8 double values in a 512-bit ZMM register. Luckily, there’s an instruction that may work with the 4 double values in a 256-bit YMM register.

CombineSums part

There’s a number of motion occurring within the CombineSums part.

• Perform a partial addition by decreasing the variety of values to be added.
• Carry out horizontal addition to get a double that incorporates the sum of all components of an enter vector.
• Arrange the XMM0 register with the routine’s return worth, specifically the double end in its decrease 64 bits.

The next subsections describe every of those bigger steps.

Shrink the variety of values to be added

The part marked CombineSums carries out the steps essential so as to add the 8 partial sums talked about within the earlier part.
Step one is to separate up the 8 double values in ZMM0 between two YMM registers, with 4 double values within the YMM0 register, and 4 extra within the YMM1 register. A part of the work is already completed as a result of the decrease half (decrease 256 bits) of ZMM0 is aliased as YMM0. The code makes use of the vextractf64x4 instruction to extract the 4 double values from the higher half (index 1) of ZMM0, inserting them into YMM1. The f64x4 suffix implies that 4 floating-point 64-bit values (i.e., 4 double values) are to be extracted. The third operand right here, 1, signifies the higher half of the ZMM register. An operand of 0 would imply that we want to extract the values within the decrease half. Determine 4 under depicts the operation of this instruction, with YMM1 as proven, and YMM0 being the decrease half of ZMM0 (in orange).

Fig. 4 Extract particular elements of a vector

After the code splits up the 8 partial sums into 4 partial sums every within the YMM0 and YMM1 registers, the following step is so as to add these two registers collectively, utilizing the vaddpd instruction. Which means this system now wants so as to add solely 4 components of a YMM register, fairly than the eight components of a ZMM register to get the identical grand complete. The left portion of Determine 5 reveals the motion of vaddpd, with the 4 components every of registers YMM1 and YMM2 being added after which saved in YMM0. The permute portion of Determine 5 shall be mentioned later.

Fig. 5 Vector addition and vector permutation

Carry out horizontal addition

At this level, the 4 partial sums that make up the grand complete are in register YMM1. This system then makes use of the vhaddpd instruction to carry out a horizontal addition of YMM1 with itself, with the outcomes saved in YMM1. The center a part of the instruction identify, hadd, is brief for horizontal addition.

Determine 6 reveals this operation. From the higher left YMM occasion, this operation takes the 2 left-most double values, provides them, and shops this worth at index 3 within the backside YMM occasion. From the higher proper YMM occasion, the operation takes the 2 left-most values, provides them, after which shops precisely the identical worth at index 2 within the backside YMM occasion. The 2 lowest values within the higher left YMM occasion are added after which saved at index 1 within the decrease YMM occasion. The 2 lowest values within the higher proper YMM occasion are added after which saved at index 0 within the decrease YMM occasion. As a result of this system is performing horizontal addition on a register with itself, the 2 values at indexes 3 and a pair of on the backside are equivalent, as are these at indexes 1 and 0.

Fig. 6 Horizontal addition of a register with itself

Let’s take a second to consider what’s taking place right here. The aim is so as to add, say, 3 + 6 + 9 + 12, which sums to 30. After the horizontal addition, YMM1 incorporates 9, 9, 21, and 21, studying from left to proper. If we will by some means add the 9 and 21, we’ll have our appropriate sum, specifically 30.

Arrange the XMM0 return worth register

Should you’re nonetheless with me, congratulations! We’re nearly completed! At this level, YMM1 incorporates <9, 9, 21, 21>, studying left to proper. The one steps remaining are:

• Copy YMM1 to YMM2, utilizing the instruction vmovapd ymm2, ymm1. Subsequently, YMM1 incorporates <9, 9, 21, 21> as does YMM2.
• Permute the bits in YMM2 in order that the double at index 2 finally ends up at index 0. The code does this by utilizing the instruction vpermpd ymm2, ymm2, 2, which is depicted in the appropriate half of Determine 4. As I’ve used it in this system, this permutation actually is extra of a shift operation. The instruction merely shifts the 64 bits at index 2 of YMM2 to index 0, leaving the opposite three double values unchanged.
  YMM1 continues to be <9, 9, 21, 21> however YMM2 is now <9, 9, 21, 9>.
• Add YMM1 and YMM2, with the outcome being written to YMM0, utilizing vaddpd ymm0, ymm1, ymm2. The results of this motion is that YMM0 is <18, 18, 42, 30>. Be aware that the sum we needed is now at index 0 in YMM0.

Do not forget that I stated that including the parts of a vector was a sophisticated operation!

The decrease half of YMM0 is also called XMM0, which is <42, 30>, with 42 within the higher half of XMM0 and 30 within the decrease half of XMM0.

When sumVector completes execution, it returns the sum of all components of the enter vector in XMM0. Nevertheless, because the caller is anticipating solely the double that’s within the decrease half of XMM0, any worth that is perhaps current within the higher half of XMM0 is invisible to the caller.

Timing comparisons – CUDA vs. AVX-512

All vectors used within the two packages consisted of 262,144 components of sort double, with each packages compiled in Launch mode.

For timing, I used the high_resolution_clock::time_point class of the chrono namespace, which has a decision as small as 1 nanosecond. I don’t present the code I used, however in case you’re , observe the hyperlink on the finish of this text to my Github repository.

I ran each packages on a Dell Precision 7820 laptop, with an Intel Xeon(R) Silver 4144 CPU.  The nominal clock fee is 2.20 GHz. The CPU has 10 cores, however I didn’t try to reap the benefits of the a number of cores.

I ran the CUDA model on the NVidia Quadro P2000 GPU (Pascal structure). This GPU has 8 multiprocessors with a complete of 1024 cores and has 5120 MB of worldwide reminiscence. The cardboard’s GPU clock fee is 1481 MHz, and its reminiscence clock fee is 3504 MHz.

Calculate term-by-term vector merchandise

The next desk incorporates the mixed instances for the term-by-term vector merchandise ##< X_iY_i>## and ##< X_i^2>## for six runs of the CUDA code and the AVX-512 code. All instances are in microseconds.

These outcomes shocked me. I believed that the AVX-512 code could be enormously outperformed by the CUDA code, however in a lot of the runs, that was not the case. The possible motive for a way a lot time the CUDA code took is that a number of the general time is expended in copying vectors from host reminiscence to system reminiscence after which copying them again from system reminiscence to host reminiscence. The bottleneck is probably going these reminiscence transfers.

Calculating vector sums

The next desk incorporates the mixed instances for the 4 vector sums ##sum X_i##,  ##sum Y_i##, ##sum X_iY_i##, and ##sum X_i^2##, for six runs of the CUDA code and the AVX-512 code. All instances are in microseconds.

The CUDA program was at a drawback in calculating the vector sums as a result of it used the CPU to calculate the sums fairly than the GPU.

Full code

For the sake of brevity, I haven’t included the whole supply code right here on this article. In case your curiosity is piqued, and you’re a glutton for punishment, you’ll find the supply code for this text, Main2.cpp, and avxTest.asm right here: https://github.com/Mark44-AVX/CUDA-vs.-AVX-512.