在fortran中使用openmp并行创建稀疏矩阵

Posted 2023-02-19

【中文标题】在fortran中使用openmp并行创建稀疏矩阵

【英文标题】：Sparse matrix parallel creation with openmp in fortran

【发布时间】：2022-01-15 03:51:39

【问题描述】：

我对fortran比较陌生，对openmp完全陌生，我有以下问题：

我想并行构建一个（大：~1% 非零元素，总共约 100 万到 10 亿个元素）稀疏矩阵（值、行、列），我没有打开 mp 的代码如下：

function M_sparse(..) result(M)

               (variables declarations)

               cnt=0
               do i=1,n
                 do j=i,n
                   v = some_function(..)
                   if (v /= 0.) then
                     cnt=cnt+1
                     ht(cnt)=v
                     it(cnt)=dble(i)
                     jt(cnt)=dble(j)
                   endif
                 end do
               enddo

              allocate(M(cnt,3))
              M(:,1)=ht(:cnt)
              M(:,2)=it(:cnt)
              M(:,3)=jt(:cnt)
              return
end function

现在我真的很困惑如何并行化它。我至少需要对 ht、it 和 jt 进行串行更新，但到目前为止，在每次尝试中，cnt 的最终值甚至在多次运行时都不稳定。

【问题讨论】：

欢迎您，请拨打tour 并阅读How to Ask。我们可能需要更多代码。 some_function() 长什么样子？是纯的吗？您是否在尝试中标记了需要为private 的变量？最好展示您的 OpenMP 尝试、完整的可编译代码 (minimal reproducible example)。

您应该可以直接从i 和j 的值计算cnt，不是吗？ cnt=cnt+1 对并行化有问题。但是需要完整的代码。

您是否担心元素存储在 M 数组中时的最终顺序？只要它们是相同的元素，它们可以按任何顺序排列吗？即使它总是相同的元素，顺序是否会随着线程的数量而变化？你可以看看***.com/questions/68404280/…的cmets@

@VladimirF cnt 是 some_function 是否返回零的函数。

@bslhrzg 当然要保持它的可读性，但要说明函数的作用以及它是否是线程安全的。写some_function(...) 至少不提及重要属性是不够的。您不妨创建一个虚拟实现，甚至调用一个随机数生成器，但显示它是否是线程安全的以及是否可以预测零结果。请务必声明您的变量。这真的很重要。我的意思是我可能已经真正看到了 十分之一 的问题，我们必须从缺少的变量声明中获取关键信息。

【参考方案1】：

这是我将如何做的一个被破解的版本 - 它本质上是@veryreverie 建议的一个版本：生成一组线程私有列表，然后将它们连接起来。注意

我假设您不在乎元素的列出顺序。如果你现在有一个这样一个本质上不平行的问题，那么你现在有一个排序问题，这将更难以解决 无法测试其结果的程序是没有意义的 - 所以我的程序会根据单线程结果检查 2、3 和 4 线程结果。请注意，因为现在是星期五晚上，我感到非常懒惰，这种检查虽然很重要，但效率却可怕，而且实际上对于大型案例而言，所花费的时间比计算本身要长得多！

这里是代码、它是如何编译的，以及我的四核笔记本电脑上的一些示例结果：

ijb@ijb-Latitude-5410:~/work/stack$ cat listing.f90
Program listing

  Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64

  Implicit None

  Type element_type
     Integer    :: i, j
     Real( wp ) :: Hij
  End Type element_type

  Type( element_type ), Dimension( : ), Allocatable :: list_of_elements_serial
  Type( element_type ), Dimension( : ), Allocatable :: list_of_elements

  Integer :: n
  Integer :: nth

  Integer( li ) :: start, finish, rate

  Logical :: worked

  Write( *, * ) 'n ?'
  Read ( *, * )  n

  nth = 1
  Call system_clock( start, rate )
  ! On a Single thread generate a reference list to check against
  Call generate_list( n, nth, list_of_elements_serial )
  Call system_clock( finish, rate )
  Write( *, * ) 'time on ', 1, ' threads = ', Real( finish - start, wp ) / rate, Size( list_of_elements_serial )

  ! On 2, 3, 4 generate the lists, compare performance, check the results are correct
  Do nth = 2, 4
     Call system_clock( start, rate )
     Call generate_list( n, nth, list_of_elements )
     Call system_clock( finish, rate )
     Write( *, * ) 'time on ', nth, ' threads = ', Real( finish - start, wp ) / rate, Size( list_of_elements )
     Call checkit( list_of_elements_serial, list_of_elements, worked )
     Write( *, '( "Checking ... ")', Advance = 'No' )
     If( .Not. worked ) Then
        Write( *, * ) 'Failed on ', nth, Size( list_of_elements )
     Else
        Write( *, * ) 'Worked'
     End If
  End Do

