1 FUNCTION pythag(a,b) 2 IMPLICIT NONE 3 INTEGER, PARAMETER :: sp = 8 4 REAL(sp), INTENT(IN) :: a,b 5 REAL(sp) :: pythag 6 !Computes (a2 + b2)1/2 without destructive under.ow or over.ow. 7 REAL(sp) :: absa,absb 8 absa=abs(a) 9 absb=abs(b) 10 if (absa > absb) then 11 pythag=absa*sqrt(1.0+(absb/absa)**2) 12 else 13 if (absb == 0.0) then 14 pythag=0.0 15 else 16 pythag=absb*sqrt(1.0+(absa/absb)**2) 17 end if 18 end if 19 END FUNCTION pythag