I'm having a weird issue with a Fortran code

Question

I am fairly new to Fortran programming so this might be an obvious issue, so bear with me. Here is the code I am working with:

program A1H

! Householder transformation

implicit none
integer,parameter::dp=selected_real_kind(15,307) ! Double precision kind

real(kind=dp), dimension(6,3)::A
real(kind=dp), dimension(6,1)::b
integer, dimension(6,6)::Pglobal ! Global identity matrix
integer::i,j,g
g = size(A,1)
do j=1,g
   do i=1,g
       Pglobal(i,j) = (i/j)*(j/i)
   end do
end do
b(:,1) = [1237,1941,2417,711,1177,475]
A(1,:) = [1,0,0]
A(2,:) = [0,1,0]
A(3,:) = [0,0,1]
A(4,:) = [-1,1,0]
A(5,:) = [-1,0,1]
A(6,:) = [0,-1,1]
call mat_print('A',A)
call mat_print('b',b)
call mat_print('Pglobal',Pglobal)
call householder(A,b)

contains

subroutine householder(A,b)
real(kind=dp), intent(in)::A(:,:),b(:,:)
real(kind=dp)::alpha,gamma,beta
real(kind=dp), dimension(6,6)::H
real(kind=dp), dimension(6,3)::y,aa
real(kind=dp), dimension(6,1)::yy,v,dglobal,ek,bb
real(kind=dp), dimension(1,6)::d
integer::k,m,n,j
m = size(A,1)
n = size(A,2)
aa = A
bb = b
do k=1,n
   dglobal(:,k) = [0,0,0,0,0,0]
   alpha = -sign(aa(k,k),aa(k,k))*norm2(aa(k:m,k))
   ek(:,1) = Pglobal(:,k) 
   dglobal(k:m,k) = aa(k:m,k)
   v(:,k) = (dglobal(:,k)) - alpha*ek(:,1) 
   d(k,:) = v(:,k)
   beta = dot_product(d(k,:),v(:,k))
   if (beta==0) then
       continue
   end if
   H = Pglobal - (2/beta)*(matmul(reshape(v(:,k),(/m,1/)),reshape(d(k,:),(/1,m/)))) 
   y = matmul(H,aa)
   yy = matmul(H,bb)
   aa = y
   bb = yy
   call mat_print('aa',y)
   call mat_print('bb',yy)
end do

end subroutine

! Matrix print subroutine

subroutine mat_print(b,a)
    character(*), intent(in)::b
    class(*), intent(in)::a(:,:)
    integer::i
    print*,' '
    print*,b
    do i=1,size(a,1)
        select type (a)
        type is (real(kind=dp)) ; print'(100f9.4)',a(i,:)
        type is (integer) ; print'(100i9  )',a(i,:)
        end select
    end do
    print*,' '
 end subroutine

 end program

The issue I'm having is that when I change the variable aa to another name, I get the wrong result for the y variable; if I leave it as is, it is correct; however, leaving the bb variable as is, the yy result is incorrect, but if I change the bb variable to another other name, I get the correct result for yy, but then my answer for y changes! I'm very confused how this can happen as my experience from coding tells me that the answer should not change based on a simple variable name change, and if you look at the code, the y and yy variables are not even connected in the algorithm. This is not the only Fortran code I have run into this issue before on. Any help would be much appreciated.

Answer 1

I was able to reproduce your bug with GNU Fortran (Homebrew GCC 8.2.0) 8.2.0. There is indeed a bug in your program. You can find this bug by compiling with the option -fbounds_check. When you run it, you will find that several of your array access don't make sense. For example, you access dglobal(:,k) = [0,0,0,0,0,0], but the second dimension of dglobal is only 1. Use this flag to help fix your code, and I'm sure this bug will disappear.

For anyone who wants to dig deep into why this bug's appearance depends on the variable name, I was able to reproduce it with the array name test_array, but not with other shorter names. I was also able to get the correct answer using the test_array name if I set -fmax-stack-var-size=100, and other values appeared with different sizes. My guess is that gfortran puts these arrays on the stack, and the order is based on the name. Certain names put it in a "safe" location, so that the values are not overwritten by the buffer overflow.