fortran向量化的極限

Question

我已經寫了一個函數，返回一個向量A，等於稀疏矩陣Sparse與另一個向量F的乘積。矩陣的非零值在Sparse（nnz），rowind（nnz）和colind（nnz）每個都包含Sparse每個特定值的行和列。通過do kx下的兩行矢量化（現在註釋的）內部循環相對簡單....我看不出如何矢量化外部循環，因爲pos對於不同的kx具有不同的大小。

問題是：外部循環（do kx = 1，nxy）是否可以被矢量化，如果是，怎麼辦？

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%

弗拉基米爾F正確地推測說我來自Python /八度世界。由於我解決的PDE變得越來越大，我已經（移回）fortran以獲得更多的性能。在半小時前，矢量化意味着擺脫do循環，fortran看起來非常擅長：節省的時間用於替換「內部循環」（do ky = 1，size（pos）..）上面兩行很驚人。我在調用-fopt-info時查看gfortran（真正的gcc？）給出的信息，並且經常使用循環修改。我將立即閱讀關於SIMD和數組表示法。請，如果有關於這個話題的好消息，請讓我知道。

在回覆Holz時，存儲稀疏矩陣的方法很多，通常會導致操作符的等級降低1：我製作的示例涉及在某個字段中的每個位置評估的強制和解矢量，因此具有等級1.那麼與（S，如在A = S.F）相關的算子是二維BUT稀疏。它以這種方式存儲，只保留非零值。如果S中有nnz個非零值，那麼與S相當的稀疏Sp就是Sp（1：nnz）。如果pos表示某個數字Sp（pos）的該序列內的位置，則原始矩陣S中的列和行位置由colind（pos）和rowind（pos）給出。

在這樣的背景下，我可能會擴大問題：什麼是最好的（通過執行時間測量）可以完成乘法？

pure function SparseMul(Sparse,F) result(A) 
    implicit none 
    integer (kind=4),allocatable :: pos(:) 
    integer (kind=4) :: kx,ky ! gp counters 
    real (kind=8),intent(in) :: Sparse(:),F(:) 
    real (kind=8),allocatable :: A(:) 
    allocate(A(nxy)) 
    do kx=1,nxy     !for each row 
    pos=pack([(ky,ky=1,nnz)],rowind==kx) 
    A(kx)=sum(Sparse(pos)*F(colind(pos))) 
!!$  A(kx)=0 
!!$  do ky=1,size(pos) 
!!$   A(kx)=A(kx)+Sparse(pos(ky))*F(colind(pos(ky))) 
!!$  end do 
    end do 
end function SparseMul 

Answer 1

我假設問題「原樣」，即：

• 我們不想改變矩陣存儲格式
• 我們不想使用外部庫來執行任務

否則，我覺得用一個外部庫應該 是解決問題的最佳方法，例如 https://software.intel.com/en-us/node/520797。

要預測「最佳」Fortran寫乘法的方法並不容易。它取決於幾個因素（編譯器，體系結構， 矩陣大小，...）。我認爲最好的策略是建議 一些（合理的）嘗試，並以現實的配置進行測試。使用包

 
do kx=1,nxy 
    pos=pack([(ky,ky=1,nnz)],rowind==kx) 
    A(kx)=0 
    do ky=1,size(pos) 
     A(kx)=A(kx)+Sparse(pos(ky))*F(colind(pos(ky))) 
    end do 
end do 

1. 保存非零的位置：

   如果我理解正確的矩陣存儲格式，我嘗試 - - 包括那些在問題報告中提供以下

2. 作爲前一個但使用Fortran陣列語法

   do kx=1,nxy 
       pos=pack([(ky,ky=1,nnz)],rowind==kx) 
       A(kx)=sum(Sparse(pos)*F(colind(pos))) 
   end do 
   

3. 使用條件來確定的組件中使用

   do kx=1,nxy 
       A(kx)=0 
       do ky=1,nnz 
        if(rowind(ky)==kx) A(kx)=A(kx)+Sparse(ky)*F(colind(ky)) 
       end do 
   end do 
   

4. 與前一個但交換回路

   A(:)=0 
   do ky=1,nnz 
       do kx=1,nxy 
       if(rowind(ky)==kx) A(kx)=A(kx)+Sparse(ky)*F(colind(ky)) 
       end do 
   end do 
   

