Annotation of rpl/lapack/blas/zrotg.f90, revision 1.1
1.1 ! bertrand 1: !> \brief \b ZROTG generates a Givens rotation with real cosine and complex sine.
! 2: !
! 3: ! =========== DOCUMENTATION ===========
! 4: !
! 5: ! Online html documentation available at
! 6: ! http://www.netlib.org/lapack/explore-html/
! 7: !
! 8: ! Definition:
! 9: ! ===========
! 10: !
! 11: ! ZROTG constructs a plane rotation
! 12: ! [ c s ] [ a ] = [ r ]
! 13: ! [ -conjg(s) c ] [ b ] [ 0 ]
! 14: ! where c is real, s is complex, and c**2 + conjg(s)*s = 1.
! 15: !
! 16: !> \par Purpose:
! 17: ! =============
! 18: !>
! 19: !> \verbatim
! 20: !>
! 21: !> The computation uses the formulas
! 22: !> |x| = sqrt( Re(x)**2 + Im(x)**2 )
! 23: !> sgn(x) = x / |x| if x /= 0
! 24: !> = 1 if x = 0
! 25: !> c = |a| / sqrt(|a|**2 + |b|**2)
! 26: !> s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2)
! 27: !> r = sgn(a)*sqrt(|a|**2 + |b|**2)
! 28: !> When a and b are real and r /= 0, the formulas simplify to
! 29: !> c = a / r
! 30: !> s = b / r
! 31: !> the same as in DROTG when |a| > |b|. When |b| >= |a|, the
! 32: !> sign of c and s will be different from those computed by DROTG
! 33: !> if the signs of a and b are not the same.
! 34: !>
! 35: !> \endverbatim
! 36: !
! 37: ! Arguments:
! 38: ! ==========
! 39: !
! 40: !> \param[in,out] A
! 41: !> \verbatim
! 42: !> A is DOUBLE COMPLEX
! 43: !> On entry, the scalar a.
! 44: !> On exit, the scalar r.
! 45: !> \endverbatim
! 46: !>
! 47: !> \param[in] B
! 48: !> \verbatim
! 49: !> B is DOUBLE COMPLEX
! 50: !> The scalar b.
! 51: !> \endverbatim
! 52: !>
! 53: !> \param[out] C
! 54: !> \verbatim
! 55: !> C is DOUBLE PRECISION
! 56: !> The scalar c.
! 57: !> \endverbatim
! 58: !>
! 59: !> \param[out] S
! 60: !> \verbatim
! 61: !> S is DOUBLE COMPLEX
! 62: !> The scalar s.
! 63: !> \endverbatim
! 64: !
! 65: ! Authors:
! 66: ! ========
! 67: !
! 68: !> \author Weslley Pereira, University of Colorado Denver, USA
! 69: !
! 70: !> \date December 2021
! 71: !
! 72: !> \ingroup single_blas_level1
! 73: !
! 74: !> \par Further Details:
! 75: ! =====================
! 76: !>
! 77: !> \verbatim
! 78: !>
! 79: !> Based on the algorithm from
! 80: !>
! 81: !> Anderson E. (2017)
! 82: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
! 83: !> ACM Trans Math Softw 44:1--28
! 84: !> https://doi.org/10.1145/3061665
! 85: !>
! 86: !> \endverbatim
! 87: !
! 88: ! =====================================================================
! 89: subroutine ZROTG( a, b, c, s )
! 90: integer, parameter :: wp = kind(1.d0)
! 91: !
! 92: ! -- Reference BLAS level1 routine --
! 93: ! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
! 94: ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 95: !
! 96: ! .. Constants ..
! 97: real(wp), parameter :: zero = 0.0_wp
! 98: real(wp), parameter :: one = 1.0_wp
! 99: complex(wp), parameter :: czero = 0.0_wp
! 100: ! ..
! 101: ! .. Scaling constants ..
! 102: real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( &
! 103: minexponent(0._wp)-1, &
! 104: 1-maxexponent(0._wp) &
! 105: )
! 106: real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( &
! 107: 1-minexponent(0._wp), &
! 108: maxexponent(0._wp)-1 &
! 109: )
! 110: real(wp), parameter :: rtmin = sqrt( safmin )
! 111: ! ..
! 112: ! .. Scalar Arguments ..
! 113: real(wp) :: c
! 114: complex(wp) :: a, b, s
! 115: ! ..
! 116: ! .. Local Scalars ..
! 117: real(wp) :: d, f1, f2, g1, g2, h2, u, v, w, rtmax
! 118: complex(wp) :: f, fs, g, gs, r, t
! 119: ! ..
! 120: ! .. Intrinsic Functions ..
! 121: intrinsic :: abs, aimag, conjg, max, min, real, sqrt
! 122: ! ..
! 123: ! .. Statement Functions ..
! 124: real(wp) :: ABSSQ
! 125: ! ..
! 126: ! .. Statement Function definitions ..
! 127: ABSSQ( t ) = real( t )**2 + aimag( t )**2
! 128: ! ..
! 129: ! .. Executable Statements ..
! 130: !
! 131: f = a
! 132: g = b
! 133: if( g == czero ) then
! 134: c = one
! 135: s = czero
! 136: r = f
! 137: else if( f == czero ) then
! 138: c = zero
! 139: if( real(g) == zero ) then
! 140: r = abs(aimag(g))
! 141: s = conjg( g ) / r
! 142: elseif( aimag(g) == zero ) then
! 143: r = abs(real(g))
! 144: s = conjg( g ) / r
! 145: else
! 146: g1 = max( abs(real(g)), abs(aimag(g)) )
! 147: rtmax = sqrt( safmax/2 )
! 148: if( g1 > rtmin .and. g1 < rtmax ) then
! 149: !
! 150: ! Use unscaled algorithm
! 151: !
! 152: ! The following two lines can be replaced by `d = abs( g )`.
! 153: ! This algorithm do not use the intrinsic complex abs.
! 154: g2 = ABSSQ( g )
! 155: d = sqrt( g2 )
! 156: s = conjg( g ) / d
! 157: r = d
! 158: else
! 159: !
! 160: ! Use scaled algorithm
! 161: !
! 162: u = min( safmax, max( safmin, g1 ) )
! 163: gs = g / u
! 164: ! The following two lines can be replaced by `d = abs( gs )`.
! 165: ! This algorithm do not use the intrinsic complex abs.
! 166: g2 = ABSSQ( gs )
! 167: d = sqrt( g2 )
! 168: s = conjg( gs ) / d
! 169: r = d*u
! 170: end if
! 171: end if
! 172: else
! 173: f1 = max( abs(real(f)), abs(aimag(f)) )
! 174: g1 = max( abs(real(g)), abs(aimag(g)) )
! 175: rtmax = sqrt( safmax/4 )
! 176: if( f1 > rtmin .and. f1 < rtmax .and. &
! 177: g1 > rtmin .and. g1 < rtmax ) then
! 178: !
! 179: ! Use unscaled algorithm
! 180: !
! 181: f2 = ABSSQ( f )
! 182: g2 = ABSSQ( g )
! 183: h2 = f2 + g2
! 184: ! safmin <= f2 <= h2 <= safmax
! 185: if( f2 >= h2 * safmin ) then
! 186: ! safmin <= f2/h2 <= 1, and h2/f2 is finite
! 187: c = sqrt( f2 / h2 )
! 188: r = f / c
! 189: rtmax = rtmax * 2
! 190: if( f2 > rtmin .and. h2 < rtmax ) then
! 191: ! safmin <= sqrt( f2*h2 ) <= safmax
! 192: s = conjg( g ) * ( f / sqrt( f2*h2 ) )
! 193: else
! 194: s = conjg( g ) * ( r / h2 )
! 195: end if
! 196: else
! 197: ! f2/h2 <= safmin may be subnormal, and h2/f2 may overflow.
! 198: ! Moreover,
! 199: ! safmin <= f2*f2 * safmax < f2 * h2 < h2*h2 * safmin <= safmax,
! 200: ! sqrt(safmin) <= sqrt(f2 * h2) <= sqrt(safmax).
! 201: ! Also,
! 202: ! g2 >> f2, which means that h2 = g2.
! 203: d = sqrt( f2 * h2 )
! 204: c = f2 / d
! 205: if( c >= safmin ) then
! 206: r = f / c
! 207: else
! 208: ! f2 / sqrt(f2 * h2) < safmin, then
! 209: ! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
! 210: r = f * ( h2 / d )
! 211: end if
! 212: s = conjg( g ) * ( f / d )
! 213: end if
! 214: else
! 215: !
! 216: ! Use scaled algorithm
! 217: !
! 218: u = min( safmax, max( safmin, f1, g1 ) )
! 219: gs = g / u
! 220: g2 = ABSSQ( gs )
! 221: if( f1 / u < rtmin ) then
! 222: !
! 223: ! f is not well-scaled when scaled by g1.
! 224: ! Use a different scaling for f.
! 225: !
! 226: v = min( safmax, max( safmin, f1 ) )
! 227: w = v / u
! 228: fs = f / v
! 229: f2 = ABSSQ( fs )
! 230: h2 = f2*w**2 + g2
! 231: else
! 232: !
! 233: ! Otherwise use the same scaling for f and g.
! 234: !
! 235: w = one
! 236: fs = f / u
! 237: f2 = ABSSQ( fs )
! 238: h2 = f2 + g2
! 239: end if
! 240: ! safmin <= f2 <= h2 <= safmax
! 241: if( f2 >= h2 * safmin ) then
! 242: ! safmin <= f2/h2 <= 1, and h2/f2 is finite
! 243: c = sqrt( f2 / h2 )
! 244: r = fs / c
! 245: rtmax = rtmax * 2
! 246: if( f2 > rtmin .and. h2 < rtmax ) then
! 247: ! safmin <= sqrt( f2*h2 ) <= safmax
! 248: s = conjg( gs ) * ( fs / sqrt( f2*h2 ) )
! 249: else
! 250: s = conjg( gs ) * ( r / h2 )
! 251: end if
! 252: else
! 253: ! f2/h2 <= safmin may be subnormal, and h2/f2 may overflow.
! 254: ! Moreover,
! 255: ! safmin <= f2*f2 * safmax < f2 * h2 < h2*h2 * safmin <= safmax,
! 256: ! sqrt(safmin) <= sqrt(f2 * h2) <= sqrt(safmax).
! 257: ! Also,
! 258: ! g2 >> f2, which means that h2 = g2.
! 259: d = sqrt( f2 * h2 )
! 260: c = f2 / d
! 261: if( c >= safmin ) then
! 262: r = fs / c
! 263: else
! 264: ! f2 / sqrt(f2 * h2) < safmin, then
! 265: ! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
! 266: r = fs * ( h2 / d )
! 267: end if
! 268: s = conjg( gs ) * ( fs / d )
! 269: end if
! 270: ! Rescale c and r
! 271: c = c * w
! 272: r = r * u
! 273: end if
! 274: end if
! 275: a = r
! 276: return
! 277: end subroutine
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zrotg.f90 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas