Questions on "forrtl: error (65): floating invalid"

Question

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.

If any of the values of `C` are around `-tiny` then you still can end up with zero in the denominator. — SethMMorton, Aug 27 '14 at 21:27
Yep, that's true. But I guess the probability of that case is very small. Or, do you have ideas to avoid such situation? — elfsummer, Aug 27 '14 at 21:29
Did you print out the values of `C`? Do small values correlate to the infinities in `const`? — SethMMorton, Aug 27 '14 at 21:30
Well, there are "-Infinity" values in C(:,i), but C(:,i) is in denominator, which makes me confused... — elfsummer, Aug 27 '14 at 21:48
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. — Alexander Vogt, Aug 27 '14 at 23:08
update: after I change an error which is found on other supercomputer but not on this one (make the pointer to =>null() at declaration), then the "forrtl: error (65): floating invalid" error goes away...but the infinity still exists...= =# — elfsummer, Aug 28 '14 at 01:32
update2: but if I change the compiler option "-fpe3" to "-fpe0", then the error caused by the same line in source code (this time is floating overflow) still appears... — elfsummer, Aug 28 '14 at 03:14
Good idea! But I need to consider accuracy and it seems that there are infinity values passing into this subroutine which causes the error...so the problem might not be here but somewhere else. However, the -fpe0 might not be able to catch the Inf in the earliest place. I am still trying to figure out which line really causes the "floating overflow/invalid" problem. Thanks all the same =) — elfsummer, Aug 28 '14 at 22:41
Then, one last hint: Try to compile with `gfortran`. Different compilers create different compile-time/run-time warnings and/or errors and maybe `gfortran` detects your problem earlier. — Stefan, Aug 29 '14 at 05:30

score 1 · Answer 1 · edited Mar 17 '17 at 10:37

(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 in const?

@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 with gfortran. Different compilers create different compile-time/run-time warnings and/or errors and maybe gfortran detects your problem earlier.

Questions on "forrtl: error (65): floating invalid"

1 Answers1

Linked