Fortran递归与GEMM

Question

我正在尝试在Fortran中实现CARMA GEMM算法，该算法在本文CARMA Publication中找到。可以在以下位置找到用C编写的Github上源代码的链接：CARMA CILK SGEMM

该算法使用分而治之的方法，将m，n，k中最大的矩阵维数分成两个子问题，通过递归并行执行。结束条件对原始矩阵的一部分执行GEMM。这段代码使用了Cilk线程，但我对并行性还不是很感兴趣。

我在理解如何在Fortran中使用递归实现就地矩阵乘法时遇到了一个问题。下面的代码是我仅拆分维度M的尝试，它仍处于测试阶段。代码可以使用ifort 16.0进行编译，但不能正确执行。在所提供的示例中，来自主程序的C1和C2将具有不同的结果。不会出现编译器警告或运行时错误。只是错误的最终答案。

! ******************************
!         MODULE 
! ******************************
MODULE RECURSION_GEMM_MOD

IMPLICIT NONE
CONTAINS

RECURSIVE SUBROUTINE GEMM_RECURSION(m,n,k,alpha,beta,A,B,C,depth,max_depth)

integer, intent(in) :: m, n, k, depth, max_depth
real, intent(in) :: alpha, beta
real, intent(in) :: A(m,k), B(k,n)
real, intent(inout) :: C(m,n)

integer :: m_half, m_remain

! End condition for recursion - compute C <- alpha*AB + beta*C
IF( depth >= max_depth ) THEN
  CALL SGEMM('N','N', m, n, k, alpha, A, m, B, k, beta, C, m)
  RETURN
END IF

m_half = m/2
m_remain = m - m_half

! Divide GEMM into two block matrix sub-problems (top and bot)

! A_top, C_top
CALL GEMM_RECURSION(m_half,n,k,alpha,beta, A,B,C,depth+1,max_depth)

! A_bot, C_bot
CALL GEMM_RECURSION(m_remain,n,k,alpha,beta,A(m_half+1,1),B,C(m_half+1,1),depth+1,max_depth)

END SUBROUTINE GEMM_RECURSION

END MODULE RECURSION_GEMM_MOD


! **********************************
!         MAIN PROGRAM 
! **********************************
PROGRAM RECURSION_GEMM

USE RECURSION_GEMM_MOD
IMPLICIT NONE

INTEGER :: m,n,k,m_half, m_remain
INTEGER :: depth, max_depth
REAL :: C1(2,4), C2(2,4)

REAL :: A(2,3) = TRANSPOSE( RESHAPE( (/ 1.0, 2.0, 3.0,&
                                    4.0, 5.0, 6.0 /), (/ 3,2 /) ) )
REAL :: B(3,4) =  TRANSPOSE( RESHAPE( (/ 5.0, 6.0, 7.0, 8.0, &
                                     1.0, 2.0, 3.0, 4.0, &
                                     9.0, 10.0, 11.0, 12.0/), (/ 4,3 /) ) )
REAL :: C(2,4) =  TRANSPOSE( RESHAPE( (/ 1.0, 0.0, 1.0, 1.0, &
                                     0.0, 1.0, 0.0, 1.0 /), (/ 4,2 /) ) )            
REAL :: alpha = 1.0, beta = 1.0

m = 2
k = 3
n = 4

! Copy C into C1 for later
C1 = C

depth = 0
max_depth = 1

! Recursion GEMM stored in C1
CALL GEMM_RECURSION(m,n,k,alpha,beta,A,B,C1,depth,max_depth)

! Reference GEMM stored in C2
C2 = alpha*MATMUL(A,B) + beta*C

PRINT*, 'C1 = '
PRINT*, C1
PRINT*, ''

PRINT*, 'C2 = '
PRINT*, C2

END PROGRAM RECURSION_GEMM

我一直在查找有关Fortran子例程和递归的每个文档，但没有那么多。我已经检查了使用指针，但我看到的一切都表明它们的性能很糟糕，这对GEMM来说并不是很好。任何帮助都是非常感谢的。

Answer 1

问题是您通过引用传递m_half、m_remain、depth和max_depth，然后在例程中对它们进行赋值。这会更改例程的每个实例中的值。

修复方法是将VALUE属性添加到这些参数的声明中。这将传递该值的一个可写的匿名副本。