In the code below I am having the following issue which arose while trying to write a wrapper (which I’m still struggling with and could use some help) to zgemm in BLAS so that I can use it more easily:
program zgemm_wrap
implicit none
integer, parameter :: N = 2
double complex, dimension(N, N) :: A, B
real, dimension(N, N) :: x
real, dimension(2) :: y
y = (/ 1.0, 2.0 /)
A = reshape((/ 1.0, 2.0, 3.0, 4.0 /), (/ 2, 2/))
call ZGEMM(‘N’, ‘N’, N, N, N, 1.0d0, A, N, A, N, 1.0d0, B, N)
write(*, *) B
write(*, *) matmul(A, A)
end program zgemm_wrap
If I comment out line 8 (i.e. the initialization of y) the results from zgemm are correct, namely (7, 10, 15, 22).
However, if I uncomment it, the results are incorrect; (7 + i14, 10 + i20, 15 + i30, 22 + i44). That is, the imaginary component becomes non-zero (which it shouldn’t) and curiously has a value twice that of the real component. I’ve had a similar problem in the past where the presence of a write statement fixed or destroyed the calculations
I’m using PGI 19.10 community edition on windows 10
I’m not sure what elso to include regarding my machine.
final remarks:
(1) Compiling with gcc (and pgi) on a linux machine at my university does not reproduce the error
(2) Running with the option -C for the pgi compiler seems to fix the error, so I suspect it has to do with array bounds.
(3) Due to the above I initially thought it may be an issue with BLAS on my machine, since I am just using what came with the pgi program.
(4) I just realized that writing out the variable y at the end of the program solved the issue, but I have no idea why. Curiously, it doesn’t fix the problem if I wrote it out at the beginning. I suspect the issue lies with y
(5) Just to be sure, I compile as: pgf90 zgemm_wrap.f90 -lblas, and then run the following executable.
thanks for any help!