One of the core operations of my code is sum of N velocity fields, each one weighted by a gaussian random variable, into a single velocity field. In LaTeX notation it would be
U(x,y,z) = \sum_{n=1}^{N} \phi_{n}(x,y,z) \xi_{n}
where \phi are the velocity fields and U is the result. I am currently considering different options of implementation to gain as much speed-up as possible, since my matrix can be really large: each field \phi has dimension (100,100,30), and they are in total 720.
I am currently investigating two ways of storing this data for an efficient use:
- The first has all the \phi fields piled one over the other, meaning in a single big matrix of dimension (100,100,30 * 720)
- The second is reshaping each \phi into a vector, and concatenate them to form a matrix of dimension (100 * 100 * 30,720)
The convenience of the second one is that once I have my random vector \xi, I can use BLAS dgemv to compute the matrix-vector product, and results show that the second option is more efficient.
However, since we are talking about a matrix of dimension (300.000,720), I thought that using OpenMp would bring a little benefit. I have set up a manual scheduling, where I divide the matrix into sub matrices, use private pointers to point at the portion of the matrix and call dgemv for each sub matrix. The result is that it takes more time than the original BLAS operation on the full matrix, even with 2 threads (which should not add too much communication overhead) so I am clearly missing some point here that I don't see.
I am compiling with a simple
gfortran -o -02 main.out MWE.F90 -llapack -lblas -fopenmp
and the code I am using is:
PROGRAM new_modes
IMPLICIT NONE
INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307)
interface
function OMP_get_wtime()
INTEGER, PARAMETER :: wp = SELECTED_REAL_KIND(12,307)
real(wp) :: OMP_get_wtime
end function OMP_get_wtime
function OMP_get_thread_num()
integer(kind = 4) :: OMP_get_thread_num
end function OMP_get_thread_num
end interface
INTEGER, PARAMETER :: nn_tlu_nmod = 720
INTEGER, PARAMETER :: jpi = 100! 32
INTEGER, PARAMETER :: jpj = 100! 17
INTEGER, PARAMETER :: jpk = 30! 31
INTEGER, PARAMETER :: np = jpk * jpj * jpi
!
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: A
REAL(wp), ALLOCATABLE, TARGET, DIMENSION(:,:) :: B
REAL(wp), POINTER, DIMENSION(:,:) :: sub_B
REAL(wp), POINTER, DIMENSION(:) :: sub_B_vec
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: U_ref
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: U
REAL(wp), ALLOCATABLE, TARGET, DIMENSION(:) :: vec_U
REAL(wp), POINTER, DIMENSION(:) :: sub_vec_U
REAL(wp), ALLOCATABLE, TARGET, DIMENSION(:) :: tcoeff
REAL(wp), POINTER, DIMENSION(:) :: tcoef_pt
!
REAL(wp) :: tic, toc
INTEGER :: ji, jj, jk, jm
INTEGER :: m_idx
INTEGER :: cnt
!
INTEGER :: ierr
!
REAL(wp) :: ddot
!
!-------------------------------------------------------------------------------------
! Declaration for the Parallel parameters
!-------------------------------------------------------------------------------------
INTEGER, PARAMETER :: num_threads=1
INTEGER, SAVE :: myID
!$OMP THREADPRIVATE(myID)
!-------------------------------------------------------------------------------------
! Declaration for the row scheduling
!-------------------------------------------------------------------------------------
INTEGER :: sub_n
INTEGER :: istart,inext
INTEGER, ALLOCATABLE :: work_sharing(:)
WRITE(6,*) 'Number of threads employed in the solver:', num_threads
ALLOCATE( A(jpi,jpj,jpk * nn_tlu_nmod), stat=ierr )
ALLOCATE( B(jpi*jpj*jpk, nn_tlu_nmod), stat=ierr )
ALLOCATE( Bt( nn_tlu_nmod, jpi*jpj*jpk), stat=ierr )
ALLOCATE( tcoeff(nn_tlu_nmod), stat=ierr )
ALLOCATE( U(jpi,jpj,jpk), &
& U_ref(jpi,jpj,jpk), &
& vec_U(jpi*jpj*jpk), stat=ierr )
cnt = 0
DO jm = 1, nn_tlu_nmod
!
m_idx = ( jm - 1 ) * jpk
!
DO jk = 1, jpk
DO jj = 1, jpj
DO ji = 1, jpi
cnt = cnt + 1
A(ji,jj,m_idx + jk) = cnt * jm
END DO
END DO
END DO
B(:,jm) = RESHAPE(A(:,:,m_idx + 1: m_idx + jpk), (/jpi*jpj*jpk/) )
Bt(jm,:) = RESHAPE(A(:,:,m_idx + 1: m_idx + jpk), (/jpi*jpj*jpk/) )
END DO
tcoeff = 1._wp
! First option, normal matrix summation (already tested faster than
! with explicit do loop)
tic = omp_get_wtime()
U = 0._wp
!
DO jm = 1, nn_tlu_nmod
! Define zero-th indexed mode
m_idx = ( jm - 1 ) * jpk
!
U(:,:,:) = U(:,:,:) + A(:,:,m_idx + 1 : m_idx + jpk ) * tcoeff(jm)
!
ENDDO
!
toc = omp_get_wtime()
print '("Time, operation without for loops = ",f13.7," seconds.")', toc-tic
U_ref = U
! Second option, exploit a different shape of the matrix in order to use BLAS
! (Usage of pointers speeds-up a lot)
tic = omp_get_wtime()
U = 0._wp
sub_B => B(:, :)
sub_n = np
sub_vec_U => vec_U(:)
tcoef_pt => tcoeff(:)
! CALL DGEMV( trans, m, n, alpha, A, ldA, x, incx, beta, y, incy )
CALL DGEMV( 'n', sub_n, nn_tlu_nmod, 1._wp, sub_B, sub_n, tcoef_pt, 1, 0._wp, sub_vec_U, 1 )
U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) )
!
toc = omp_get_wtime()
print '("Time, 2D matrix with BLAS = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref))
! Third option, Divide the job in a static way between M threads, use pointers to point at the
! portion of the matrix assigned to each processor and use BLAS
tic = omp_get_wtime()
vec_U = 0._wp
!$OMP PARALLEL PRIVATE(istart, sub_n)
myID = omp_get_thread_num()+1
istart = work_sharing(myID)
sub_n = work_sharing(myID+1)-istart
! CALL DGEMV( trans, m, n, alpha, A, ldA, x, incx, beta, y, incy )
CALL DGEMV( 'n', sub_n, nn_tlu_nmod, 1._wp, B(istart,1), np, tcoeff, 1, 0._wp, vec_U(istart), 1 )
!$OMP END PARALLEL
U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) )
toc = omp_get_wtime()
print '("Time, Static scheduling BLAS OpenMP = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref))
! Fourth option, run a parallel Do and use DDOT
tic = omp_get_wtime()
vec_U = 0._wp
!$OMP PARALLEL PRIVATE(ji)
!$OMP DO SCHEDULE (STATIC)
DO ji=1,np
! res = DDOT( n, x, inc_x, y, inc_y )
vec_U(ji) = DDOT( nn_tlu_nmod, B(ji, :), 1, tcoeff, 1 )
END DO
!$OMP END DO
!$OMP END PARALLEL
U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) )
toc = omp_get_wtime()
print '("Time, Parallelized for loop DDOT = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref))
! Fifth option, Divide the job in a static way between M threads, but using a transposed matrix
tic = omp_get_wtime()
vec_U = 0._wp
!$OMP PARALLEL PRIVATE(istart, sub_n)
myID = omp_get_thread_num()+1
istart = work_sharing(myID)
sub_n = work_sharing(myID+1) - istart
! CALL DGEMV( trans, m, n, alpha, A, ldA, x, incx, beta, y, incy )
CALL DGEMV( 't', nn_tlu_nmod, sub_n, 1._wp, Bt(1,istart), nn_tlu_nmod, tcoeff, 1, 0._wp, vec_U(istart), 1 )
!$OMP END PARALLEL
U = RESHAPE(vec_U, (/ jpi, jpj, jpk /) )
toc = omp_get_wtime()
print '("Time, Static scheduling tran OpenMP = ",f13.7," seconds.", f25.7)', toc-tic, MAXVAL(ABS(U-U_ref))
END PROGRAM
Thanks to all
EDIT: I have modified the initial code with all the suggestions of PierU, and added a fifth possibility, that basically transpose the matrix from the start (feasible in my case), in order to do a scheduling of contiguous memory positions. For some reasons this last one is anyway slower than the third one, which is counter-intuitive for me.