2
$\begingroup$

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
Improve this question
$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$
    – M.K. aka Grisu
    Apr 7, 2017 at 13:56
  • $\begingroup$ 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. $\endgroup$
    – trblnc
    Apr 7, 2017 at 14:16
  • $\begingroup$ I read your comment again @M.K.akaGrisu, now I understand what do you mean, and I agree. I will reshape my question accordingly. $\endgroup$
    – trblnc
    Apr 7, 2017 at 16:42

1 Answer 1

Reset to default
5
$\begingroup$

You say that you want an MPI version. Then you need to study the literature, as the distributed memory variant of matrix-matrix product are not a simple parallellization of the sequential version.

  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.

  2. The Summa algorithm is more or less a bunch of rank k updates. This is more general in shapes of the matrices and processor grids.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.