# Handling underflow errors in fortran 90

Tags:
1. Jul 10, 2017

### meteo student

I am using fortran 90 to solve a elliptic partial differential equation numerically(Poisson's equation) and I have encountered a underflow error. All the arrays involved are stored as real data types. The original expression is highlighted in the comments and the my attempt at refactoring into simpler expressions is shown as expr0, expr1, expr2,expr3. I am not sure about my attempt to handle the underflow error because while the FPE problem goes away the solution does not converge and I get all zero values.

Any suggestions to improve my coding attempts ?
Code (Fortran):

use,intrinsic :: ieee_arithmetic
use, intrinsic :: iso_fortran_env, only: compiler_version, compiler_options
implicit none
logical :: underflow_support, underflow

do k=1,nz-1
do i=1,nx-1
do j=1,ny-1
expr0 = omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*&
(psi(i+1,j,k)+psi(i-1,j,k)))
expr1 = b(i,j,k)*psi(i,j+1,k)+b(i,j-1,k)&
*psi(i,j-1,k)+c(i,j,k)*&
psi(i,j,k+1)+c(i,j,k-1)&
*psi(i,j,k-1)
expr1 = expr0 + expr1
expr2 = (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k-1)+c(i,j,k))

call ieee_get_flag(ieee_underflow,underflow)
if (underflow) then
end if
if (expr1 == 0.0) then
expr3 = expr1/expr2
expr3 = expr3 + (1.-omega)*psi(i,j,k)
psi(i,j,k) = expr3
end if

!
!             psi(i,j,k)=omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*
!   >              (psi(i+1,j,k)+psi(i-1,j,k))
!   >              +b(i,j,k)*psi(i,j+1,k)+b(i,j-1,k)
!   >              *psi(i,j-1,k)+c(i,j,k)*
!   >              psi(i,j,k+1)+c(i,j,k-1)
!   >              *psi(i,j,k-1))/
!   >              (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)
!   >              +c(i,j,k-1)+
!   >              c(i,j,k))
!                 +(1.-omega)*psi(i,j,k)

enddo
enddo
enddo

Last edited: Jul 11, 2017
2. Jul 11, 2017

### Staff: Mentor

I don't see how this could possibly work.
1. There is no PROGRAM section. Possibly you just copied a short section of your code.
2. Your three DO loop run from 1 to, respectively, nz-1, nx-1, and ny-1. None of the nx, ny, nz variables has been initialized, at least in the section of code that you posted.
3. I don't see any declarations for any variables, including omega and dxyz, or for any of the arrays.
When I see obvious errors like these, I stop looking.

Since you asked how you might improve your coding attempts, here are some suggestions.
• Add spaces in your code. The compiler doesn't care, but jamming everything together makes it harder for humans to read and understand what you're trying to do. Here is one of your assignment statements, written on one line:
Code (Fortran):
expr0 = omega*(-dxyz*rhoref(2*k)*pv(i,j,k)+a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)))
Granted, you wrote it on two lines, using &, but it's almost unreadable with everything jammed together.
• Use more informative names -- rhoref means what? And what is dxyz?
Single-letter names are OK for loop control variables (i, j, and k in your code), but you really should have better names than nx, ny, nz.
• Where's the rest of the code? The only declarations I see are for underflow_support and underflow, neither of which has been initialized.
• If your code isn't working correctly, use a debugger, and insert a breakpoint to inspect the values of the variables being manipulated. If you don't know how to use a debugger, invest some time in doing this. Any time you spend will be amply rewarded. In the meantime, use the tried and true method of inserting write statements to print the values of variables. When your code is working as expected, you can delete or comment out the extra write statements.

3. Jul 11, 2017

### meteo student

Thank you Mark for the response. I apologize for not adding the code tag. I added the tag now.
Thanks also for the suggestions on the readability.
Yes I already used a debugger (gdb) and could not make much progress with it.I am using gfortran and I added the breakpoint on the offending line. As the values for the arrays are very large gdb outputs these statements -
<error reading variable: value requires 250000 bytes, which is more than max-value-size>.
So I set max-value unlimited. After that I tried to run gdb it hangs without doing anything.
Then I removed following lines during compilation fpe=overflow,underflow,zero,denormal and then ran it and got an IEEE underflow error. I tried printing the values and there were some very small values - 1.E-45 etc.

I could post the full code but I wanted to avoid it as I had a specific question i.e. how to handle a underflow. Is there a way to flush the really low values to zero and then move on ?

Last edited: Jul 11, 2017
4. Jul 11, 2017

### Tom.G

Code (Fortran):

if (expr1 == 0.0) then
expr3 = expr1/expr2

Here you are trying to divide Zero by something. Even with pencil and paper this does not work well.

5. Jul 11, 2017

### meteo student

Tom - I agree. If you can read my response to Mark I have summarized my efforts so far. My question was how to handle an underflow error and move on.

So if expr1 is a really low value such as -1.0E-45 is it reasonable to flush it to zero and move on ? That is the question I am asking here.

6. Jul 11, 2017

### Staff: Mentor

As a minor nit, I'd reformat this piece of code to make it even more readable:

Code (Fortran):

expr0 = omega * ( (-dxyz) * rhoref(2*k) * pv(i, j, k) + a(i, j, k) * ( psi(i+1, j, k) + psi(i-1, j, k) ) )

expr1 =    b(i, j, k)    * psi(i, j+1, k)    &
+ b(i ,j-1, k)  * psi(i, j-1, k)    &
+ c(i, j, k)    * psi(i, j, k+1)    &
+ c(i, j, k-1)  * psi(i, j, k-1)

or

Code (Fortran):

expr0 = omega * ( (-dxyz) * rhoref(2*k) * pv(i, j, k)
+ a(i, j, k) * ( psi(i+1, j, k) + psi(i-1, j, k) ) )

expr1 = b(i, j, k)    * psi(i, j+1, k)    &
+ b(i ,j-1, k)  * psi(i, j-1, k)    &
+ c(i, j, k)    * psi(i, j, k+1)    &
+ c(i, j, k-1)  * psi(i, j, k-1)

If you spend some time now reformatting it then you'll simplify your debugging effort later on especially if you have to revisit it 6 months later to fix something.

7. Jul 11, 2017

### Staff: Mentor

Sure it does. There's no problem dividing 0.0 by a nonzero value -- you get 0.0 as the answer. The problem comes when you divide by zero.

8. Jul 11, 2017

### Staff: Mentor

250KB isn't all that big. Can you show the line of code that is causing this problem? Be sure to include declarations for all variables, and their values before that line is hit.

gdb shouldn't hang. Once you hit a breakpoint, the program will stop, and you can single-step through succeeding lines of code.
Seems reasonable to me.