MIC アプリケーションの SIMD ベクトル化ループでベクトル強度 0.0 がなる問題

この記事は、インテル® デベロッパー・ゾーンに掲載されている「vectorization intensity 0.0 in SIMD vectorized loop for MIC application」の日本語参考訳です。

皆さんからアドバイスをいただきたくて投稿しました。

私は、インテル® Xeon Phi™ コプロセッサー向けの N 体シミュレーションを作成しています。現在作業中のサブルーチンは、ペアを合計して N 体系のエネルギーを計算します。以下の Fortran コードは O(n2) アルゴリズムの実装で、単純にすべての粒子のペアと np (粒子数) の二重和を計算しています。2 つの粒子が距離 rcut2 よりも近い場合はそのエネルギーを加算します。

subroutine nearest_int
    implicit none
    double precision :: dx,dy,dz
    double precision :: x1,y1,z1
    double precision :: x2,y2,z2
    double precision :: dr2,dr2i,dr6i,dr12i
    integer  :: i,j
    integer :: T1,T2,clock_rate,clock_max
    
    potential =  0.0d0


    call system_clock(T1,clock_rate,clock_max)

    !dir$ offload begin target(mic:0) in(position)

    !$omp parallel do schedule(dynamic) reduction(+:potential) default(private),&
    !$omp& shared(position,rcut2,np)
    do i = 1,np
       
       x1 = position(i)%x; y1 = position(i)%y; z1 = position(i)%z

       !dir$ simd reduction(+:potential)  
       do j = i+1,np
                    
          x2 = position(j)%x; y2 = position(j)%y; z2 = position(j)%z
          dx = x2-x1
          dy = y2-y1
          dz = z2-z1

          dr2 = dx*dx + dy*dy + dz*dz
      
          if(dr2.lt.rcut2)then
             
          
             dr2i = 1.0d0/dr2
             dr6i = dr2i*dr2i*dr2i
             dr12i = dr6i*dr6i
             
        
             potential = potential + 4.0d0*(dr12i-dr6i)
          endif
       enddo
    enddo
    !$omp end parallel do
   !dir$ end offload


    call system_clock(T2,clock_rate,clock_max)
    print*,'elapsed time nint:',real(T2-T1)/real(clock_rate),potential
  end subroutine nearest_int

ここで、position は次の構造体配列です。

type atom
    double precision :: x,y,z
end type atom

type(atom) :: position

この O(n2) アルゴリズムのコードを実行すると、ベクトル強度は 6.51 になります。ギャザー/スキャッター操作を行っていることを考慮すれば、悪くないと言えるでしょう。以下は、インテル® VTune™ Amplifier の [Summary] 画面のスクリーンショットです。

O(n2) は大規模なシステムではスケーリングが制限されるため、O(N) アルゴリズムに変更します。近い粒子 (1.2*rcut) とそれに隣接する粒子数を配列に格納し、ポインターを使って position 配列にアクセスするように、O(n2) アルゴリズムを変更します。

subroutine nearest_int
    implicit none
    double precision :: dx,dy,dz
    double precision :: x1,y1,z1
    double precision :: x2,y2,z2
    double precision :: dr2,dr2i,dr6i,dr12i
    integer  :: i,j
    integer :: T1,T2,clock_rate,clock_max
    integer :: neigh
    
    potential =  0.0d0


    call system_clock(T1,clock_rate,clock_max)

    !dir$ offload begin target(mic:0) in(position,vlistl,numneigh)

    !$omp parallel do schedule(dynamic) reduction(+:potential) default(private),&
    !$omp& shared(position,neigh_alloc,vlistl,numneigh,rcut2,np)
    do i = 1,np
       
       x1 = position(i)%x; y1 = position(i)%y; z1 = position(i)%z

       !dir$ simd reduction(+:potential)  
       do j = 1,numneigh(i)
          
         neigh = vlistl(j + neigh_alloc*(i-1))
         x2 = position(neigh)%x; y2 = position(neigh)%y; z2 = position(neigh)%z
          
          dx = x2-x1
          dy = y2-y1
          dz = z2-z1

          dr2 = dx*dx + dy*dy + dz*dz
      
          if(dr2.lt.rcut2)then
             
          
             dr2i = 1.0d0/dr2
             dr6i = dr2i*dr2i*dr2i
             dr12i = dr6i*dr6i
             
        
             potential = potential + 4.0d0*(dr12i-dr6i)
          endif
       enddo
    enddo
    !$omp end parallel do
   !dir$ end offload


    call system_clock(T2,clock_rate,clock_max)
    print*,'elapsed time nint:',real(T2-T1)/real(clock_rate),potential
  end subroutine nearest_int

ここで、vlistl は次のように割り当てています。

neigh_alloc = 500
allocate(vlistl(500*np))

-vec-report6 オプションを指定してコンパイルすると、内側のループがベクトル化されたことを示すメッセージが表示されますが、ベクトル強度は 0.0 で、パフォーマンスが大幅に低下します (シリアルよりもわずかに良いだけで、ベクトル強度が 0.0 であれば当然と言えるでしょう)。以下は、インテル® VTune™ Amplifier の [Summary] 画面のスクリーンショットです。

以下が、私が直面している問題です。

1. ベクトル化されたループの内側のベクトル強度が 0.0 になったのはなぜでしょうか? パフォーマンスを向上するにはどうしたら良いのでしょうか?

2. この O(n) アルゴリズムのコードを MIC で実行する場合、行 27 (neigh = vlistl(j + neigh_alloc*(i-1)) でレイテンシーの問題が発生することが予想されます。ここでプリフェッチを利用する場合、ギャザー操作は 16 キャッシュラインで 8 つの倍精度データを処理できますが、レイテンシーを隠蔽するにはどのようにプリフェッチを利用すべきでしょうか? 行 27 の直後に次のループを配置してみましたが、パフォーマンスは変わりませんでした。

do k = 0, 15
     call mm_prefecth(position(vlistl(j + neigh_alloc*(i-1) + k+8)%x,1)
enddo

3. すべての配列を 64 バイト境界でアライメントするため、-align array64byte オプションを指定してコンパイルしましたが、これによりキャッシュラインのサイズの倍数で配列がパディングされるのでしょうか?もしそうでなければ、どうしたらそうなりますか?

simd プラグマによりループはベクトル化されるはずです。粒子はソートして、空間で近くにあるものは、メモリーでも近くに配置しています。ベクトル化には配列構造体のほうが構造体配列よりも適していることは分かっていますが、構造体配列を使用した O(n2) アルゴリズムでは良いパフォーマンスが得られました。また、インテルでは、このアルゴリズムを利用するさまざまなソフトウェアでこのデータ構造を実装しています (例えば、lammps など。ただし、Fortran と C++ という言語の違いがあります)。

添付のコードを利用して完全なコードをコンパイルできます。問題のサブルーチンは、mod_force.f90 にあります。numneigh 配列と vlistl 配列は、mod_neighbor.f90 にある build_neighbor_n2 サブルーチンで作成し、mod_neighbor.f90 にある init_list サブルーチンで割り当てています。すべての配列は、module global.f90 でグローバルとして定義しています。また、次のコマンドを実行してコンパイルしました。

ifort -align array64byte -openmp global.f90 get_started.f90 mod_init_posit.f90 mod_neighbor.f90 mod_force.f90 MD.f90 -O3 -o new.out

長くなりましたが、これまで私が試したことがお分かりいただけたかと思います。皆さんからのアドバイスをお待ちしています。