1

I am having trouble with floating point exceptions in a subroutine that converts from geographic to geocentric coordinates. The variable geo(x) is fed into the subroutine as pairs of latitude and longitude. The variable xyz(x) is output as a triplet of components (x -- Greenwich meridian and equator; y -- 90 degree longitude and equator; z -- north pole).

subroutine geo2xyz(geo,xyz)
IMPLICIT REAL*8  (A-H,O-Z)
dimension geo(2),xyz(3)
rr=6367443.5
xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
xyz(2)=rr*sind(90.-geo(1))*sind(geo(2))
xyz(3)=rr*cosd(90.-geo(1))
return
end

I realize that sind and cosd are non-standard functions, but they are linked to an object file with the appropriate conversions. I have tested this before and it works for other codes that use degrees instead of radians: 

real function sind(x)
IMPLICIT REAL*8 (A-H,O-Z)
sind=sin(x*3.141592653589793d0/180.0d0)
return
end

real function cosd(x)
IMPLICIT REAL*8 (A-H,O-Z)
cosd=cos(x*3.141592653589793d0/180.0d0)
return
end

I tried going through the program with GDB, and couldn't figure out the problem. It seems that all the variables in the equation are OK, but xyz(1) = 1, which is not the case if I do that calculation with a calculator. 

Program received signal SIGFPE, Arithmetic exception.
0x0000000100001cb8 in geo2xyz (geo=..., xyz=...) at disslip.f:758
758     xyz(1)=rr*sind(90.-geo(1))*cosd(geo(2))
(gdb) print xyz(1)
$1 = 1
(gdb) print rr
$2 = 6367443.5
(gdb) print geo(1)
$3 = 50.350000000000001
(gdb) print geo(2)
$4 = -127.55800000000001
(gdb)

What am I missing here that is causing the floating point exception? I'm sure it's pretty simple, I'm very new to this.

Vladimir F Героям слава
  • 57,977
  • 4
  • 76
  • 119
tfinley
  • 21
  • 4
  • How are you calling the subroutine? What are the values of the arguments? Use `implicit none` and modules. It will save you from a lot of problems. Implicit none is the almost required thing for decent style for decades. – Vladimir F Героям слава Feb 09 '17 at 06:40
  • Are you sure the debugger gives the correct location? I would check it woth some print statements, this is strange. We must know which function is actually called. Add `external sind` and `external cosd` in the suvroutine or better use modules. – Vladimir F Героям слава Feb 09 '17 at 06:45
  • @VladimirF Unfortunately this is part of a large package of codes written by many people since the mid 80s, hence the obsolete style. It would take a monumental effort for me to update the style. The subroutine is called elsewhere like this: `call geo2xyz(rgeo,rxyz)` The values for rr, geo(1), and geo(2) are correct, but should not result in xyz(1)=1, so I assume the debugger has correctly identified a problem line. – tfinley Feb 09 '17 at 07:30
  • And whout about the `external`? I am suggesting those things to find out what is happening, not just for the style, they are a part of debugging. BTW if you delete all the `return`s the program remains the same, just shorter and easier to read. – Vladimir F Героям слава Feb 09 '17 at 07:33

1 Answers1

4

The sind and cosd functions do return a REAL (single precision) floating point.

But subroutine geo2xyz has no clue about that.
With the IMPLICIT REAL*8 it will assume that those functions return a double precision floating point.

Therefore, you must properly declare sind and cosd as real before using them, and/or change their return type to real*8.

The funny thing is that this program might have worked on old x86 architecture because the single precision value returned were promoted to double in the ST0 register. It won't work on 64 bit intel which is using SSE2 and register XMM0: there is no promotion to double in such architecture.

aka.nice
  • 9,100
  • 1
  • 28
  • 40
  • Bingo, I assumed it is some interface problem, modules prevent these errors. – Vladimir F Героям слава Feb 09 '17 at 08:46
  • Thanks @aka.nice, both of your solutions worked, as did simply compiling the object file with the `-freal-4-real-8` promotion flag. The FPE flags stopped signalling when I run the program containing `subroutine geo2xyz`. Unfortunately, this did not solve the greater issue at hand, which is that the main program is still calculating zero-values for several of the desired columns in the output file. But it seems that is a problem for another thread... – tfinley Feb 10 '17 at 03:06