在 Fortran 中的数组子集上使用 BLAS ?gemm

Question

BLAS ?gemm 函数的各种 LDx 参数可以对较大数组的切片进行操作。例如，这个小型 C 程序对 (200,200) 矩阵的左上角和右上角 (100,100) 子矩阵进行矩阵乘法，并将结果存储在左下角 (100,100) 子矩阵中。

#include <stdlib.h>
#include <stdio.h>
#include <cblas.h>
int main()
{
    double * a = malloc(200*200*sizeof(double));
    for(int i = 0; i < 200*200; i++) a[i] = 1;
    cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, 100, 100, 100, 1, a, 200, a+100*200, 200, 0, a+100, 200);
    printf("%f\n", a[100]);
    return 0;
}

$ gcc -o cblas_test{,.c} -lcblas
$ ./cblas_test
100.000000

输出为 100，符合预期。由于 BLAS 最初是一个 fortran 库，我希望这也可以在那里实现。但我看不到如何在那里指定数组偏移量。使用切片不起作用，导致以下程序出现分段错误：

program fblas_example
    implicit none
    real(8), allocatable :: a(:,:)
    allocate(a(200,200))
    a = 1
    call dgemm('n','n',100,100,100,1d0,a,200,a(1:100,101:200),200,0d0,a(101:200,1:100),200)
    write(*,*) a(101,1)
end program

$ gfortran -o fblas_example{,.f90} -lblas
$ ./fblas_example
*** Error in `./fblas_example': double free or corruption (out): 0x00000000011351c0 ***

我做错了什么？

编辑：如果我先将 a 按摩到一维数组中，这将起作用：

program fblas_example2
    use, intrinsic :: iso_c_binding
    implicit none
    real(8), target, allocatable :: a(:,:)
    real(8), pointer :: fa(:)
    allocate(a(200,200))
    a = 1
    call c_f_pointer(c_loc(a), fa, [size(a)])
    call dgemm('n','n',100,100,100,1d0,fa,200,fa(100*200+1:),200,0d0,fa(101:),200)
    write(*,*) fa(101)
end program

$ gfortran -o fblas_example2{,.f90} -lblas
$ ./fblas_example2
   100.00000000000000

但奇怪的是，必须通过 ISO C 绑定才能从 fortran 调用 fortran 库。

明确一点：我知道可以复制块，处理它们，然后将它们复制回来。但我在这里要问的是如何使用 BLAS 自己的支持来处理数组的子集。

Answer 1

当您给出 LDx 时，您声明了矩阵的实际大小（前导维度）。然后BLAS使用这个值从乘法中跳过未使用的数据。

我认为您应该像 C 一样使用 dgemm。我的意思是你不应该传递子矩阵 a(1:100,101:200) 而是传递子矩阵的第一个值 a(1,101) 的位置。

我测试了以下内容，似乎工作正常。

program fblas_example
  implicit none
  real(8), allocatable :: a(:,:)
  allocate(a(200,200))
  a = 1
  call dgemm('n','n',100,100,100,1d0,a,200,a(1,101),200,0d0,a(101,1),200)
  write(*,*) a(101,1)
end program

$ gfortran dgemm.f -lblas
$ ./a.out 
  100.00000000000000

Answer 2

如果要传递 Fortran >=90 中的数组切片（而不是某个数组元素的地址），我认为 Fortran 代码中的以下行

call dgemm('n','n', 100,100,100, 1d0, a,200, a(1:100,101:200),200, 0d0, a(101:200,1:100),200 )

应该修改为

call dgemm('n','n', 100,100,100, 1d0, a,200, a(1:100,101:200),100, 0d0, a(101:200,1:100),100 )

其中后两个数组切片（或子数组）的前导维度已从 200 更改为 100。

更具体地说，a(1:100,101:200) 和 a(101:200,1:100) 通过隐式接口传递给 dgemm，其中数组切片在内存中是不连续的。在这种情况下，编译器首先会为每个子数组准备两个大小为 100x100 的临时数组，将原始数组 a 的矩阵元素复制到临时数组中，并传递临时数组的第一个元素的地址。一旦 dgemm 的计算完成，临时对象的内容将适当地复制回原始数组 a。如果发生这种情况，dgemm 会收到两个前导维度为 100（而不是 200）的数组。

可以通过 -fcheck-array-temporaries（在 gfortran 中）和 -check arg_temp_created（在 ifort 中）选项检查临时数组的生成。例如，后者给出

forrtl: warning (402): fort: (1): In call to DGEMM, an array temporary was created for argument #9
forrtl: warning (402): fort: (1): In call to DGEMM, an array temporary was created for argument #12
100.000000000000

在 Fortran 中的数组子集上使用 BLAS ?gemm

Using BLAS ?gemm on a subset of an array in fortran

arrays

fortran

slice

blas