Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.

Vectorization and OpenMP threadprivate

AONym
New Contributor II
908 Views

I'm trying to vectorize this loop in Fortran XE2019 u5:

DO iRadius = 1, nRadii
    Grev(iRadius, 1)=Dradius(iRadius-1)*dt*A1rev(iRadius, 1)+Brev(iRadius, 1)-sRadius(iRadius-1)*omSq(iRun)*dt*A2rev(iRadius, 1)
    Grev(iRadius, 2)=Dradius(iRadius)*dt*A1rev(iRadius, 2)+Brev(iRadius, 2)-sRadius(iRadius)*omSq(iRun)*dt*A2rev(iRadius, 2)
    Grev(iRadius, 3)=Dradius(iRadius+1)*dt*A1rev(iRadius, 3)+Brev(iRadius, 3)-sRadius(iRadius+1)*omSq(iRun)*dt*A2rev(iRadius, 3)
END DO L_R2_Serial

All the variables except iRadius are REAL(8). Grev is declared in this subroutine as DIMENSION(10000, 3), INTENT(OUT). A1rev, A2rev and Brev are declared in a separate module as DIMENSION(10000, 3). There is also the OpenMP line

!$OMP THREADPRIVATE(Brev, A1rev, A2rev) 

Compiling with /O3 in VisualStudio 2019 gives, in the optimization report

LOOP BEGIN at ...MyProgram.f90(414,1)
...MyProgram.f90(488,9):remark #25084: Preprocess Loopnests: Moving Out Store 
...MyProgram.f90(414,1):remark #17104: loop was not parallelized: existence of parallel dependence
...MyProgram.f90(414,1):remark #17106: parallel dependence: assumed FLOW dependence between GREV(IRADIUS,1) (420:13) and A1REV (427:13)
...MyProgram.f90(414,1):remark #17106: parallel dependence: assumed ANTI dependence between A1REV (427:13) and GREV(IRADIUS,1) (420:13)
...MyProgram.f90(414,1):remark #15344: loop was not vectorized: vector dependence prevents vectorization
...MyProgram.f90(414,1):remark #15346: vector dependence: assumed FLOW dependence between GREV(IRADIUS,1) (420:13) and A1REV (427:13)
...MyProgram.f90(414,1):remark #15346: vector dependence: assumed ANTI dependence between A1REV (427:13) and GREV(IRADIUS,1) (420:13)
...MyProgram.f90(414,1):remark #25015: Estimate of max trip count of loop=10002
LOOP END

However, if I remove the OMP threadprivate line, I get

LOOP BEGIN at ...MyProgram.f90(414,1)
...MyProgram.f90(488,9):remark #25084: Preprocess Loopnests: Moving Out Store 
...MyProgram.f90(420,13):remark #15388: vectorization support: reference GREV(IRADIUS,1) has aligned access
...MyProgram.f90(420,13):remark #15388: vectorization support: reference DRADIUS(IRADIUS-1) has aligned access
...MyProgram.f90(420,48):remark #15388: vectorization support: reference A1REV(IRADIUS,1) has aligned access
...MyProgram.f90(420,51):remark #15388: vectorization support: reference BREV(IRADIUS,1) has aligned access
...MyProgram.f90(420,69):remark #15388: vectorization support: reference SRADIUS(IRADIUS-1) has aligned access
...MyProgram.f90(420,116):remark #15388: vectorization support: reference A2REV(IRADIUS,1) has aligned access
...MyProgram.f90(423,13):remark #15388: vectorization support: reference GREV(IRADIUS,2) has aligned access
...MyProgram.f90(423,13):remark #15388: vectorization support: reference DRADIUS(IRADIUS) has aligned access
...MyProgram.f90(423,46):remark #15388: vectorization support: reference A1REV(IRADIUS,2) has aligned access
...MyProgram.f90(423,49):remark #15388: vectorization support: reference BREV(IRADIUS,2) has aligned access
...MyProgram.f90(423,67):remark #15388: vectorization support: reference SRADIUS(IRADIUS) has aligned access
...MyProgram.f90(423,112):remark #15388: vectorization support: reference A2REV(IRADIUS,2) has aligned access
...MyProgram.f90(427,13):remark #15388: vectorization support: reference GREV(IRADIUS,3) has aligned access
...MyProgram.f90(427,13):remark #15388: vectorization support: reference DRADIUS(IRADIUS+1) has aligned access
...MyProgram.f90(427,48):remark #15388: vectorization support: reference A1REV(IRADIUS,3) has aligned access
...MyProgram.f90(427,51):remark #15388: vectorization support: reference BREV(IRADIUS,3) has aligned access
...MyProgram.f90(427,69):remark #15388: vectorization support: reference SRADIUS(IRADIUS+1) has aligned access
...MyProgram.f90(427,116):remark #15388: vectorization support: reference A2REV(IRADIUS,3) has aligned access
...MyProgram.f90(414,1):remark #15305: vectorization support: vector length 4
...MyProgram.f90(414,1):remark #15309: vectorization support: normalized vectorization overhead 0.400
...MyProgram.f90(414,1):remark #15300: LOOP WAS VECTORIZED
...MyProgram.f90(414,1):remark #15448: unmasked aligned unit stride loads: 15 
...MyProgram.f90(414,1):remark #15449: unmasked aligned unit stride stores: 3 
...MyProgram.f90(414,1):remark #15475: --- begin vector cost summary ---
...MyProgram.f90(414,1):remark #15476: scalar cost: 65 
...MyProgram.f90(414,1):remark #15477: vector cost: 11.250 
...MyProgram.f90(414,1):remark #15478: estimated potential speedup: 5.770 
...MyProgram.f90(414,1):remark #15488: --- end vector cost summary ---
...MyProgram.f90(414,1):remark #25015: Estimate of max trip count of loop=2500
LOOP END

