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