SGEMM 的意外结果

Question

我正在尝试在计算 A*B' 的代码中替换以下 matmul 实例：

print *, mat%A
print *, mat%B
print *, mat%dimsA
print *, mat%dimsB
C = matmul(mat%A,transpose(mat%B))
print *, C

使用 SGEMM：

print *, mat%A
print *, mat%B
print *, mat%dimsA
print *, mat%dimsB
call SGEMM('N', 'T', mat%dimsA(1), mat%dimsB(1), mat%dimsA(2), 1, mat%A, &
    mat%dimsA(1), mat%B, mat%dimsB(1), 0, C, mat%dimsA(1))
print *, C

其中 mat 是包含以下相关数据的类型：

integer :: dimsA(2), dimsB(2)
real, pointer :: A(:,:),B(:,:)

运行小测试的代码，其中 A 是：

1 1 1 
2 2 2

B 是：

1 1 1

我得到以下 matmul 输出：

1.00000000       2.00000000       1.00000000       2.00000000       1.00000000       2.00000000    
1.00000000       1.00000000       1.00000000    
2           3
1           3
3.00000000       6.00000000

和 SGEMM：

1.00000000       2.00000000       1.00000000       2.00000000       1.00000000       2.00000000    
1.00000000       1.00000000       1.00000000    
2           3
1           3
4.20389539E-45   8.40779079E-45

我假设这两个代码块应该有相同的结果是不是错了？

Answer 1

这是因为您将 alpha 和 beta 的整数传递给 lapack 子例程。编译器没有可用的 SGEMM 接口，因此无法警告您由此产生的类型混淆...

只需使用（默认精度）实数的表示法 - 1. 而不是 1，0. 而不是 0：

call SGEMM('N', 'T', mat%dimsA(1), mat%dimsB(1), mat%dimsA(2), 1., mat%A, &
    mat%dimsA(1), mat%B, mat%dimsB(1), 0., C, mat%dimsA(1))

SGEMM 的意外结果

Unexpected results from SGEMM

fortran

gfortran

blas