使用 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

假设你有一个大矩阵 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 就知道它必须跳过内存块的哪些部分。

Answer 2

对 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 ) 如果我们是。所以矩阵的布局是 'first index moving fastest'。要定义二维对象的这种布局，我们只需要告诉程序您需要跳过 A 的多少个元素才能从 A( i, j ) 到 A( i, j + 1 ) - 就是这样number 即文档中的 'leading dimension'、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

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

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

这种将两个较大矩阵的一部分相乘的用法出奇地普遍，尤其是在分块线性代数算法中，数学维度和布局维度的分离让您可以做到这一点

使用 DGEMM 的 BLAS LDB

BLAS LDB using DGEMM

matrix

fortran

hpc

lapack

blas