Return to drotg.f90 CVS log

Up to [local] / rpl / lapack / blas

Ajout de fichiers de lapack 3.11

    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