Fortran - making subroutine more numerically stable

Author

making subroutine more numerically stable

Bart Vandewoestyne

2006-07-25, 8:00 am

I have the following subroutine which -for the ones that want to
know- represents a certain variable transformation to be used in
numerical integration:

subroutine f_mori(t, y, jacobian_det)
real(kind=qp), dimension(:), intent(in) :: t
real(kind=qp), dimension(:), intent(out) :: y
real(kind=qp), intent(out) :: jacobian_det

real(kind=qp) :: a, b

a = pi/2.0_qp
b = pi/4.0_qp

y = 0.5_qp*tanh( a*sinh(b*(1.0_qp/(1.0_qp-t) - 1.0_qp/t)) ) + 0.5_qp

jacobian_det = product( 0.5_qp*(1 - tanh(a*sinh(b * (1/(1-t) - 1/t)))**2) &
* a * cosh(b*(1/(1-t) - 1/t)) &
* b * (1/(1 - t)**2 + 1/t**2) )

end subroutine f_mori

Trying to evaluate this function for t=-1 to 2 in steps of 3.0/1000,
results in arithmetic exceptions with g95 and the F compiler. This happens
when the argument for the sinh function becomes too large.

Now my question is: what would be the best way to make this subroutine more
robust against such arithmetic exceptions? I assume something like

if the argument for the sinh function becomes larger than MAX_ARG then
set result to MAX
else
simply evaluate the sinh function, there is no danger for overflow
end if

Could this be the way to do it? And if `yes', then what do i take for MAX_ARG
and MAX? Something like MAX_ARG = log(huge(1.0_qp)/2) and MAX = huge(1.0_qp)/2
or so?

Any suggestions welcome!
Bart

--
"Share what you know. Learn what you don't."

highegg

2006-07-25, 7:00 pm

> y = 0.5_qp*tanh( a*sinh(b*(1.0_qp/(1.0_qp-t) - 1.0_qp/t)) ) + 0.5_qp
>
> jacobian_det = product( 0.5_qp*(1 - tanh(a*sinh(b * (1/(1-t) - 1/t)))**2) &
> * a * cosh(b*(1/(1-t) - 1/t)) &
> * b * (1/(1 - t)**2 + 1/t**2) )

The computation of y shouldn't overflow if the implementation of
TANH is smart enough to deal with infinities -
I'm sure that's the case of g95 and probably most compilers (where Inf
is available).
As for the derivatives, you-re getting NaNs because of multiplying
zeros
(1-tanh(...)**2) by Infs cosh(...)
you may, however, notice that
1-tanh(x)**2 = cosh(x)**(-2)
for any x and simplify the derivatives to the form
0.5_qp* a / cosh(b*(1/(1-t) - 1/t)) * b * (1/(1 - t)**2 + 1/t**2) )

which will probably behave better. But I guess you may
still end with PRODUCT multiplying an array of numbers with zeros and
Infs present,
thus the result will be NaN anyway.

Jaroslav

Terence

2006-07-27, 7:00 pm

Check your source.

> jacobian_det = product( 0.5_qp*(1 - tanh(a*sinh(b * (1/(1-t) - 1/t)))**2) &
> * a * cosh(b*(1/(1-t) - 1/t)) &
> * b * (1/(1 - t)**2 + 1/t**2) )

I don't like the fact you are mutltiplying the expression of line 1 by
the expressions of line 2 and line 3 to calculate thevariable "
jacobian_det".

One possibility is that one or both of those line-starting asterisks
should be a plus sign, because you have

( (term) *a* (term) * b * (term) .)

It would be unusual to separate the a and b multipliers.

Possibility number two is, are you sure you should not have:
jacobian_det = product( term1,term2,term3)

and product( ) is a real function?
In which case those leading asterisks should be commas.

Also you use terms 1/(1-t) and 1/t. So in ALL cases integer variable t
should NOT be either 0 or 1, which are two of your step values for this
stepping variable.

As to finding how to improve whatever you are doing, try forming
intermediate simpler terms and seeing where the exceptions occur.

Terence Wright

Sponsored Links