Return to drotg.f90 CVS log

Up to [local] / rpl / lapack / blas

Annotation of rpl/lapack/blas/drotg.f90, revision 1.1

1.1     ! bertrand    1: !> \brief \b DROTG
        !             2: !
        !             3: !  =========== DOCUMENTATION ===========
        !             4: !
        !             5: ! Online html documentation available at
        !             6: !            http://www.netlib.org/lapack/explore-html/
        !             7: !
        !             8: !  Definition:
        !             9: !  ===========
        !            10: !
        !            11: !  DROTG constructs a plane rotation
        !            12: !     [  c  s ] [ a ] = [ r ]
        !            13: !     [ -s  c ] [ b ]   [ 0 ]
        !            14: !  satisfying c**2 + s**2 = 1.
        !            15: !
        !            16: !> \par Purpose:
        !            17: !  =============
        !            18: !>
        !            19: !> \verbatim
        !            20: !>
        !            21: !> The computation uses the formulas
        !            22: !>    sigma = sgn(a)    if |a| >  |b|
        !            23: !>          = sgn(b)    if |b| >= |a|
        !            24: !>    r = sigma*sqrt( a**2 + b**2 )
        !            25: !>    c = 1; s = 0      if r = 0
        !            26: !>    c = a/r; s = b/r  if r != 0
        !            27: !> The subroutine also computes
        !            28: !>    z = s    if |a| > |b|,
        !            29: !>      = 1/c  if |b| >= |a| and c != 0
        !            30: !>      = 1    if c = 0
        !            31: !> This allows c and s to be reconstructed from z as follows:
        !            32: !>    If z = 1, set c = 0, s = 1.
        !            33: !>    If |z| < 1, set c = sqrt(1 - z**2) and s = z.
        !            34: !>    If |z| > 1, set c = 1/z and s = sqrt( 1 - c**2).
        !            35: !>
        !            36: !> \endverbatim
        !            37: !
        !            38: !  Arguments:
        !            39: !  ==========
        !            40: !
        !            41: !> \param[in,out] A
        !            42: !> \verbatim
        !            43: !>          A is DOUBLE PRECISION
        !            44: !>          On entry, the scalar a.
        !            45: !>          On exit, the scalar r.
        !            46: !> \endverbatim
        !            47: !>
        !            48: !> \param[in,out] B
        !            49: !> \verbatim
        !            50: !>          B is DOUBLE PRECISION
        !            51: !>          On entry, the scalar b.
        !            52: !>          On exit, the scalar z.
        !            53: !> \endverbatim
        !            54: !>
        !            55: !> \param[out] C
        !            56: !> \verbatim
        !            57: !>          C is DOUBLE PRECISION
        !            58: !>          The scalar c.
        !            59: !> \endverbatim
        !            60: !>
        !            61: !> \param[out] S
        !            62: !> \verbatim
        !            63: !>          S is DOUBLE PRECISION
        !            64: !>          The scalar s.
        !            65: !> \endverbatim
        !            66: !
        !            67: !  Authors:
        !            68: !  ========
        !            69: !
        !            70: !> \author Edward Anderson, Lockheed Martin
        !            71: !
        !            72: !> \par Contributors:
        !            73: !  ==================
        !            74: !>
        !            75: !> Weslley Pereira, University of Colorado Denver, USA
        !            76: !
        !            77: !> \ingroup single_blas_level1
        !            78: !
        !            79: !> \par Further Details:
        !            80: !  =====================
        !            81: !>
        !            82: !> \verbatim
        !            83: !>
        !            84: !>  Anderson E. (2017)
        !            85: !>  Algorithm 978: Safe Scaling in the Level 1 BLAS
        !            86: !>  ACM Trans Math Softw 44:1--28
        !            87: !>  https://doi.org/10.1145/3061665
        !            88: !>
        !            89: !> \endverbatim
        !            90: !
        !            91: !  =====================================================================
        !            92: subroutine DROTG( a, b, c, s )
        !            93:    integer, parameter :: wp = kind(1.d0)
        !            94: !
        !            95: !  -- Reference BLAS level1 routine --
        !            96: !  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
        !            97: !  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
        !            98: !
        !            99: !  .. Constants ..
        !           100:    real(wp), parameter :: zero = 0.0_wp
        !           101:    real(wp), parameter :: one  = 1.0_wp
        !           102: !  ..
        !           103: !  .. Scaling constants ..
        !           104:    real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( &
        !           105:       minexponent(0._wp)-1, &
        !           106:       1-maxexponent(0._wp) &
        !           107:    )
        !           108:    real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( &
        !           109:       1-minexponent(0._wp), &
        !           110:       maxexponent(0._wp)-1 &
        !           111:    )
        !           112: !  ..
        !           113: !  .. Scalar Arguments ..
        !           114:    real(wp) :: a, b, c, s
        !           115: !  ..
        !           116: !  .. Local Scalars ..
        !           117:    real(wp) :: anorm, bnorm, scl, sigma, r, z
        !           118: !  ..
        !           119:    anorm = abs(a)
        !           120:    bnorm = abs(b)
        !           121:    if( bnorm == zero ) then
        !           122:       c = one
        !           123:       s = zero
        !           124:       b = zero
        !           125:    else if( anorm == zero ) then
        !           126:       c = zero
        !           127:       s = one
        !           128:       a = b
        !           129:       b = one
        !           130:    else
        !           131:       scl = min( safmax, max( safmin, anorm, bnorm ) )
        !           132:       if( anorm > bnorm ) then
        !           133:          sigma = sign(one,a)
        !           134:       else
        !           135:          sigma = sign(one,b)
        !           136:       end if
        !           137:       r = sigma*( scl*sqrt((a/scl)**2 + (b/scl)**2) )
        !           138:       c = a/r
        !           139:       s = b/r
        !           140:       if( anorm > bnorm ) then
        !           141:          z = s
        !           142:       else if( c /= zero ) then
        !           143:          z = one/c
        !           144:       else
        !           145:          z = one
        !           146:       end if
        !           147:       a = r
        !           148:       b = z
        !           149:    end if
        !           150:    return
        !           151: end subroutine