I perform matrix matrix multiplications (between rank-3 and rank-2 arrays) in fortran using following subroutine,
subroutine ab_c_dir1(n1,n2,n3,a,b,c) implicit none integer :: n1,n2,n3 real(rk) :: a(n1,n1),b(n1,n2,n3),c(n1,n2,n3) call dgemm('n','n',n1,n2*n3,n1,1.d0,a,n1,b,n1,0.d0,c,n1) end subroutine ab_c_dir1
I like to have MPI version of it. I understand that it can became less efficient.
How this can be done in an easy and clean way?
I think about the options using
matmatmult(may be it is using
PDGEMMas well) function. I am already using PETSC in the same code to solve a linear system, also I find
PDGEMM harder to grasp at first place.
Is there any other (better) approach?
Using OpenMP accelerated BLAS is working quite well.
Alternatively following subroutine can be used as well, please see the warning in the comment section from M.K. aka Grisu.
subroutine ab_c_dir1_openMP(n1,n2,n3,a,b,c) !$ use OMP_LIB implicit none integer :: n1,n2,n3 real(rk) :: a(n1,n1),b(n1,n2,n3),c(n1,n2,n3) integer :: k !$OMP PARALLEL PRIVATE(k) !$OMP DO SCHEDULE(RUNTIME) do k=1,n3 call dgemm('n','n',n1,n2,n1,1.d0,a,n1,b(1,1,k),n1,0.d0,c(1,1,k),n1) enddo !$OMP END DO !$OMP END PARALLEL end subroutine ab_c_dir1_openMP