使用 DGEMM 的 BLAS LDB

Question

我想将矩阵乘以 D*W'，其中 W' 是 W 的转置版本。

虽然我将使用 DGEMM，但我在 @IanBush 的帮助下发现，这种情况下的 LDB 应该是矩阵 W 的行数而不是列数。本例的代码是

  Call dgemm('n', 't', N1, M1, N1, 1.0_wp, D, N1, W, M1, 0.0_wp, c, n1)

其中 n1 和 m1 是我的矩阵的维度

Dimensions of the matrices:
W = M1*N1
D = N1*N1

正如官方文档中所说

      LDB is INTEGER
       On entry, LDB specifies the first dimension of B as declared
       in the calling (sub) program. When  TRANSB = 'N' or 'n' then
       LDB must be at least  max( 1, k ), otherwise  LDB must be at
       least  max( 1, n ).

为什么会出现这种情况？

Answer 1

对 dgemm 的调用包含两组整数。第一个是定义数学描述的矩阵的维度。查看http://www.netlib.org/lapack/explore-html/d7/d2b/dgemm_8f.html上的 BLAS 文档，我们可以看到例程的参数定义为（并以有助于我的方式格式化它记住它是如何工作的）：

         Subroutine dgemm (TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, &
                                                           B, LDB, & 
                                                     BETA, C, LDC)

定义数学维度的整数集M, N, K在该集中

• M 是结果矩阵的第一个数学维度，C
• N 是结果矩阵的第二个数学维度C
• K 是生成结果矩阵的单个元素所需的点积的长度C

所以M, N, K纯粹与问题的数学有关，没有别的。但我们是在计算机上，所以有第二个问题——每个二维矩阵在计算机内存中是如何布局的，而计算机内存本质上是一维对象？现在，Fortran 按列主顺序存储其对象。这意味着在内存中，如果您有元素 A( i, j )，则下一个内存位置将保存 A( i + 1, j ) 如果我们不是列的末尾，或者 A( 1, j + 1 ) 如果我们不是列的末尾我们是。因此，矩阵的布局是“第一个索引移动最快”。要定义二维对象的布局，我们需要告诉程序的是您需要跳过 A 的多少个元素才能从 A( i, j ) 到 A( i, j + 1 ) - 就是这样数字是“主维度”，LDA, LDB, LDC文档中的因此，它只是声明或分配的矩阵中的行数。

现在让我们看看您引用的文档

  LDB is INTEGER
   On entry, LDB specifies the first dimension of B as declared
   in the calling (sub) program. When  TRANSB = 'N' or 'n' then
   LDB must be at least  max( 1, k ), otherwise  LDB must be at
   least  max( 1, n ).

LDB 是矩阵 B 中声明的行数。因此

1. 如果 B 未转置，B 的每一列将参与与 A 的行或列的点积，具体取决于 A 的转置选项。数学维度表示点积的长度为 k。因此，B 中的行数必须至少为 k 才能使数学正确，因此 LDB 必须至少为 k
2. 如果 B 转置，B 的数学第一维也将是结果矩阵的数学第二维，并且由 N 给出。因此在这种情况下 LDB 必须至少为 N

因此，这回答了大部分问题，如果您总是将声明的一个矩阵的所有内容与第二个矩阵的所有内容相乘，则前导维度将只是声明的每个矩阵的第一个维度。但“至少”意味着您可以将数组的位相乘。考虑以下内容，并通过仔细区分定义数学的整数和定义内存布局的整数，您能弄清楚为什么它会这样做吗？请注意，参数列表中的 a( 2, 2 ) 意味着我们从该元素“开始”矩阵 - 现在仔细思考主要维度告诉您什么，您应该能够弄清楚它是如何工作的。

ian@eris:~/work/stack$ cat matmul2.f90
Program sirt

  Integer, Parameter :: wp = Kind( 1.0d0 )

  Real( wp ), Dimension( 1:5, 1:5 ) :: A
  Real( wp ), Dimension( 1:3, 1:3 ) :: B
  Real( wp ), Dimension( 1:4, 1:4 ) :: C

  Integer :: i

  A = 0.0_wp
  Do i = 1, 5
     A( i, i ) = Real( i, wp )
  End Do
  Write( *, * ) 'A = '
  Do i = 1, Size( A, Dim = 1 )
     Write( *, '( 100( f4.1, 1x ) )' ) A( i, : )
  End Do

  B = 0.0_wp
  B( 1, 1 ) = 1.0_wp
  B( 3, 2 ) = 1.0_wp
  B( 2, 3 ) = 1.0_wp
  Write( *, * ) 'B = '
  Do i = 1, Size( B, Dim = 1 )
     Write( *, '( 100( f4.1, 1x ) )' ) B( i, : )
  End Do

  ! Lazy - should really just initialise only the bits of C that are NOT touched
  ! by the dgemm
  C = 0.0_wp

  Call dgemm('N', 'N', 3, 3, 3, 1.0_wp, A( 2, 2 ), Size( A, Dim = 1 ), &
                                        B        , Size( B, Dim = 1 ), &
                                0.0_wp, C        , Size( C, Dim = 1 ) )

  Write( *, * ) 'C after dgemm'
  Write( *, * ) 'B = '
  Do i = 1, Size( C, Dim = 1 )
     Write( *, '( 100( f4.1, 1x ) )' ) C( i, : )
  End Do


End Program sirt
ian@eris:~/work/stack$ gfortran-8  -std=f2008 -fcheck=all -Wall -Wextra -pedantic -O matmul2.f90 -lblas
/usr/bin/ld: warning: libgfortran.so.4, needed by //usr/lib/x86_64-linux-gnu/libopenblas.so.0, may conflict with libgfortran.so.5
ian@eris:~/work/stack$ ./a.out
 A = 
 1.0  0.0  0.0  0.0  0.0
 0.0  2.0  0.0  0.0  0.0
 0.0  0.0  3.0  0.0  0.0
 0.0  0.0  0.0  4.0  0.0
 0.0  0.0  0.0  0.0  5.0
 B = 
 1.0  0.0  0.0
 0.0  0.0  1.0
 0.0  1.0  0.0
 C after dgemm
 B = 
 2.0  0.0  0.0  0.0
 0.0  0.0  3.0  0.0
 0.0  4.0  0.0  0.0
 0.0  0.0  0.0  0.0

这种将两个较大矩阵的一部分相乘的用法非常常见，尤其是在块线性代数算法中，并且数学维度和布局维度的分离使您可以做到这一点