Search code examples
fortranlapack

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


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?


Solution

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