1!> \brief \b CLARTG 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 CLARTG( 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!> CLARTG 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!> When F and G are real, the formulas simplify to C = F/R and 42!> S = G/R, and the returned values of C, S, and R should be 43!> identical to those returned by CLARTG. 44!> 45!> The algorithm used to compute these quantities incorporates scaling 46!> to avoid overflow or underflow in computing the square root of the 47!> sum of squares. 48!> 49!> This is a faster version of the BLAS1 routine CROTG, except for 50!> the following differences: 51!> F and G are unchanged on return. 52!> If G=0, then C=1 and S=0. 53!> If F=0, then C=0 and S is chosen so that R is real. 54!> 55!> Below, wp=>sp stands for single precision from LA_CONSTANTS module. 56!> \endverbatim 57! 58! Arguments: 59! ========== 60! 61!> \param[in] F 62!> \verbatim 63!> F is COMPLEX(wp) 64!> The first component of vector to be rotated. 65!> \endverbatim 66!> 67!> \param[in] G 68!> \verbatim 69!> G is COMPLEX(wp) 70!> The second component of vector to be rotated. 71!> \endverbatim 72!> 73!> \param[out] C 74!> \verbatim 75!> C is REAL(wp) 76!> The cosine of the rotation. 77!> \endverbatim 78!> 79!> \param[out] S 80!> \verbatim 81!> S is COMPLEX(wp) 82!> The sine of the rotation. 83!> \endverbatim 84!> 85!> \param[out] R 86!> \verbatim 87!> R is COMPLEX(wp) 88!> The nonzero component of the rotated vector. 89!> \endverbatim 90! 91! Authors: 92! ======== 93! 94!> \author Edward Anderson, Lockheed Martin 95! 96!> \date August 2016 97! 98!> \ingroup OTHERauxiliary 99! 100!> \par Contributors: 101! ================== 102!> 103!> Weslley Pereira, University of Colorado Denver, USA 104! 105!> \par Further Details: 106! ===================== 107!> 108!> \verbatim 109!> 110!> Anderson E. (2017) 111!> Algorithm 978: Safe Scaling in the Level 1 BLAS 112!> ACM Trans Math Softw 44:1--28 113!> https://doi.org/10.1145/3061665 114!> 115!> \endverbatim 116! 117subroutine CLARTG( f, g, c, s, r ) 118 use LA_CONSTANTS, & 119 only: wp=>sp, zero=>szero, one=>sone, two=>stwo, czero, & 120 rtmin=>srtmin, rtmax=>srtmax, safmin=>ssafmin, safmax=>ssafmax 121! 122! -- LAPACK auxiliary routine (version 3.10.0) -- 123! -- LAPACK is a software package provided by Univ. of Tennessee, -- 124! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 125! February 2021 126! 127! .. Scalar Arguments .. 128 real(wp) c 129 complex(wp) f, g, r, s 130! .. 131! .. Local Scalars .. 132 real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w 133 complex(wp) :: fs, gs, t 134! .. 135! .. Intrinsic Functions .. 136 intrinsic :: abs, aimag, conjg, max, min, real, sqrt 137! .. 138! .. Statement Functions .. 139 real(wp) :: ABSSQ 140! .. 141! .. Statement Function definitions .. 142 ABSSQ( t ) = real( t )**2 + aimag( t )**2 143! .. 144! .. Executable Statements .. 145! 146 if( g == czero ) then 147 c = one 148 s = czero 149 r = f 150 else if( f == czero ) then 151 c = zero 152 g1 = max( abs(real(g)), abs(aimag(g)) ) 153 if( g1 > rtmin .and. g1 < rtmax ) then 154! 155! Use unscaled algorithm 156! 157 g2 = ABSSQ( g ) 158 d = sqrt( g2 ) 159 s = conjg( g ) / d 160 r = d 161 else 162! 163! Use scaled algorithm 164! 165 u = min( safmax, max( safmin, g1 ) ) 166 uu = one / u 167 gs = g*uu 168 g2 = ABSSQ( gs ) 169 d = sqrt( g2 ) 170 s = conjg( gs ) / d 171 r = d*u 172 end if 173 else 174 f1 = max( abs(real(f)), abs(aimag(f)) ) 175 g1 = max( abs(real(g)), abs(aimag(g)) ) 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 if( f2 > rtmin .and. h2 < rtmax ) then 185 d = sqrt( f2*h2 ) 186 else 187 d = sqrt( f2 )*sqrt( h2 ) 188 end if 189 p = 1 / d 190 c = f2*p 191 s = conjg( g )*( f*p ) 192 r = f*( h2*p ) 193 else 194! 195! Use scaled algorithm 196! 197 u = min( safmax, max( safmin, f1, g1 ) ) 198 uu = one / u 199 gs = g*uu 200 g2 = ABSSQ( gs ) 201 if( f1*uu < rtmin ) then 202! 203! f is not well-scaled when scaled by g1. 204! Use a different scaling for f. 205! 206 v = min( safmax, max( safmin, f1 ) ) 207 vv = one / v 208 w = v * uu 209 fs = f*vv 210 f2 = ABSSQ( fs ) 211 h2 = f2*w**2 + g2 212 else 213! 214! Otherwise use the same scaling for f and g. 215! 216 w = one 217 fs = f*uu 218 f2 = ABSSQ( fs ) 219 h2 = f2 + g2 220 end if 221 if( f2 > rtmin .and. h2 < rtmax ) then 222 d = sqrt( f2*h2 ) 223 else 224 d = sqrt( f2 )*sqrt( h2 ) 225 end if 226 p = 1 / d 227 c = ( f2*p )*w 228 s = conjg( gs )*( fs*p ) 229 r = ( fs*( h2*p ) )*u 230 end if 231 end if 232 return 233end subroutine 234