While trying to use a type for blockdiagonal matrices in Fortran code, I stumbled over a surprising bug in the following piece of code using the following compilers:
GNU Fortran (SUSE Linux) 7.4.0
ifort (IFORT) 18.0.5 20180823
ifort (IFORT) 16.0.1 20151021
If I compile
gfortran -Wall -Werror --debug ifort_bug.f && valgrind ./a.out
I have no errors reported from valgrind.
If I compile
ifort -warn all,error -debug -stacktrace ifort_bug.f && valgrind ./a.out
I get a segmentation fault for ifort_18 in my code and "only" memory leaks for ifort_16.
Is this a bug in the intel compilers, or does gfortran silently fix my bad code?
module blockdiagonal_matrices
implicit none
private
public :: t_blockdiagonal, new, delete,
& blocksizes, operator(.mult.), mult
save
integer, parameter :: dp = kind(1.d0)
type :: t_blockdiagonal
real(dp), allocatable :: block(:, :)
end type
interface new
module procedure block_new
end interface
interface delete
module procedure block_delete
end interface
interface operator (.mult.)
module procedure mult_blocks
end interface
contains
subroutine block_new(blocks, blocksizes)
type(t_blockdiagonal), intent(out) :: blocks(:)
integer, intent(in) :: blocksizes(:)
integer :: i, L
do i = 1, size(blocks)
L = blocksizes(i)
allocate(blocks(i)%block(L, L))
end do
end subroutine
subroutine block_delete(blocks)
type(t_blockdiagonal) :: blocks(:)
integer :: i
do i = 1, size(blocks)
deallocate(blocks(i)%block)
end do
end subroutine
function blocksizes(A) result(res)
type(t_blockdiagonal), intent(in) :: A(:)
integer :: res(size(A))
integer :: i
res = [(size(A(i)%block, 1), i = 1, size(A))]
end function
function mult_blocks(A, B) result(C)
type(t_blockdiagonal), intent(in) :: A(:), B(:)
type(t_blockdiagonal) :: C(size(A))
integer :: i
call new(C, blocksizes=blocksizes(A))
do i = 1, size(A)
C(i)%block = matmul(A(i)%block, B(i)%block)
end do
end function
subroutine mult(A, B, C)
type(t_blockdiagonal), intent(in) :: A(:), B(:)
type(t_blockdiagonal) :: C(:)
integer :: i
do i = 1, size(A)
C(i)%block = matmul(A(i)%block, B(i)%block)
end do
end subroutine
end module blockdiagonal_matrices
program time_blockdiagonal
use blockdiagonal_matrices
integer, parameter :: n_blocks = 2, L_block = 10**2
type(t_blockdiagonal) :: A(n_blocks), B(n_blocks), C(n_blocks)
integer :: i
integer :: start, finish, rate
call system_clock(count_rate=rate)
call new(A, blocksizes=[(L_block, i = 1, n_blocks)])
call new(B, blocksizes=[(L_block, i = 1, n_blocks)])
call new(C, blocksizes=[(L_block, i = 1, n_blocks)])
do i = 1, n_blocks
call random_number(A(i)%block)
call random_number(B(i)%block)
end do
call system_clock(start)
C = A .mult. B
call system_clock(finish)
write(6,*) 'Elapsed Time in seconds:',
& dble(finish - start) / dble(rate)
call system_clock(start)
call mult(A, B, C)
call system_clock(finish)
write(6,*) 'Elapsed Time in seconds:',
& dble(finish - start) / dble(rate)
call delete(A)
call delete(B)
call delete(C)
end program time_blockdiagonal
The ifort_18 segmentation fault is:
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
a.out 00000000004134BD Unknown Unknown Unknown
libpthread-2.26.s 00007FF720EB6300 Unknown Unknown Unknown
a.out 000000000040ABB8 Unknown Unknown Unknown
a.out 000000000040B029 Unknown Unknown Unknown
a.out 0000000000407113 Unknown Unknown Unknown
a.out 0000000000402B4E Unknown Unknown Unknown
libc-2.26.so 00007FF720B0AF8A __libc_start_main Unknown Unknown
a.out 0000000000402A6A Unknown Unknown Unknown
When I went into it with gdb it turned out, that the segfault is raised upon returning from function mult_blocks.
blockdiagonal_matrices::mult_blocks (c=..., a=..., b=...) at ifort_bug.f:63
63 do i = 1, size(A)
(gdb) s
64 C(i)%block = matmul(A(i)%block, B(i)%block)
(gdb) s
s
s
65 end do
(gdb) s
64 C(i)%block = matmul(A(i)%block, B(i)%block)
(gdb) s
65 end do
(gdb) s
66 end function
(gdb) s
Program received signal SIGSEGV, Segmentation fault.
0x000000000040abc4 in do_deallocate_all ()
(gdb) q
A debugging session is active.
Even with this information I cannot find a bug in my code.
I found a fix, although I don't understand why it works.
If I use -heap-arrays for compilation it works. So at first glance it seems to run into a Stackoverflow. If I do ulimit -s unlimited it does not solve the problem though.
If I explicitly allocate in the code it solves the issue.
subroutine new(blocks, blocksizes)
type(t_blockdiagonal), allocatable, intent(out) :: blocks(:)
integer, intent(in) :: blocksizes(:)
integer :: i, L
allocate(blocks(size(blocksizes)))
do i = 1, size(blocks)
L = blocksizes(i)
allocate(blocks(i)%block(L, L))
end do
end subroutine
subroutine delete(blocks)
type(t_blockdiagonal), allocatable :: blocks(:)
integer :: i
do i = 1, size(blocks)
deallocate(blocks(i)%block)
end do
deallocate(blocks)
end subroutine
function mult_blocks(A, B) result(C)
type(t_blockdiagonal), intent(in) :: A(:), B(:)
type(t_blockdiagonal), allocatable :: C(:)
integer :: i
call new(C, blocksizes=blocksizes(A))
do i = 1, size(A)
C(i)%block = matmul(A(i)%block, B(i)%block)
end do
end function
This way of writing is IMHO actually better than before and not a "dirty hack".
It is not anymore possible to call new with a differing size of blocksize and blocks.
Speaking C the type(t_blockdiagonal) :: blocks(n)
Should be a float** blocks[n] so just a vector of pointers to pointers.
The allocation of the actual blocks happened also in the first version on the heap. Hence I don't get the Stackoverflow for a vector that contains ca. 10 pointers.
According to this article, ulimit -s unlimited does not mean that the stack will be literally "unlimited":
The size of "unlimited" varies by Linux configuration, so you may need to specify a larger, specific number to ulimit (for example, 999999999). On Linux also note that many 32bit Linux distributions ship with a pthread static library (libpthread.a) that at runtime will fix the stacksize to 2093056 bytes regardless of the ulimit setting.
It is most likely that you do indeed have a stack overflow, given that -heap-arrays solves the issue.