1

I have a very weird problem when I compile Fortran code using the zgerc subroutine from BLAS. Basically, this subroutine calculate the outer product of vector x with the conjugate of vector y. More about that function here. My simple code is as following:

program main
    implicit none
    integer :: i
    complex(8), dimension(10) :: a = [(i, i=0,9)]
    complex(8), dimension(10) :: b = [(i, i=0,9)]
    complex(8), dimension(10, 10) :: c
    c = 0
    CALL zgerc(10, 10, 1.D0, a, 1, b, 1, c, 10)
    WRITE(*, *) c
end program main

I have here 2 complex vectors, a and b, both goes from 0 to 9 and their imaginary part is 0.

Now for the weird part. If I compile that code with absolute path: gfortran -c /home/myUser/Fortran/tests/main.f90 -o main.o I get the correct result, but if I compile with gfortran -c main.f90 -o main.o (of course I'm in the right directory, and I've also tried ./main.f90) the result of the real part is correct, but for the imaginary part I get numbers like 1E+225 (and if I use ./main.f90 I get numbers like 1E+163).

I can't understand why the path to my code change the result of the imaginary part... I'll be glad for your help.

I use Ubuntu 20.04.2 with the default gfortran (9.3.0)

P.S, my final goal is to use this as part of a more complex subroutine in Python with f2py.

EDIT: my full commands:

#gfortran -c /home/myUser/Fortran/tests/main.f90 -o main.o
gfortran -c main.f90 -o main.o
gfortran -o test main.o /home/myUser/PycharmProjects/GSIE_2D/fortran_scripts/libblas.a /home/myUser/PycharmProjects/GSIE_2D/fortran_scripts/liblapack.a
rm ./main.o
./test

line 1 and 2 are the 2 cases, so I run only one of them each time.

Ian Bush
  • 6,996
  • 1
  • 21
  • 27
dorzv
  • 31
  • 4
  • Welcome, please take the [tour]. The first thing is to double-check that there are no old versions of the executable around. First go through your directories and delete all `a.out`. Make absolutely sure that there is only one resulting binary around and that you absolutely know which you are executing. Your description does not contain a complete sequence of your commands. If the check does not resolve the problem, please supply the *full and exact sequence of your commands* including your running command. – Vladimir F Героям слава Apr 15 '21 at 08:00
  • Thanks for your reply. I remove the old a.o files each time. I edited the question with my full commands. – dorzv Apr 15 '21 at 08:12
  • 1
    The constant alpha (the third argument) should be "double complex", not "double precision". If you fix this does this solve your problem? Also note complex(8) is not portable and is poor practice - see https://stackoverflow.com/questions/838310/fortran-90-kind-parameter and note that the kind of a complex is the same as the kind of the reals that compose it. I use the iso_fortran_env route personally nowadays. – Ian Bush Apr 15 '21 at 08:39
  • `1.d0` is not a complex, which the documentation says you should pass, so you can't expect any particular result at all. – francescalus Apr 15 '21 at 08:39

1 Answers1

1

You supply 1d0 which is a double precision literal whereas zgerc assumes a double complex value.

call zgerc(10, 10, cmplx(1, kind=8), a, 1, b, 1, c, 10)

By including explicit interfaces (via some kind of blas module) you would have gotten compile time errors when supplying arguments of wrong data type. Intel's mkl provides such explicit interfaces in its blas95 module as well as generic routines (gerc instead of {c,z}gerc). There is also this open-source module providing explicit interfaces to standard blas routines.

Also make use of the portable types defined in iso_fortran_env.

program main
  use blas95,          only: gerc
  use iso_fortran_env, only: real64
  implicit none

  integer, parameter :: n = 10
  integer            :: i
  complex(real64)    :: a(n) = [(i, i=0,n-1)], b(n) = [(i, i=0,n-1)], c(n,n)

  c = 0
  ! call zgerc(10, 10, cmplx(1, kind=8), a, 1, b, 1, c, 10) ! standard blas 
  call gerc(c, a, b, alpha=cmplx(1, kind=real64))           !  generic blas95 
  print *, c
end program
jack
  • 1,658
  • 1
  • 7
  • 18
  • Is there a specification for the BLAS95 interface anywhere? Is it only available for MKL or is it more widespread than that? Google is not showing up much to clarify – Ian Bush Apr 15 '21 at 10:21
  • Thanks, using parameter `(1.D0, 0.D0)` actually solved the problem. I (mistakenly) believed the compiler will know how to handle this. I'll try the BLAS95 also. – dorzv Apr 15 '21 at 10:22
  • 1
    @IanBush I always thought BLAS95 would be some open standard but it might just be intel specific. I edited the answer accordingly. – jack Apr 15 '21 at 12:37