mat_pythag Function

public function mat_pythag(a, b)

Arguments

Type IntentOptional Attributes Name
doubleprecision :: a
doubleprecision :: b

Return Value doubleprecision


Source Code

doubleprecision function mat_pythag(a,b)
doubleprecision :: a
doubleprecision :: b

doubleprecision :: p
doubleprecision :: q
doubleprecision :: r
doubleprecision :: s
doubleprecision :: t

   p = dmax1(dabs(a),dabs(b))
   q = dmin1(dabs(a),dabs(b))

   if (q .ne. 0.0d0) then

      infinite : do
         r = (q/p)**2
         t = 4.0d0 + r
         if (t .eq. 4.0d0) exit infinite
         s = r/t
         p = p + 2.0d0*p*s
         q = q*s
      enddo infinite

   endif

   mat_pythag = p
end function mat_pythag