What really is vectorization and how to implement it?

Uncategorized
1

Sorry for this perhaps stupid question.
I always feel not exactly clear what is vectorization, never really and fully understand it. I know there are many materials I can find on internet, google. But still not very clear what is vectorization.

Could anyone give some small examples? Thanks in advance!

Like, all the below code are using -O3 -xHost (or -Ofast -march=native in gfortran) optimization, and we are using a AVX2 CPU (256 bit register or something whatever). Actually all the below examples are using Intel Fortran, I am not sure if gfortran behave simialrly but I guess so.

Anyway,
I know in the following code, sorry for the sloppy code, see example1

example1

integer, parameter :: r8=selected_real_kind(15,9)
subroutine f(a,b,c)
real(kind=r8) :: a(1000), b(1000), c(1000)
c(:)=a(:)+b(:)
return
end

The above code should be vectorized,right? The real is real 8 type, it should took 64 bit in the ram. For AVX2 it is with 256 bit register or something, The above code is like doing c(1:4) = a(1:4) + b(1:4) in one clock cycle or something, and 4 is because 256/64=4.

But why the hell many times I found that actually doing the loop is faster? see example 2,

example2

integer, parameter :: r8=selected_real_kind(15,9)
subroutine f(a,b,c)
real(kind=r8) :: a(1000), b(1000), c(1000)
integer :: i
do i=1,1000
    c(i)=a(i)+b(i)
endddo
return
end

I notice that the compiler in example1 may complain that the array are not aligned or something. Therefore actually example2 is faster than example1.
But I mean in principle both example1 and 2 should be the same speed right? In example2, the compiler should be smart enough to vectorize the code correspondingly.

Finally, in the below code, example3, is there any vectorization?

example3

integer, parameter :: r8=selected_real_kind(15,9)
subroutine f(a,b,c)
real(kind=r8) :: a(1000), b(1000), c(1000)
c(:) = exp(a(:)) + log(b(:))
return
end

or slightly complicated things like, example4
example4

integer, parameter :: r8=selected_real_kind(15,9)
subroutine f(a,b,c)
real(kind=r8) :: a(1000), b(1000), c(1000)
c(:) = exp(a(:)) + log(b(:)) + exp(a(:)) * log(b(:))
return
end

I mean is vectorization is only for those very basic operatons like, op = + - * /

a(:) = b(:) op c(:)

However, for something more complicated, like

 a(:) = log(b(:) + c(:)*exp(a(:))*ln(b(:)))

Then vectorization may not work?

It seems in many case using a do loop is faster than writing things like a(:)=b(:)+c(:). It seem the compiler are doing very good job or highly optimized in just doing the do loops.

I also noticed an author wrote a guide for optimization Fortran code, is what he said in the manual still true?

2 Likes
2

It is possible to write vectorized versions of elementary functions like exp and log that are faster than the scaler versions. That said calling them from a language without metaprogramming will be ugly and annoying.

3 
subroutine f(a,b,c)
    real(kind=r8),intent(in) :: a(:), b(:)
    real(kind=r8),intent(in) :: c(:)
    c = a + b
end

Assuming you have explicit interfaces (e.g. subroutine is in a module) I would have made the example as above. I have seen examples of compilers getting confused on optimisation so make things as clear as possible. If it is a whole array operation do not specify a slice for example and add intent. That said you simple examples should vectorise. The do loop example is simple if there were other operations in the loop the compiler might have to work hard or be quite clever to work out if it can be vectorised. Maybe some of the values change in more than one place in the loop etc. Give the compiler the best chance possible to optimise is my rule. Also look at the optimisation reports (if any) to see what was actually done…

3 Likes
4

You can experiment with the code compilers produce, here: Compiler Explorer

3 Likes
5

On first reading, there are some good points and some misunderstandings in there. Perhaps I will have time to comment further.

1 Like
6

That subroutine has a bug: c must be declared with intent(out) not in.

1 Like
7

@CRquantum One other thing that is probably contributing to the confusion here, is that there are two completely different things that are sometimes referred to as “vectorization”.

  1. Vectorization can mean the use of SIMD instructions like AVX, which allow the CPU to perform multiple operations in the time it would otherwise take to perform one, as you know.

  2. Vectorization can mean using array statements instead of loops. For example, some people refer to c(:) = a(:) + b(:) as “vectorized”, where an equivalent do loop (like in your example) would be “non-vectorized”. This difference doesn’t matter that much in Fortran, where do loops are still fast, but it matters a lot in interpreted languages like matlab and python, where a for loop might be tens or hundreds of times slower than using a numpy array statement.

But the key thing that I think is causing some of your confusion is that those two different uses of the word “vectorization” are completely different. Using an array statement doesn’t necessarily mean that the compiler is more likely to emit SIMD instructions than if you use a do loop.

As for why you sometimes find that do loops are faster than array statements, I don’t know for sure but I can speculate. At the end of the day, array statements are just syntactic shorthand for loops - in the simple example you posted, I think a compiler could very well emit the exact same assembly code for the two versions. But in a more complicated example, there could be differences. Suppose, for example, you have a series of array statements, and you’re comparing that with a do-loop that does a series of operations in the body of the loop. If the compiler emits a separate loop for each array statement, then the array statement version might well be slower than that do-loop version, because it has to walk over the data multiple times. See for example section 2.3.8 in the optimization guide that you linked (I haven’t read the whole guide yet so I won’t comment on its overall accuracy).

7 Likes
8