Contains

  Subroutine generate_list( n, nth, list_of_elements )

    ! Generate a list of the non-zero elements

    Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64

    Use omp_lib, Only : omp_get_thread_num

    Implicit None

    Integer                                          , Intent( In    ) :: n                 ! Size of matrix
    Integer                                          , Intent( In    ) :: nth               ! number of threads
    Type( element_type ), Dimension( : ), Allocatable, Intent(   Out ) :: list_of_elements  ! The list of elements

    Real( wp ), Parameter :: tol = 1.0e-16_wp
    
    Integer, Parameter :: n_chunk = 16384

    Type( element_type ), Dimension( : ), Allocatable :: private_list
    Type( element_type ), Dimension( : ), Allocatable :: temp_list
    
    Real( wp ) :: v

    Integer, Dimension( : ), Allocatable :: counts
    
    Integer :: private_count
    Integer :: my_start
    Integer :: i, j

    Interface
       Pure Function func( n, i, j ) Result( v )
         Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64
         Real( wp ) :: v
         Integer, Intent( In ) :: n
         Integer, Intent( In ) :: i
         Integer, Intent( In ) :: j
       End Function func
    End Interface

    !$omp parallel num_threads( nth ) default( none ) &
    !$omp private( private_count, private_list, temp_list, my_start, v, i, j ) &
    !$omp shared( n, nth, counts, list_of_elements )

    ! Generate a subset of the elements local to this thread
    Allocate( private_list( 1:n_chunk ) )

    private_count = 0

    !$omp do 
    Do i = 1, n
       Do j = 1, n
          v = func( n, i, j )
          If( Abs( v ) > tol ) Then
             private_count = private_count + 1
             If( private_count > Ubound( private_list, Dim = 1 ) ) Then
                Allocate( temp_list( 1:Ubound( private_list, Dim = 1 ) + n_chunk ) )
                temp_list( 1:Ubound( private_list, Dim = 1 ) ) = private_list
                Call move_alloc( temp_list, private_list )
             End If
             private_list( private_count )%i   = i
             private_list( private_count )%j   = j
             private_list( private_count )%Hij = v
          End If
       End Do
    End Do

    ! Concatenate the private lists into one shared list
    
    !$omp single
    Allocate( counts( 0:nth - 1 ) )
    !$omp end single

    counts( omp_get_thread_num() ) = private_count
    !$omp barrier

    !$omp single
    Allocate( list_of_elements( 1:Sum( counts ) ) )
    !$omp end single

    my_start = Sum( counts( 0:omp_get_thread_num() - 1 ) ) + 1
    list_of_elements( my_start:my_start + private_count - 1 ) = private_list( 1:private_count )

    !$omp end parallel
    
  End Subroutine generate_list

  Pure Subroutine checkit( list_ref, list, worked )

    ! Check whether the given list is just a rearrangement of the reference list
    ! HORRIBLY inefficient, should really use sorting - can't be bothered.

    Implicit None
    
    Type( element_type ), Dimension( : ), Intent( In    ) :: list_ref
    Type( element_type ), Dimension( : ), Intent( In    ) :: list
    Logical                             , Intent(   Out ) :: worked

    Type( element_type ), Dimension( : ), Allocatable :: temp

    Integer :: i, j
    
    worked = .True.

    If( Size( list_ref ) /= Size( list ) ) Then
       worked = .False.
    End If
    
    Allocate( temp, Source = list )

    Do i = 1, Size( list_ref )
       Do j = 1, Size( list )
          ! Search for element i of the reference list in the list being checked
          If( list_ref( i )%i == temp( j )%i .And. &
              list_ref( i )%j == temp( j )%j .And. &
              Abs( list_ref( i )%Hij - temp( j )%Hij ) < 1e-15_wp ) Then
             Exit
          End If
       End Do
       If( j == Size( list ) + 1 ) Then
          worked = .False.
          Return
       End If
       ! Mark it as used already
       temp( j )%i   = -1
       temp( j )%j   = -1
       temp( j )%Hij = Huge( temp( j )%Hij )
    End Do

  End Subroutine checkit
    
End Program listing

Pure Function func( n, i, j ) Result( v )

  ! silly function for sparse matrix
  
  Use, Intrinsic :: iso_fortran_env, Only : wp => real64, li => int64

  Real( wp ) :: v

  Integer, Intent( In ) :: n
  Integer, Intent( In ) :: i
  Integer, Intent( In ) :: j

  If( 100 * i < n .And. 100 * j < n ) Then
     v = 1.0_wp
  Else
     v = 0.0_wp
  End If
  
End Function func
ijb@ijb-Latitude-5410:~/work/stack$ gfortran-11 --version
GNU Fortran (GCC) 11.1.0
Copyright © 2021 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

ijb@ijb-Latitude-5410:~/work/stack$ gfortran-11 -std=f2008 -Wall -Wextra -O3 -g -fopenmp  listing.f90 -o gen_list
ijb@ijb-Latitude-5410:~/work/stack$ ./gen_list
 n ?
10000
 time on            1  threads =    6.7302687000000000E-002        9801
 time on            2  threads =    2.6817233999999999E-002        9801
Checking ...  Worked
 time on            3  threads =    1.5919547999999999E-002        9801
Checking ...  Worked
 time on            4  threads =    1.1952938000000000E-002        9801
Checking ...  Worked
ijb@ijb-Latitude-5410:~/work/stack$ ./gen_list
 n ?
30000
 time on            1  threads =   0.44568265400000001            89401
 time on            2  threads =   0.21186449299999999            89401
Checking ...  Worked
 time on            3  threads =   0.14133034500000000            89401
Checking ...  Worked
 time on            4  threads =   0.12390519100000000            89401
Checking ...  Worked
ijb@ijb-Latitude-5410:~/work/stack$ ./gen_list
 n ?
60000
 time on            1  threads =    1.7274770189999999           358801
 time on            2  threads =   0.85456061200000000           358801
Checking ...  Worked
 time on            3  threads =   0.57058082499999996           358801
Checking ...  Worked
 time on            4  threads =   0.42949695500000001           358801
Checking ...  Worked
ijb@ijb-Latitude-5410:~/work/stack$ 

【讨论】：

非常感谢，它适用于我的情况，现在我需要再次仔细阅读以确保理解。对于具有 100 万个非零元素的测试用例，我在笔记本电脑上得到 9s（6 个线程）而不是 25s（单线程）。 （我想现在在可分配数组上使用追加必须减慢一点进程，但优点是它对内存限制更有弹性），再次感谢！

【参考方案2】：

另一个想法：将密集数组分成块，每个线程负责一个块。让每个线程从其自己的密集数组部分生成稀疏数组的一部分，然后在必要时将这些部分连接在一起。

【讨论】：

做到了，我需要帮助来实现它，谢谢

【参考方案3】：

这里有一个解决方案：创建一个矩阵大小的数组，计算所有i,j 的some_function 并记录有多少js 给出非零结果，将其存储在位置i。这是完全平行的。

现在您知道需要多少空间，并且可以为每个线程指定其在存储中的起点。再次遍历some_function 并实际填充元素，cnt 是每个线程本地的。

好的，所以这使标量工作量加倍。但是你让它完全平行，所以你真的不在乎，对吧？

【讨论】：

感谢您的提议，我会尝试一些不同的方法（使用非零元素创建 i,j 的映射），看看它是否比下面的解决方案更快（使用不可分配的数组)，但我喜欢一开始不必声明巨大的数组，因为我也很快面临内存问题

不确定您所说的“i,j 地图”是什么意思。这就是整个问题：您不知道需要存储多少 i,j 元素。但是你知道有多少 i，所以你可以用它来计算每个 i 的 j。这里有一个想法：也许你可以写一个廉价版本的some_function 来确定 i,j 元素是否非零，但不经过全值计算。