From 1ce7be9da747a06a27ed8c880c4be740d27c5331 Mon Sep 17 00:00:00 2001 From: allen Date: Mon, 24 Jan 2000 08:31:55 +0000 Subject: Fix from Thomas Radke for floating point exception on the Alphas caused by exponention by a real number .... here's Thomas's description: The include file generated by Mathematica contains an exponentiation operation with a REAL exponent. My Fortran book says that y**x = exp(x * ln(y)) On O2K this intrinsic function checks the exponent against 0, on Alphas it doesn't and you get a FP exception if y <= 0. I just changed the type of the exponent to be INTEGER instead of REAL which was probably also the intention of the programmer. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/IDAxiBrillBH/trunk@8 0a4070d5-58f5-498f-b6c0-2693e757fa0f --- src/IDAxiBrillBH.F | 1 + src/bhbrill.x | 10 ++++++++-- 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/src/IDAxiBrillBH.F b/src/IDAxiBrillBH.F index d064efc..01000ef 100644 --- a/src/IDAxiBrillBH.F +++ b/src/IDAxiBrillBH.F @@ -62,6 +62,7 @@ c Perhaps this and others should go into cctk.h real*8 o70,o71,o72,o73,o74,o75,o76,o77,o78,o79 real*8 o80,o81,o82,o83,o84,o85,o86,o87,o88,o89 real*8 o90,o91,o92,o93,o94,o95,o96,o97,o98,o99 + integer i22 real*8 pi real*8 adm integer :: nx,ny,nz diff --git a/src/bhbrill.x b/src/bhbrill.x index b7721c9..270519b 100644 --- a/src/bhbrill.x +++ b/src/bhbrill.x @@ -19,8 +19,14 @@ o19 = -(o15*o18) o20 = o16 + o17 + o19 o21 = exp(o20) - o22 = -2.00000000000000d0 + n - o23 = o8**o22 + i22 = -2 + n +c The next lines causes a floating point exception on Alphas, +c if o22 is 0.0 but o8 is negative. +c This is because of y**x = exp(x*ln(y)). +c Changing the type of the exponent to integer fixes the problem. +c o22 = -2.00000000000000d0 + n +c o23 = o8**o22 + o23 = o8**i22 o24 = n**2 o25 = 2.00000000000000d0*etagrd(i)*eta0*o15 o26 = o16 + o19 + o25 -- cgit v1.2.3