When you’re talking “vectorization” in terms of the compiler emitting vector SIMD instructions (e.g. AVX2), AFAIK it is only turned on at higher levels of compiler optimization (-O2 or above for nvfortran; -O3 or above for ifort) or has to be enabled manually (as in specifying -march and -mtune for gfortran).

1 Like
9

These are my comments on the PDF document

Page 7, last para:
In fortran, data in arrays are stored in such an order that the most left-hand side digit is growing first.

Not specified in modern Fortran but the practice is universal among compilers. “Array element order” is what is specified in the Standard and it makes sense to make array layout in linearly addressed memory as close to array element order as possible.

Page 8, penultimate para:
Unfortunately that implicit shaped declarations also instruct the compiler that these dummy arrays are not necessarily contiguous in memory.

Correct because the subroutine can now be called with actual arguments being sections of arrays (which may include vector subscripts).

Page 8
Implicit shaped declarations should not be used

Disagree (and the term is assumed-shape)! Instead you can use intrinsic IS_CONTIGUOUS to determine at runtime if the actual argument is contiguous, and dispatch to more specific code:

subroutine add(a,b,c)
  Real, Intent(Out), Target:: a(:)
  Real, Intent(In), Target:: b(:),c(:)
  If (IS_CONTIGUOUS(a).AND.IS_CONTIGUOUS(b).AND. IS_CONTIGUOUS(c)) Then
    Block 
      Real, Pointer, Contiguous :: p(:),q(:),r(:)
      p=>a;q=>b;r=>c
      p(:)=q(:)+r(:)
    End Block
  Else
    a(:)=b(:)+c(:)
  End If
End Subroutine

Advanced compilers might insert such code speculatively. You can easily check the generated code from some compilers in Compiler Explorer website. Intel compiler produced AVX2 SIMD code for the above loop.

Page 9
Do not use array syntax in computational loops

Disagree! Advanced compilers could perform loop fusion. Do not underestimate their power!

4 Likes
10

GCC 12 Enables Auto-Vectorization For -O2 Optimization Level

https://www.phoronix.com/scan.php?page=news_item&px=GCC-12-Auto-Vec-O2

5 Likes
11

It is -O2 for ifort.

12

There can be several reasons for this, and some have already been mentioned in this thread. Another possibility is that the compiler might not realize that a given expression can be evaluated without intermediate work arrays, so it allocates temporary memory, evaluates the expression using those temporary arrays, and then deallocates the temporary arrays after the expression has been evaluated. That allocation/deallocation process is relatively expensive, thus resulting in nonoptimal array syntax code. In contrast, when the the compiler sees the loop code, the programmer has already determined that no array workspace is required and has arranged the expression in such a way that the compiler can see the correct order of operations; then given that, the compiler can evaluate using scalar registers only, or it can compile to pipelined instructions, or it can optimize to avx vector instructions, or maybe even optimize to gpu code, depending on what is the underlying hardware.

This situation is a little ironic, in a sadly disappointing way. One of the reasons that fortran 88, fortran 8x, and eventually f90 took so long to approve was that many vendors had invested much effort in compiler optimization technology to compile f77 do-loop codes into efficient vector instructions (here, “vector instructions” includes actual vector instructions like on Cray, Cyber-205, Convex, or ETA-10 hardware and also pipelined instructions on VLIW and RISC hardware of that era). To give one type of example, many f77 compilers would recognize three nested do-loops as a matrix-matrix product, rearrange the loops automatically into the best order for the underlying hardware, and then compile to the optimal instruction sequence. Some vendors thought that adding array syntax to the language would obviate that compiler optimization effort, so they slow-walked their approval of the language, hoping to exploit their optimizing compiler features for the maximum amount of time. The irony is that array syntax was not the big benefit that they had feared, and now, some 40 years later, we are still talking about do-loop based code running faster than array syntax code and trying to understand what happened.

2 Likes
13

If one looks at the compilers for the earliest vector and SIMD machines, e.g., the Burroughs BSP and Illiac-IV, the Texas Instruments ASC, the CDC STAR-100, they all had various forms of array syntax. A colleague of mine who worked at CDC on the 7600 and CYBER-205 projects before moving to Cray told me of programs being ‘converted’ that had little in common with the original version “except the PROGRAM statement”… In the 1980s, folks like Thinking Machines also used array syntax.

Meanwhile over at Cray Research, for the CRAY-1 (ca 1976) the compiler group instead decided to focus on looking at normal DO loops and attempting to generate vector instructions for them. For a long time, CFT was limited to vectorizing loops that fit in a single basic block - e.g., no IF/GOTO statements, no non-intrinsic procedure calls, no I/O, etc. This worked for the CRAY-1 because it allowed users to dump old code into the compiler and get speedups compared to previous machines like the CDC 7600 - both due to the faster scalar processing, and also any additional speedups provided by whatever vector instructions the compiler could muster.

Scalar code on the CRAY-1 was ~2x faster than the 7600. If a loop could be vectorized, that portion of the code was 10x faster. Amdahl’s Law comes into play here. So there was then a game to profile codes to find the ‘hot’ DO loops and make sure they were coded in a way that the compiler could vectorize them, arrays were accessed so that no bank conflicts occurred, and so on.

CFT and its successor CFT77 stayed close to the Fortran 77 Standard and never did implement array syntax. They both did get very good at deconstructing DO loops (including those with IFs and GOTOs). It wasn’t until cft90 that array syntax was implemented - because F90 required it.

It is interesting to look at in retrospect. I was one who was really looking forward to F90 array syntax. However it turned out to have fairly minor impact on the language compared to modules, derived types, explicit interfaces, dynamic memory allocation (allocatables and pointers), free-form, long names, etc.

3 Likes