      DOUBLE PRECISION FUNCTION PYTHAG (A, B)
      IMPLICIT NONE
C
      DOUBLE PRECISION A, B
C
      INCLUDE MATLAB$KOM:IOP.KOM
C
      DOUBLE PRECISION P, Q, R, S, T
C
C
      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)
1     FORMAT (' PYTHAG', 1P2D23.15)
      IF (DDT.EQ.25) WRITE (WTE, 2) P, Q
2     FORMAT (1X, 1P2D23.15)
C
10    CONTINUE
      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
C
20    CONTINUE
      PYTHAG = P
C
      RETURN
      END