5. 使用intrisic sum與掩碼參數

   do kx=1,nxy  
       A(kx)=sum(Sparse*F(colind), mask=(rowind==kx)) 
   enddo 
   

6. 作爲前一個，但使用隱含的do循環

   A =[(sum(Sparse*F(colind), mask=(rowind==kx)), kx=1,nxy)] 
   

這些是使用具有 33％的非零值的1000×1000矩陣的結果。該機器是英特爾至強 和我的測試使用英特爾v17和GNU 6.1編譯器 執行使用沒有優化，高優化但沒有矢量化 和高優化。

  V1 V2 V3 V4 V5 V6 
-O0 
ifort 4.28 4.26 0.97 0.91 1.33 2.70 
gfortran 2.10 2.10 1.10 1.05 0.30 0.61 

-O3 -no-vec 
ifort 0.94 0.91 0.23 0.22 0.23 0.52 
gfortran 1.73 1.80 0.16 0.15 0.16 0.32 

-O3 
ifort 0.59 0.56 0.23 0.23 0.30 0.60 
gfortran 1.52 1.50 0.16 0.15 0.16 0.32 

對結果短短評論：

• 版本3-4-5通常是最快的人
• 編譯器優化的作用是任何版本
• 矢量化關鍵似乎只對 非最佳版本發揮重要作用
• 版本4是最好的兩個編譯器
• gfortran V4是「最佳」版本
• 優雅並不意味着良好的性能（V6不是很好）
• 附加註釋可以做分析 編譯器優化的報告。

如果我們有多核機器，我們可以嘗試使用 所有核心。這意味着處理代碼並行化，其中 是一個廣泛的問題，但只是給一些提示讓我們測試 兩種可能的OpenMP並行化。我們在 系列最快版本上工作（儘管不能保證它 也是最好的並行版本）。

的OpenMP 1.

!$omp parallel 
!$omp workshare 
    A(:)=0 
!$omp end workshare 
!$omp do 
    do ky=1,nnz 
     do kx=1,nxy     !for each row 
     if(rowind(ky)==kx) A(kx)=A(kx)+Sparse(ky)*F(colind(ky)) 
     end do 
    end do 
!$omp end do 
!$omp end parallel 
</pre> 

的OpenMP 2.附加FIRSTPRIVATE只讀矢量，以改善存儲器訪問

!$omp parallel firstprivate(Sparse, colind, rowind) 
    ... 
!$omp end parallel 

這些是16個核的結果爲多達16個線程：

#threads 1  2  4  8 16 
OpenMP v1 
ifort 0.22 0.14 0.088 0.050 0.027 
gfortran 0.155 0.11 0.064 0.035 0.020 
OpenMP v2 
ifort 0.24 0.12 0.065 0.042 0.029 
gfortran 0.157 0.11 0.052 0.036 0.029 

考慮到的可擴展性（在16個線程處大約8個）是合理的它是一個內存限制計算。第一個私有優化 僅對少量線程有優勢。 gfortran使用 16個線程是「最佳」 OpenMP解決方案。

Answer 2

我很難看到COLIND在哪裏以及它在做什麼......還有KX和KY。對於你想要向量化的內部循環，對於我使用OpenMP SIMD REDUCTION來說這似乎是最簡單的。我特意看這裏：

!!$  A(kx)=0 
!!$  do ky=1,size(pos) 
!!$   A(kx)=A(kx)+Sparse(pos(ky))*F(colind(pos(ky))) 
!!$  end do 

如果你必須收集（PACK），那麼它可能沒有多大幫助。如果F中有超過7/8的零，那麼F可能比PACK更好。否則，最好將矢量乘以一切（包括零和）。

主要規則是數據需要是連續的，所以你不能在第二維矢量化......如果感覺像Sparse和F是rank = 2，但它們顯示爲RANK = 1。即使它們確實是一個rank = 2的數組，也可以作爲向量來使用。 UNION/MAP也可用於實現二維陣列，也是一維矢量。

稀疏和F真的等級= 1？ nax，nay，nxy和colind用於什麼？並且其中許多沒有定義（例如，nn，nnz和colind）