时间戳

首页 > 解决方案 > 使用 DGEMM 的 BLAS LDB

matrix - 使用 DGEMM 的 BLAS LDB

问题描述

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

为什么会这样？

标签： matrixfortranhpclapackblas

解决方案

对 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 中声明的行数。因此

  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

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

推荐阅读