我使用以下子程序在 fortran 中执行矩阵矩阵乘法（在 rank-3 和 rank-2 数组之间），

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

我喜欢它的 MPI 版本。我知道它的效率可能会降低。

如何以一种简单而干净的方式做到这一点？

我考虑使用PDGEMMor （可能也在使用）函数PETSc的选项。我已经在相同的代码中使用 PETSC 来求解线性系统，而且我发现一开始更难掌握。matmatmultPDGEMMPDGEMM

还有其他（更好的）方法吗？

使用 OpenMP 加速的 BLAS 运行良好。

或者，也可以使用以下子程序，请参阅MK 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