Annotation of rpl/lapack/lapack/zlartg.f90, revision 1.1
1.1 ! bertrand 1: !> \brief \b ZLARTG generates a plane 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: ! SUBROUTINE ZLARTG( F, G, C, S, R )
! 12: !
! 13: ! .. Scalar Arguments ..
! 14: ! REAL(wp) C
! 15: ! COMPLEX(wp) F, G, R, S
! 16: ! ..
! 17: !
! 18: !> \par Purpose:
! 19: ! =============
! 20: !>
! 21: !> \verbatim
! 22: !>
! 23: !> ZLARTG generates a plane rotation so that
! 24: !>
! 25: !> [ C S ] . [ F ] = [ R ]
! 26: !> [ -conjg(S) C ] [ G ] [ 0 ]
! 27: !>
! 28: !> where C is real and C**2 + |S|**2 = 1.
! 29: !>
! 30: !> The mathematical formulas used for C and S are
! 31: !>
! 32: !> sgn(x) = { x / |x|, x != 0
! 33: !> { 1, x = 0
! 34: !>
! 35: !> R = sgn(F) * sqrt(|F|**2 + |G|**2)
! 36: !>
! 37: !> C = |F| / sqrt(|F|**2 + |G|**2)
! 38: !>
! 39: !> S = sgn(F) * conjg(G) / sqrt(|F|**2 + |G|**2)
! 40: !>
! 41: !> Special conditions:
! 42: !> If G=0, then C=1 and S=0.
! 43: !> If F=0, then C=0 and S is chosen so that R is real.
! 44: !>
! 45: !> When F and G are real, the formulas simplify to C = F/R and
! 46: !> S = G/R, and the returned values of C, S, and R should be
! 47: !> identical to those returned by DLARTG.
! 48: !>
! 49: !> The algorithm used to compute these quantities incorporates scaling
! 50: !> to avoid overflow or underflow in computing the square root of the
! 51: !> sum of squares.
! 52: !>
! 53: !> This is the same routine ZROTG fom BLAS1, except that
! 54: !> F and G are unchanged on return.
! 55: !>
! 56: !> Below, wp=>dp stands for double precision from LA_CONSTANTS module.
! 57: !> \endverbatim
! 58: !
! 59: ! Arguments:
! 60: ! ==========
! 61: !
! 62: !> \param[in] F
! 63: !> \verbatim
! 64: !> F is COMPLEX(wp)
! 65: !> The first component of vector to be rotated.
! 66: !> \endverbatim
! 67: !>
! 68: !> \param[in] G
! 69: !> \verbatim
! 70: !> G is COMPLEX(wp)
! 71: !> The second component of vector to be rotated.
! 72: !> \endverbatim
! 73: !>
! 74: !> \param[out] C
! 75: !> \verbatim
! 76: !> C is REAL(wp)
! 77: !> The cosine of the rotation.
! 78: !> \endverbatim
! 79: !>
! 80: !> \param[out] S
! 81: !> \verbatim
! 82: !> S is COMPLEX(wp)
! 83: !> The sine of the rotation.
! 84: !> \endverbatim
! 85: !>
! 86: !> \param[out] R
! 87: !> \verbatim
! 88: !> R is COMPLEX(wp)
! 89: !> The nonzero component of the rotated vector.
! 90: !> \endverbatim
! 91: !
! 92: ! Authors:
! 93: ! ========
! 94: !
! 95: !> \author Weslley Pereira, University of Colorado Denver, USA
! 96: !
! 97: !> \date December 2021
! 98: !
! 99: !> \ingroup OTHERauxiliary
! 100: !
! 101: !> \par Further Details:
! 102: ! =====================
! 103: !>
! 104: !> \verbatim
! 105: !>
! 106: !> Based on the algorithm from
! 107: !>
! 108: !> Anderson E. (2017)
! 109: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
! 110: !> ACM Trans Math Softw 44:1--28
! 111: !> https://doi.org/10.1145/3061665
! 112: !>
! 113: !> \endverbatim
! 114: !
! 115: subroutine ZLARTG( f, g, c, s, r )
! 116: use LA_CONSTANTS, &
! 117: only: wp=>dp, zero=>dzero, one=>done, two=>dtwo, czero=>zzero, &
! 118: safmin=>dsafmin, safmax=>dsafmax
! 119: !
! 120: ! -- LAPACK auxiliary routine --
! 121: ! -- LAPACK is a software package provided by Univ. of Tennessee, --
! 122: ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 123: ! February 2021
! 124: !
! 125: ! .. Scalar Arguments ..
! 126: real(wp) c
! 127: complex(wp) f, g, r, s
! 128: ! ..
! 129: ! .. Local Scalars ..
! 130: real(wp) :: d, f1, f2, g1, g2, h2, u, v, w, rtmin, rtmax
! 131: complex(wp) :: fs, gs, t
! 132: ! ..
! 133: ! .. Intrinsic Functions ..
! 134: intrinsic :: abs, aimag, conjg, max, min, real, sqrt
! 135: ! ..
! 136: ! .. Statement Functions ..
! 137: real(wp) :: ABSSQ
! 138: ! ..
! 139: ! .. Statement Function definitions ..
! 140: ABSSQ( t ) = real( t )**2 + aimag( t )**2
! 141: ! ..
! 142: ! .. Constants ..
! 143: rtmin = sqrt( safmin )
! 144: ! ..
! 145: ! .. Executable Statements ..
! 146: !
! 147: if( g == czero ) then
! 148: c = one
! 149: s = czero
! 150: r = f
! 151: else if( f == czero ) then
! 152: c = zero
! 153: if( real(g) == zero ) then
! 154: r = abs(aimag(g))
! 155: s = conjg( g ) / r
! 156: elseif( aimag(g) == zero ) then
! 157: r = abs(real(g))
! 158: s = conjg( g ) / r
! 159: else
! 160: g1 = max( abs(real(g)), abs(aimag(g)) )
! 161: rtmax = sqrt( safmax/2 )
! 162: if( g1 > rtmin .and. g1 < rtmax ) then
! 163: !
! 164: ! Use unscaled algorithm
! 165: !
! 166: ! The following two lines can be replaced by `d = abs( g )`.
! 167: ! This algorithm do not use the intrinsic complex abs.
! 168: g2 = ABSSQ( g )
! 169: d = sqrt( g2 )
! 170: s = conjg( g ) / d
! 171: r = d
! 172: else
! 173: !
! 174: ! Use scaled algorithm
! 175: !
! 176: u = min( safmax, max( safmin, g1 ) )
! 177: gs = g / u
! 178: ! The following two lines can be replaced by `d = abs( gs )`.
! 179: ! This algorithm do not use the intrinsic complex abs.
! 180: g2 = ABSSQ( gs )
! 181: d = sqrt( g2 )
! 182: s = conjg( gs ) / d
! 183: r = d*u
! 184: end if
! 185: end if
! 186: else
! 187: f1 = max( abs(real(f)), abs(aimag(f)) )
! 188: g1 = max( abs(real(g)), abs(aimag(g)) )
! 189: rtmax = sqrt( safmax/4 )
! 190: if( f1 > rtmin .and. f1 < rtmax .and. &
! 191: g1 > rtmin .and. g1 < rtmax ) then
! 192: !
! 193: ! Use unscaled algorithm
! 194: !
! 195: f2 = ABSSQ( f )
! 196: g2 = ABSSQ( g )
! 197: h2 = f2 + g2
! 198: ! safmin <= f2 <= h2 <= safmax
! 199: if( f2 >= h2 * safmin ) then
! 200: ! safmin <= f2/h2 <= 1, and h2/f2 is finite
! 201: c = sqrt( f2 / h2 )
! 202: r = f / c
! 203: rtmax = rtmax * 2
! 204: if( f2 > rtmin .and. h2 < rtmax ) then
! 205: ! safmin <= sqrt( f2*h2 ) <= safmax
! 206: s = conjg( g ) * ( f / sqrt( f2*h2 ) )
! 207: else
! 208: s = conjg( g ) * ( r / h2 )
! 209: end if
! 210: else
! 211: ! f2/h2 <= safmin may be subnormal, and h2/f2 may overflow.
! 212: ! Moreover,
! 213: ! safmin <= f2*f2 * safmax < f2 * h2 < h2*h2 * safmin <= safmax,
! 214: ! sqrt(safmin) <= sqrt(f2 * h2) <= sqrt(safmax).
! 215: ! Also,
! 216: ! g2 >> f2, which means that h2 = g2.
! 217: d = sqrt( f2 * h2 )
! 218: c = f2 / d
! 219: if( c >= safmin ) then
! 220: r = f / c
! 221: else
! 222: ! f2 / sqrt(f2 * h2) < safmin, then
! 223: ! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
! 224: r = f * ( h2 / d )
! 225: end if
! 226: s = conjg( g ) * ( f / d )
! 227: end if
! 228: else
! 229: !
! 230: ! Use scaled algorithm
! 231: !
! 232: u = min( safmax, max( safmin, f1, g1 ) )
! 233: gs = g / u
! 234: g2 = ABSSQ( gs )
! 235: if( f1 / u < rtmin ) then
! 236: !
! 237: ! f is not well-scaled when scaled by g1.
! 238: ! Use a different scaling for f.
! 239: !
! 240: v = min( safmax, max( safmin, f1 ) )
! 241: w = v / u
! 242: fs = f / v
! 243: f2 = ABSSQ( fs )
! 244: h2 = f2*w**2 + g2
! 245: else
! 246: !
! 247: ! Otherwise use the same scaling for f and g.
! 248: !
! 249: w = one
! 250: fs = f / u
! 251: f2 = ABSSQ( fs )
! 252: h2 = f2 + g2
! 253: end if
! 254: ! safmin <= f2 <= h2 <= safmax
! 255: if( f2 >= h2 * safmin ) then
! 256: ! safmin <= f2/h2 <= 1, and h2/f2 is finite
! 257: c = sqrt( f2 / h2 )
! 258: r = fs / c
! 259: rtmax = rtmax * 2
! 260: if( f2 > rtmin .and. h2 < rtmax ) then
! 261: ! safmin <= sqrt( f2*h2 ) <= safmax
! 262: s = conjg( gs ) * ( fs / sqrt( f2*h2 ) )
! 263: else
! 264: s = conjg( gs ) * ( r / h2 )
! 265: end if
! 266: else
! 267: ! f2/h2 <= safmin may be subnormal, and h2/f2 may overflow.
! 268: ! Moreover,
! 269: ! safmin <= f2*f2 * safmax < f2 * h2 < h2*h2 * safmin <= safmax,
! 270: ! sqrt(safmin) <= sqrt(f2 * h2) <= sqrt(safmax).
! 271: ! Also,
! 272: ! g2 >> f2, which means that h2 = g2.
! 273: d = sqrt( f2 * h2 )
! 274: c = f2 / d
! 275: if( c >= safmin ) then
! 276: r = fs / c
! 277: else
! 278: ! f2 / sqrt(f2 * h2) < safmin, then
! 279: ! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
! 280: r = fs * ( h2 / d )
! 281: end if
! 282: s = conjg( gs ) * ( fs / d )
! 283: end if
! 284: ! Rescale c and r
! 285: c = c * w
! 286: r = r * u
! 287: end if
! 288: end if
! 289: return
! 290: end subroutine
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlartg.f90 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack