Sparse matrix parallel creation with openmp in fortran

user_17644599 0
Viewed 150

I am relatively new to fortran, and totally new to openmp, I have the following problem :

I want to construct a (big : ~1% non zero elements over ~1 million to 1 billion elements in total) sparse matrix (value,row,column) in parallel, my code without open mp is the following :

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

Now I am really confused about how to parallelize this. I need at least the updates on ht,it and jt done in serial, but in every attempt so far, the final value of cnt is not even stable against multiples runs.

openmp fortran
3 Answers

Here is one solution: make an array of the matrix size, compute the some_function for all i,j and record how many js give a nonzero result, store that in location i. This is perfectly parallel.

Now you know how much space you need, and you can give each thread its starting point in the storage. Go through the some_function again and actually fill in the elements, with cnt local to each thread.

Ok, so this doubles the amount of scalar work. But you make it completely parallel, so you don't really care, right?

user_20444540

Another idea: split your dense array into chunks, with each thread responsible for one chunk. Have each thread generate a section of the sparse array from its own section of the dense array, and then concatenate the sections together at the end if necessary.

user_85593000

Here's a hacked together version of how I would do it - it's essentially a version of what @veryreverie suggests: generate a set of thread private lists and then concatenate them. Note

  1. I have assumed you don't care what order the elements are listed. If you do you now have an ordering an as such an inherently non-parallel problem, which will be somewhat more difficult to solve
  2. A program which can't test its results is pointless - so my program checks the 2, 3 and 4 thread results against the single threaded ones. Note because it is Friday evening and I am feeling very lazy this checking, though vital, is done horribly inefficiently, and in fact for big cases takes much longer than the calculation itself!

Anyway here is the code, how it is compiled, and some example results on my quad core laptop:

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$
user_12804390
Related