arrays - 在 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 绑定(bind)才能从 fortran 调用 fortran 库。

明确一点:我知道人们可以复制 block ，处理它们，然后将它们复制回来。但我在这里要问的是如何使用 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