Tuesday, June 14, 2022
HomePhysicsParallel Programming on a CPU with AVX-512

Parallel Programming on a CPU with AVX-512

This text is the second of a two-part sequence that presents two distinctly totally different approaches to parallel programming. Within the two articles, I take advantage of totally different approaches to resolve the identical drawback: 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 models (GPUs). In SIMT programming, a single instruction is executed concurrently on lots of of microprocessors on a graphics card.
  • Single-instruction a number of information (SIMD) as offered on x64 processors from Intel® and AMD® (this text). In SIMD programming, a single instruction operates on extensive registers that may include 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 applications.

Regression line calculations

Given a set of factors (Xi, Yi), the place i ranges from 1 to N, the slope m and y-intercept b of the regression line may 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 may be seen from these formulation, there are a lot of 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.

Much 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##.

Essential program duties

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

  • Declare exterior (i.e., assembly-coded) features and declare the arrays and variables that will probably 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 distinct 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.

Perform prototypes and variable declarations

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

The 4 arrays declared under characterize 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 substantiate that the calculated slope and y-intercept of the regression line are right. Within the code under that initializes the 2 coordinate vectors, I’ve “unrolled” the loop by calculating 4 units of factors per loop iteration, quite 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 include the term-by-term merchandise of the corresponding parts 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 in a similar way.

Slope and intercept calculations

In spite of everything 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

Quite than displaying 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 towards 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 notice on AVX-512 registers

An AVX-512 register comparable to ZMM0 is 64 bytes or 512 bits extensive and may include a vector of 8 double values, 16 float or int values, 32 brief values, or 64 char values. Its decrease half of ZMM0, which comprises 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 exhibits the relationships between ZMM, YMM, and XMM registers.

ZMM, YMM, and XMM registers

Fig. 1 ZMM, YMM, and XMM registers

Knowledge part

The information 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 kind double (64 bytes).


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 - deal with of output array, handed in RCX.
    ;     inArr1 - deal with of first enter array, handed in RDX.
    ;     inArr2 - deal with of second enter array, handed in R8.
    ;     arr_len - variety of parts of every array, handed in R9.
    ; Return worth: None.
    ; On return, RAX holds the deal with 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 reserve it.
    sub rsp, 20h 
    mov r10, rsp
    push rsi          ; Save RSI as nicely.

; Initialization 
    vzeroall          ; Zero out all ZMM registers.
    mov rax, rcx      ; Copy output array deal with 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.

    ; Iterate by means of 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         ; Alter the stack again to unique state.
    pop rdi              ; Restore RDI. 
vecProduct ENDP

How vecProduct works

Of the 2 meeting routines, vecProduct is quite a bit 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 only put aside some stack reminiscence after which later reclaim it.

Initialization part

