      DOUBLE PRECISION FUNCTION PYTHAG(A,B)
      DOUBLE PRECISION A,B
      INTEGER DDT,ERR,FMT,LCT(4),LIN(1024),LPT(6),HIO,RIO,WIO,RTE,WTE
      COMMON /IOP/ DDT,ERR,FMT,LCT,LIN,LPT,HIO,RIO,WIO,RTE,WTE
      DOUBLE PRECISION P,Q,R,S,T
      P = DMAX1(DABS(A),DABS(B))
      Q = DMIN1(DABS(A),DABS(B))
      IF (Q .EQ. 0.0D0) GO TO 20
      IF (DDT .EQ. 25) WRITE(WTE,1)
      IF (DDT .EQ. 25) WRITE(WTE,2) P,Q
    1 FORMAT(1X,'PYTHAG',1P2D23.15)
    2 FORMAT(1X,1P2D23.15)
   10 R = (Q/P)**2
      T = 4.0D0 + R
      IF (T .EQ. 4.0D0) GO TO 20
      S = R/T
      P = P + 2.0D0*P*S
      Q = Q*S
      IF (DDT .EQ. 25) WRITE(WTE,2) P,Q
      GO TO 10
   20 PYTHAG = P
      RETURN
      END

