为什么 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

returns 这个

 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 子程序工作正常。什么可能导致这种行为不端？这是一个错误，还是我误解了什么？

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
...

给出了正确的结果。

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

Why does lapack dtrmm.f seem to not work properly?

fortran

lapack