I've been running a huge Fortran code with Intel compiler version 13.1.3.192 under the debug mode (with -O0 -g -traceback -fpe3
flags being turned on). It gave me the following output message:
... ...
forrtl: warning (402): fort: (1): In call to MPI_ALLGATHER, an array temporary was created for argument #1
forrtl: error (65): floating invalid
Image PC Routine Line Source
arts 00000000016521D9 pentadiagonal_ser 385 pentadiagonal.f90
arts 0000000001644862 pentadiagonal_ 62 pentadiagonal.f90
arts 00000000004DF167 implicit_solve_x_ 1201 implicit.f90
arts 0000000000538079 scalar_bquick_inv 383 scalar_bquick.f90
arts 00000000004EFEAC scalar_step_ 190 scalar.f90
arts 0000000000401744 simulation_run_ 182 simulation.f90
arts 0000000000401271 MAIN__ 10 main.f90
arts 0000000000400FAC Unknown Unknown Unknown
arts 000000000420E444 Unknown Unknown Unknown
arts 0000000000400E81 Unknown Unknown Unknown
and the source of the error comes from the subroutine pentadiagonal_serial, which is to solve a penta-diagonal matrix:
subroutine pentadiagonal_serial(A,B,C,D,E,R,n,lot)
use precision
implicit none
integer, intent(in) :: n,lot
real(WP), dimension(lot,n) :: A ! LOWER-2
real(WP), dimension(lot,n) :: B ! LOWER-1
real(WP), dimension(lot,n) :: C ! DIAGONAL
real(WP), dimension(lot,n) :: D ! UPPER+1
real(WP), dimension(lot,n) :: E ! UPPER+2
real(WP), dimension(lot,n) :: R ! RHS - RESULT
real(WP), dimension(lot) :: const
integer :: i
if (n .eq. 1) then
! Solve 1x1 system
R(:,1) = R(:,1)/C(:,1)
return
else if (n .eq. 2) then
! Solve 2x2 system
const(:) = B(:,2)/C(:,1)
C(:,2) = C(:,2) - D(:,1)*const(:)
R(:,2) = R(:,2) - R(:,1)*const(:)
R(:,2) = R(:,2)/C(:,2)
R(:,1) = (R(:,1) - D(:,1)*R(:,2))/C(:,1)
return
end if
! Forward elimination
do i=1,n-2
! Eliminate A(2,i+1)
const(:) = B(:,i+1)/(C(:,i)+tiny(1.0_WP))
C(:,i+1) = C(:,i+1) - D(:,i)*const(:)
D(:,i+1) = D(:,i+1) - E(:,i)*const(:)
R(:,i+1) = R(:,i+1) - R(:,i)*const(:)
in which the line
const(:) = B(:,i+1)/(C(:,i)+tiny(1.0_WP))
causes the error. I tried to print out the values of const(:)
and find that there did exist Infinity
values. However, I cannot understand why it can generate infinities. As far as I see, to avoid zero denominator, the tiny(1.0_WP)
is added to C(:,i)
and now it is nearly impossible for the denominator to be zero... I also checked when this subroutine is being called, everything is initialized or given a value after declaration. So I couldn't figure out where could go wrong.
(Answers in the comments. See Question with no answers, but issue solved in the comments (or extended in chat). There is a lot of chat in the comments and its hard to extract the actual answer, but the OP indicates it is solved. )
@SethMMorton wrote:
If any of the values of
C
are around-tiny
then you still can end up with zero in the denominator.Did you print out the values of
C
? Do small values correlate to the infinities inconst
?
@Alexander Vogt wrote:
Well, adding something to infinity results in a floating point exception... And I'm not sure ifort can detect that a division by infinity is zero, I guess that is another floating point exception! You'd better make sure that no infinite values occur inside C.
@Stefan wrote:
Workaround: What about the routine
?gbsv
from MKL/Lapack? Then, one last hint: Try to compile withgfortran
. Different compilers create different compile-time/run-time warnings and/or errors and maybegfortran
detects your problem earlier.