which is what I'm trying to achieve. Unfortunately, due to parallelization in the caller of this subroutine, A1rev, A2rev and Brev need to be different for each thread.

I don't understand why the optimizer is seeing this flow and anti flow dependence. In fact, A1rev etc are not written in this subroutine, only in the caller.

Can someone suggest what is necessary to get the vectorization I am looking for?

0 Kudos
9 Replies
jimdempseyatthecove
Honored Contributor III
908 Views

Not sure why the difference. As a work around, the procedure could be passed in Brev, A1rev, A2rev as dummy arguments (tpBrev, tpA1rev, tpA2rev), then those used in the loop.

Note, it is difficult to assess OpenMP issues without seeing the full context of the !$OMP directives. In this case "Preprocess Loopnests: Moving Out Store" is suspicious of a nested loop. Can you rewrite your code sketching anything related to how the parallel region, inclusive of nesting, is organized (you can elide non-pertinent code using ...).

Jim Dempsey

0 Kudos
AONym
New Contributor II
908 Views
@jimdempseyatthecove: Here's an outline showing the OMP structure:
Subroutine GOF:

!$OMP PARALLEL DEFAULT(SHARED)
!$OMP DO PRIVATE(iRun, nRadii, line,  ...
DO iRun=1, nRuns
...
   CALL CS(..., iRun)
...
END DO ! iRun
!$OMP END DO
!$OMP BARRIER
!$OMP END PARALLEL
-----------------------------------------------
Subroutine CS:

REAL(8), DIMENSION(10000, 3) :: Grev
nThreads= GC(... , Grev)
-----------------------------------------------
Function GC:

USE CI
USE RI
REAL(8), DIMENSION(10000), INTENT(OUT) :: Grev
REAL(8), DIMENSION(0:10001) :: Dradius, sRadius
...
L_R2_Serial: &
    DO iRadius = 1, nRadii
        Grev(iRadius, 1)=Dradius(iRadius-1)*dt*A1rev(iRadius, 1)+Brev(iRadius, 1)-sRadius(iRadius-1)*omSq(iRun)*dt*A2rev(iRadius, 1)
        Grev(iRadius, 2)=Dradius(iRadius)*dt*A1rev(iRadius, 2)+Brev(iRadius, 2)-sRadius(iRadius)*omSq(iRun)*dt*A2rev(iRadius, 2)
        Grev(iRadius, 3)=Dradius(iRadius+1)*dt*A1rev(iRadius, 3)+Brev(iRadius, 3)-sRadius(iRadius+1)*omSq(iRun)*dt*A2rev(iRadius, 3)
    END DO L_R2_Serial
-----------------------------------------------
Module CI:

!DIR$ IF DEFINED(ALIGNED)
!DIR$ attributes align : 64 :: Brev, A1rev, A2rev
!DIR$ ENDIF
REAL(8), DIMENSION(10000, 3 ) :: Brev, A1rev, A2rev
!$OMP THREADPRIVATE(Brev, A1rev, A2rev)   ! **** if this is removed, vectorization is successful
REAL(8) :: dt
-----------------------------------------------
Module RI:

REAL(8), DIMENSION(32) :: omSq


 

0 Kudos
jimdempseyatthecove
Honored Contributor III
908 Views

Shape of Grev on line 23 does not match shape used on line 28.

As sketched, it is unknown as to if Grev on line 16 is a local array to subroutine CS or a dummy argument of CS or a global (module) array.

If Grev on line 16 references the same memory for all threads, then the DO loop on line 27 is not correct.

Do CS or GC contain !$OMP directives, if so, include this with your sketch.

Jim Dempsey

0 Kudos
John_Campbell
New Contributor II
908 Views

Would the following change help ? It could have the benefit of processing memory sequentially

do k = 1,3
  DO iRadius = 1, nRadii
    Grev(iRadius, k)=Brev(iRadius, k) + dt*( Dradius(iRadius+k-2)*A1rev(iRadius, k) - sRadius(iRadius+k-2)*A2rev(iRadius, k)*omSq(iRun) )
  END DO
end do

 

0 Kudos
AONym
New Contributor II
908 Views

@jimdempseyatthecove: I made a mistake in copying. The correct declaration of Grev in GC is

REAL(8), DIMENSION(10000, 3), INTENT(OUT) :: Grev

As far as !$OMP directives, CS does contain a few !$OMP CRITICAL / !$OMP END CRITICAL surrounding debugging outputs that should be contiguous for each thread, but nothing else. GC did contain a number or them, but I removed them all for this test.

@John Campbell: This is a good idea. The original code was DO iRadius / DO k, then it was changed to reverse the subscripts of Grev, then the k loop was unrolled, then the subscripts were reversed back to the original!

0 Kudos
AONym
New Contributor II
908 Views

In addition to implementing John Campbell's suggestion, I used jimdempseyatthecove's idea of changing A1rev, A2rev and Brev from threadprivate in module CI to local in CS, and passed as arguments to all subroutines that used them, including GC. Now the code is

DO iRadius = 1, nRadii
    DO j = 1, 3
        k = j - 2 ! -1, 0, +1
        Grev(iRadius, j)=Dradius(iRadius+k)*dt*A1rev(iRadius, j) + Brev(iRadius, j) - sRadius(iRadius+k)*omSq(iRun)*dt*A2rev(iRadius, j)
    END DO
END DO

and the optimization report shows vectorization as shown in the second report in my initial post.

0 Kudos
John_Campbell
New Contributor II
908 Views

Is there any difference if you have the alternative DO order:

DO j = 1, 3
   k = j - 2 ! -1, 0, +1
   DO iRadius = 1, nRadii

This would process sequentially and could be faster for AVX for large nRadii, although could already be changed by optimisation ?

0 Kudos
AONym
New Contributor II
908 Views

@John Campbell: The loop order you suggested (j, then iRadius) is slightly faster (9:06 vs 8:48 m:ss), but there is something strange in the optimization report.

For DO iRadius / DO j, I get

LOOP BEGIN at ...MyProgram.f90(416,1)
...MyProgram.f90(504,5):remark #25084: Preprocess Loopnests: Moving Out Store 
...MyProgram.f90(425,13):remark #15388: vectorization support: reference GREV(IRADIUS,J) has aligned access
...MyProgram.f90(425,13):remark #15388: vectorization support: reference DRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(425,52):remark #15388: vectorization support: reference A1REV(IRADIUS,J) has aligned access
...MyProgram.f90(425,72):remark #15388: vectorization support: reference BREV(IRADIUS,J) has aligned access
...MyProgram.f90(425,70):remark #15388: vectorization support: reference SRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(425,124):remark #15388: vectorization support: reference A2REV(IRADIUS,J) has aligned access
...MyProgram.f90(425,13):remark #15388: vectorization support: reference GREV(IRADIUS,J) has aligned access
...MyProgram.f90(425,13):remark #15388: vectorization support: reference DRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(425,52):remark #15388: vectorization support: reference A1REV(IRADIUS,J) has aligned access
...MyProgram.f90(425,72):remark #15388: vectorization support: reference BREV(IRADIUS,J) has aligned access
...MyProgram.f90(425,70):remark #15388: vectorization support: reference SRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(425,124):remark #15388: vectorization support: reference A2REV(IRADIUS,J) has aligned access
...MyProgram.f90(425,13):remark #15388: vectorization support: reference GREV(IRADIUS,J) has aligned access
...MyProgram.f90(425,13):remark #15388: vectorization support: reference DRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(425,52):remark #15388: vectorization support: reference A1REV(IRADIUS,J) has aligned access
...MyProgram.f90(425,72):remark #15388: vectorization support: reference BREV(IRADIUS,J) has aligned access
...MyProgram.f90(425,70):remark #15388: vectorization support: reference SRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(425,124):remark #15388: vectorization support: reference A2REV(IRADIUS,J) has aligned access
...MyProgram.f90(416,1):remark #15305: vectorization support: vector length 4
...MyProgram.f90(416,1):remark #15309: vectorization support: normalized vectorization overhead 0.400
...MyProgram.f90(416,1):remark #15300: LOOP WAS VECTORIZED
...MyProgram.f90(416,1):remark #15448: unmasked aligned unit stride loads: 15 
...MyProgram.f90(416,1):remark #15449: unmasked aligned unit stride stores: 3 
...MyProgram.f90(416,1):remark #15475: --- begin vector cost summary ---
...MyProgram.f90(416,1):remark #15476: scalar cost: 65 
...MyProgram.f90(416,1):remark #15477: vector cost: 11.250 
...MyProgram.f90(416,1):remark #15478: estimated potential speedup: 5.770 
...MyProgram.f90(416,1):remark #15488: --- end vector cost summary ---

But for DO j / DO iRadius, I get

LOOP BEGIN at ...MyProgram.f90(419,9)
...MyProgram.f90(421,9):remark #25084: Preprocess Loopnests: Moving Out Store 
...MyProgram.f90(419,9):remark #17108: loop was not parallelized: insufficient computational work
...MyProgram.f90(420,13):remark #15389: vectorization support: reference GREV(IRADIUS,J) has unaligned access
...MyProgram.f90(420,13):remark #15388: vectorization support: reference DRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(420,52):remark #15389: vectorization support: reference A1REV(IRADIUS,J) has unaligned access
...MyProgram.f90(420,72):remark #15389: vectorization support: reference BREV(IRADIUS,J) has unaligned access
...MyProgram.f90(420,70):remark #15388: vectorization support: reference SRADIUS(IRADIUS+J-2) has aligned access
...MyProgram.f90(420,124):remark #15389: vectorization support: reference A2REV(IRADIUS,J) has unaligned access
...MyProgram.f90(419,9):remark #15381: vectorization support: unaligned access used inside loop body
...MyProgram.f90(419,9):remark #15305: vectorization support: vector length 4
...MyProgram.f90(419,9):remark #15399: vectorization support: unroll factor set to 4
...MyProgram.f90(419,9):remark #15309: vectorization support: normalized vectorization overhead 0.171
...MyProgram.f90(419,9):remark #15300: LOOP WAS VECTORIZED
...MyProgram.f90(419,9):remark #15448: unmasked aligned unit stride loads: 2 
...MyProgram.f90(419,9):remark #15450: unmasked unaligned unit stride loads: 3 
...MyProgram.f90(419,9):remark #15451: unmasked unaligned unit stride stores: 1 
...MyProgram.f90(419,9):remark #15475: --- begin vector cost summary ---
...MyProgram.f90(419,9):remark #15476: scalar cost: 22 
...MyProgram.f90(419,9):remark #15477: vector cost: 10.250 
...MyProgram.f90(419,9):remark #15478: estimated potential speedup: 2.140 
...MyProgram.f90(419,9):remark #15488: --- end vector cost summary ---
...MyProg

One difference is the aligned vs unaligned arrays. However, the arrays are allocated aligned to 64 bytes and the first dimension is 10000, which is divisible by 8, so the columns are also 64-byte aligned. The only difference between the two builds is the reversal of the DO loops.

The other difference is the "unit stride loads"; not sure what this means.

0 Kudos
jimdempseyatthecove
Honored Contributor III
908 Views

Try this (untested code):

DO j = 1, 3
    k = j - 2 ! -1, 0, +1
    ! use ASSOCIATE to re-interpret 2D arrays as a 1D slice
    ASSOCIATE(qGrev=>Grev(:,j), qA1Rev=>A1rev(:,j), qBrev=>Brev(:,j), qA2Rev=>A2rev(:,j))
    ! inform the compiler of alignment (that you asserted was true)
    !DIR$ ASSUME_ALIGNED qGrev:64, qA1rev:64, qBrev:64, qA2rev:64
    DO iRadius = 1, nRadii
        qGrev(iRadius)=Dradius(iRadius+k)*dt*qA1rev(iRadius) + qBrev(iRadius) - sRadius(iRadius+k)*omSq(iRun)*dt*qA2rev(iRadius)
    END DO
    !DIR$ END ASSOCIATE
END DO

"unit stride" means an array slice has elements that are adjacent in memory A(1,J), A(2), A(3,J),... as opposed to A(I,1), A(I,2), A(I,3),...

Jim Dempsey

0 Kudos
Reply