Inconsistent results from zgemm and direct computation

Viewed 53

In the school's cluster, I find using a lapack subroutine zgemm can get different results from direct computation when performing matrix multiplication. Here is a sample code:

program main
    implicit none 

    integer, parameter :: Dcut=64

    complex(8), dimension(Dcut,Dcut) :: Lmat
    complex(8), dimension(Dcut,Dcut**2) :: Tmat, outmat1, outmat2

    integer :: a,b,c
    real(8) :: Lr(Dcut,Dcut), Li(Dcut,Dcut), Tr(Dcut,Dcut**2), Ti(Dcut,Dcut**2)


    ! create Lmat and Tmat
    call random_number( Lr ) 
    call random_number( Li )
    call random_number( Tr )
    call random_number( Ti )

    Lmat = cmplx( Lr, Li, kind=8 )
    Tmat = cmplx( Tr, Ti, kind=8 )

    ! use zgemm to compute Lmat.Tmat
    call zgemm( 'N','N',Dcut,Dcut**2,Dcut,(1.0D0,0.0D0),Lmat,Dcut,Tmat,Dcut,(0.0D0,0.0D0),&
                 outmat1,Dcut )

    ! direct computation of Lmat.Tmat
    outmat2 = ( 0.0D0,0.0D0 )
    do b=1, Dcut**2
        do a=1, Dcut
            do c=1, Dcut
                outmat2(a,b) = outmat2(a,b) + Lmat( a,c ) * Tmat( c,b )
            end do 
        end do 
    end do 


    print *,'diff', maxval( abs( outmat1-outmat2 ) )


end program main

And the output is: ' diff 4.32958511410900557E-014 '

Why two methods give different results, although the difference is very small?

I use gfortran and the version on the cluster is:

' gcc version 4.4.7 20120313 (Red Hat 4.4.7-11) (GCC) '

I use intel mkl and the version seems to be ' 2018.5.274 '

The code is run on CentOS release 6.6 (Final).

But when running on my own PC (wsl Ubuntu), the output is exactly 0.

For my PC, the version of gfortran is 'gcc version 9.4.0'

0 Answers
Related