The Initialization part zeroes all the 512-bit ZMM registers, and copies the output and enter vector addresses from the registers used to go them into totally different registers. It additionally units up the RCX register to the variety of iterations of a loop that will probably be used to course of all the values within the enter vectors. The loop will multiply 8 values of kind double from every of the 2 enter vectors in every iteration, and retailer the end result 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. (Observe that the principle 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 exhibits the outcomes of the vmulpd instruction, multiplying every corresponding double within the two supply registers, with the outcomes being saved in ZMM0.

Vector multiplication

Fig. 2 Vector multiplication

The ends in ZMM0 are then copied to the reminiscence for the output vector utilizing vmovapd once more, the place the output vector’s base vacation spot deal with 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 remains to be nonzero. When RCX lastly will get to zero, management drops by means of to the following line after the loop instruction.

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

sumVector routine

The sumVector routine sums all the parts 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 - deal with of an array of doubles, handed in RCX
    ;     array_len - variety of parts 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 reserve it.
    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 deal with to RAX.
    mov rcx, rdx ; Copy aspect rely to RCX.
    shr rcx, 4 ; RCX <- Variety of 2 X 8-double chunks.
    xor rsi, rsi ; Zero RSI. 
    ; In every iteration, accumulate 8 partial sums of doubles.
    ; When loop terminates, ZMM0 will maintain 8 doubles that, 
    ; when added, would be the total 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 

    ; Use vextractf64x4 to extract 4 doubles (256 bits) from the higher half of ZMM0 to YMM1.
    ; YMM0 already comprises the decrease 256 bits (4 doubles) of ZMM0.
    vextractf64x4 ymm1, zmm0, 1 
    vaddpd ymm1, ymm1, ymm0 ; YMM0 <-- YMM1 + YMM0.
    ; Assuming YMM0 comprises y3, y2, y1, y0, and YMM1 comprises x3, x2, x1, x0, 
    ; YMM1 now comprises 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 comprises as qwords y3+x3+y2+x2 (two occasions), y1+x1+y0+x0 (two occasions).

    vmovapd ymm2, ymm1 ; Copy qwords in YMM1 to YMM2

    ; Permute the bits in ymm2, basically shifting the quadword at index 2 to index 0, 
    ; with all different indexes left unchanged. After permutation, 
    ; YMM2 will include 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, comprises the sum of all parts of the enter vector.
; Epilogue
    pop rsi
    add rsp, 20h      ; Alter the stack again to unique state
    pop rdi           ; Restore RDI 
    sumVector ENDP

How sumVector works

The sumVector routine is kind of a bit harder to explain as a result of including the weather of a ZMM register isn’t 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 just like these for vecProduct.

Initialization part

The part with the label Initialization takes care of crucial chores, comparable to 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 by means of 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 an alternative 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 potential. Determine 3 under supplies a visible clarification of how the 16 double values are added to the register that’s performing because the accumulator. The loop is completed when the RCX register lastly reaches zero.

Vector addition

Fig. 3 Vector addition

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

CombineSums part

There’s loads of motion happening 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 comprises the sum of all parts 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 crucial 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, putting 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 components 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 that this system now wants so as to add solely 4 parts of a YMM register, quite than the eight parts of a ZMM register to get the identical grand complete. The left portion of Determine 5 exhibits the motion of vaddpd, with the 4 parts every of registers YMM1 and YMM2 being added after which saved in YMM0. The permute portion of Determine 5 will probably be mentioned later.

Vector addition, and vector permutation

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 exhibits 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 couple of on the backside are equivalent, as are these at indexes 1 and 0.

Horizontal addition

Fig. 6 Horizontal addition of a register with itself

Let’s take a second to consider what’s occurring right here. The purpose is so as to add, say, 3 + 6 + 9 + 12, which sums to 30. After the horizontal addition, YMM1 comprises 9, 9, 21, and 21, studying from left to proper. If we will in some way add the 9 and 21, we’ll have our right sum, specifically 30.

Arrange the XMM0 return worth register

Should you’re nonetheless with me, congratulations! We’re virtually completed! At this level, YMM1 comprises <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 comprises <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 through the use of the instruction vpermpd ymm2, ymm2, 2, which is depicted in the fitting 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 remains to be <9, 9, 21, 21> however YMM2 is now <9, 9, 21, 9>.
  • Add YMM1 and YMM2, with the end result being written to YMM0, utilizing vaddpd ymm0, ymm1, ymm2. The results of this motion is that YMM0 is <18, 18, 42, 30>. Observe that the sum we wished is now at index 0 in YMM0.

Keep in mind that I stated that including the elements of a vector was an advanced 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 parts of the enter vector in XMM0. Nevertheless, for the reason that caller is anticipating solely the double that’s within the decrease half of XMM0, any worth that could be current within the higher half of XMM0 is invisible to the caller.

Timing comparisons – CUDA vs. AVX-512

All vectors used within the two applications consisted of 262,144 parts of kind double, with each applications 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 for those who’re , comply with the hyperlink on the finish of this text to my Github repository.

I ran each applications on a Dell Precision 7820 laptop, with an Intel Xeon(R) Silver 4144 CPU.  The nominal clock charge is 2.20 GHz. The CPU has 10 cores, however I didn’t try and 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 world reminiscence. The cardboard’s GPU clock charge is 1481 MHz, and its reminiscence clock charge is 3504 MHz.

Calculate term-by-term vector merchandise

The next desk comprises the mixed occasions 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 occasions are in microseconds.

CUDA (usec.) AVX-512 (usec)
5688 4264
5468 6873
2804 3959
5586 3597
2745 3990
5632 4194
Common: 4653.8 microseconds Common: 4479.5 microseconds

These outcomes stunned me. I assumed that the AVX-512 code can be drastically outperformed by the CUDA code, however in many of the runs, that was not the case. The seemingly cause for the way a lot time the CUDA code took is that loads 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 comprises the mixed occasions 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 occasions are in microseconds.

CUDA (time in usec.) AVX512 (time in usec)
1463 654
1624 852
698 640
1519 616
724 640
1615 665
Common: 1273.8 microseconds Common 677.8 microseconds

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

Full code

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



Please enter your comment!
Please enter your name here

Most Popular

Recent Comments