I decided to use lapack subroutine dtrmm
instead of matmul
to multiply lower triangular (d,d) matrix and general (d,n) matrix. However, it doesn't seem to work correctly. The following code compares results of matlum
(top) and 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 this
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
Other lapack suboutines I used work fine. What could be causing this to misbehave? Is this a bug, or have I just badly misunderstood something?
It is a simple data type error. The factor alpha
must be of type double precision whereas you supplied an integer of default kind.
Thus
...
! 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
...
gives the correct result.