为什么lapack dtrmm.f似乎无法正常工作

Question

我决定使用lapack子程序dtrmm而不是matmul来乘以下三角(d，d(矩阵和一般(d，n(矩阵。然而，它似乎不能正常工作。以下代码比较matlum(顶部(和dtrmm的结果

Program test
implicit none
integer, parameter :: n = 3, d = 3 ! arbitrary numbers, no luck with other values
integer :: i
real(8) :: a(d,d), b(d,n), c(d,n)
call random_number(a)
call random_number(b)
do i=2,d
a(1:i-1,i) = 0
end do
c = matmul(a,b)
call dtrmm('L','L','N','N',d,n,1,a,d,b,d) ! documentation linked in the question
print*, 'correct result : '
do i=1,d
print*, c(i,:)
end do
print*, 'dtrmm yields : '
do i=1,d
print*, b(i,:)
end do
End Program test

返回此

correct result : 
0.75678922130735249       0.51830058846965921       0.51177304237548271     
1.1974740765385026       0.46115246753697681       0.98237114765741340     
0.98027798241945430       0.53718796235743815        1.0328498570683342     
dtrmm yields : 
6.7262070844500211E+252   4.6065628207770121E+252   4.5485471599475983E+252
1.0642935166471391E+253   4.0986405551607272E+252   8.7311388520015068E+252
8.7125351741793473E+252   4.7744304178222945E+252   9.1797845822711462E+252

我使用的其他lapack子例程运行良好。是什么导致了这种行为不端？这是一个bug，还是我只是严重误解了什么？

Answer 1

这是一个简单的数据类型错误。因子alpha的类型必须是双精度，而您提供的是默认类型的整数。

因此

...
! call dtrmm('L','L','N','N',d,n,1,a,d,b,d)  WRONG
call dtrmm('L','L','N','N',d,n,1d0,a,d,b,d)  ! note 1d0 instead of 1
...

给出了正确的结果。