0

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'

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
Jame
  • 1
  • 1
  • 5
    The difference is only 10^-14, that is kind of expected for double precision, isn't it? The result depends on the precise ordering of all operations. What is your actual question? You did not ask one. – Vladimir F Героям слава Sep 14 '22 at 15:36
  • @Vladimir F Героям слава The question is, why the difference is not exactly zero? I ask the question because I have a code of thousands of lines, and I find simply replacing some direct computation of matrix multiplication by zgemm can give a different result, with difference 10^(-2), which is large to me. – Jame Sep 14 '22 at 16:09
  • See https://stackoverflow.com/questions/14994521/two-functions-different-result-suspected-fortran-feature or https://stackoverflow.com/questions/64790273/why-the-c-and-fortran-versions-of-this-same-program-produce-different-results – Vladimir F Героям слава Sep 14 '22 at 16:27
  • If some particular computation gives a large difference, we would need the specific details. – Vladimir F Героям слава Sep 14 '22 at 16:30
  • 3
    In your text the printed difference is 10^-14, not 10^-2. Floating point operations are approximate by design and for instance the result of `(a+b)+c` can be different from the result of `a+(b+c)` – PierU Sep 14 '22 at 20:04
  • 2
    @PierU nitpick - floating point maths is exact, but subtly different from the maths of the real numbers that we are used to. The approximation is that we are using floating point maths to represent real number maths, not in the maths of the floating point numbers themselves. – Ian Bush Sep 15 '22 at 07:29

1 Answers1

0

I copied your program as it is and compiled with the following command

gfortran file_name.f90 -lblas

and thereafter on execution (./a.out) found the difference to be zero as expected.

AITIJHYA SAHA
  • 1
  • 1
  • 1
    This is more a commant rather than an answer. Did you use Intel MKL as did the OP? It seems to me that you used Netlib BLAS instead. – Vladimir F Героям слава Jun 04 '23 at 11:01
  • See I think the issue is due to precision. In my machine when I write complex(8) it takes the real and imaginary part as double precision. For complex(4) it takes as single precision. However in the documentation, complex(16) is basically double precision. I think when you use it in the cluster using complex(8) with older version of gfortran it takes it as single precision and when you run it in your own machine it becomes double precision as in my case. I was having trouble with ZGEMM with complex(16), when I took the matrices as complex(8) the results were correct. – AITIJHYA SAHA Jun 05 '23 at 13:46
  • Do not confuse complex(16) and complex*16, they are different! Also, complex(8) is different from complex*8. If available at all, `complex(16)` is normally the quadruple precision, not double precision. See https://stackoverflow.com/questions/32935743/why-change-complex16-to-complex16-cause-the-runtime-increased-unreasonabl?noredirect=1&lq=1 https://stackoverflow.com/questions/838310/fortran-90-kind-parameter https://stackoverflow.com/questions/3170239/fortran-integer4-vs-integer4-vs-integerkind-4?noredirect=1&lq=1 – Vladimir F Героям слава Jun 05 '23 at 16:37
  • However, it is best to not use the magic numbers 4, 8 and 16, which are not portable, but use named constants with values such as `kind(1.0)` for single precision and `kind(1.d0)` for double precision. If you need the exact storage size, use `real32`, `real64` from `iso_fortran_env`. – Vladimir F Героям слава Jun 05 '23 at 16:44
  • Then complex(8) and complex*16 are same? – AITIJHYA SAHA Jun 06 '23 at 12:57
  • Only provided complex(8) actually exists. But yes, usually they are. – Vladimir F Героям слава Jun 06 '23 at 13:29