# Distributed (MPI) matrix matrix multiplication

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 PDGEMM or PETSc's 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

• How large are n1, n2, n3? Having in mind, that for most OpenMP accelerated BLAS libraries the number of threads for a single BLAS call are reduced to 1 if they are called from an OpenMP parallelized region. For this reason it is only worth using your OpenMP idea if n1 and n2 are relatively small. – M.K. aka Grisu Apr 7 '17 at 13:56
• At this moment they are pretty small such as (32,32,32). I tested OpenMP way with the sizes around (100,100,100), it is okay in terms of efficiency (but limited with the node size). To scale, I am using a kind of domain decomposition(DD), so that part is on MPI. If I can have MPI version of the routines above, I plan to replace DD with it, and dreaming to solve full size of approximately (6000, 500, 500), via few thousand cores. – trblnc Apr 7 '17 at 14:16
• I read your comment again @M.K.akaGrisu, now I understand what do you mean, and I agree. I will reshape my question accordingly. – trblnc Apr 7 '17 at 16:42

1. The Cannon algorithm is pretty cute if you're on a square processor grid. In each step you rotate the input matrix rows and columns, so that in the end each processor contains a sum of A(i,k)B(k,j) products.