使用 DGEMM 的 BLAS LDB答案

【问题标题】：BLAS LDB using DGEMM使用 DGEMM 的 BLAS LDB
【发布时间】：2020-03-31 22:33:52
【问题描述】：

我想将矩阵乘以 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 ).

为什么会这样？

【问题讨论】：

标签： matrix fortran hpc lapack blas

【解决方案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 以 Column Major 顺序存储其对象。这意味着如果您有元素 A( i, j ) 在内存中，如果我们不是列的末尾，则下一个内存位置将保存 A( i + 1, j )，或者如果我们不是列的末尾，则为 A( 1, j + 1 )我们是。因此，矩阵布局为“第一个索引移动最快”。要定义二维对象的这种布局，我们只需要告诉程序从 A(i, j) 到 A(i, j + 1) 需要跳过 A 的多少元素 - 就是这样数字，这是文档中的“领先维度”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 中声明的行数。因此

如果 B 未转置，则 B 的每一列都将包含在具有 A 行或列的点积中，具体取决于 A 的转置选项。数学维度表示点积的长度为 k。因此 B 中的行数必须至少为 k 才能使数学正确，因此 LDB 必须至少为 k
如果 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

这种将两个较大矩阵的一部分相乘的用法非常普遍，尤其是在块状线性代数算法中，并且数学和布局维度的分离允许您这样做

【讨论】：

【解决方案2】：

假设你有一个大矩阵B，看起来像

b b b b B B B
b b b b B B B
b b b b B B B
B B B B B B B
B B B B B B B
B B B B B B B

B 是一个6x7 矩阵。但是在您的程序中，您不使用 B 作为矩阵，而只是作为一些存储位置。您感兴趣的真正矩阵是 b，它位于 B 的第一部分。现在，由于您将它存储在表示 B 的二维数组中，因此这里有一种内存布局形式。矩阵 b 在内存中不是连续的，因为矩阵 B 是。在记忆中，B看起来像：

b b b B B B b b b B B B b b b B B B b b b B B B B B B B B B B B B B B B B B B B B B 
  col 1    |  col 2    |  col 3    |  col 4    |  col 5    |  col 6    |  col 7

如果要将矩阵 b 传递给 LAPACK 或 BLAS，则需要告知其内存布局。这是使用大小N 和M 以及前导维度LDA 来完成的。前导维度LDA 是存储它的较大矩阵的总行数。因此，在这里您将传递N=3、M=4 和LDA=6。这样，BLAS 和 LAPACK 就知道它必须跳过内存块的哪些部分。

【讨论】：