1!> \brief \b CROTG 2! 3! =========== DOCUMENTATION =========== 4! 5! Online html documentation available at 6! http://www.netlib.org/lapack/explore-html/ 7! 8! Definition: 9! =========== 10! 11! CROTG constructs a plane rotation 12! [ c s ] [ a ] = [ r ] 13! [ -conjg(s) c ] [ b ] [ 0 ] 14! where c is real, s ic 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!> When a and b are real and r /= 0, the formulas simplify to 28!> r = sgn(a)*sqrt(|a|**2 + |b|**2) 29!> c = a / r 30!> s = b / r 31!> the same as in CROTG when |a| > |b|. When |b| >= |a|, the 32!> sign of c and s will be different from those computed by CROTG 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 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 COMPLEX 50!> The scalar b. 51!> \endverbatim 52!> 53!> \param[out] C 54!> \verbatim 55!> C is REAL 56!> The scalar c. 57!> \endverbatim 58!> 59!> \param[out] S 60!> \verbatim 61!> S is REAL 62!> The scalar s. 63!> \endverbatim 64! 65! Authors: 66! ======== 67! 68!> \author Edward Anderson, Lockheed Martin 69! 70!> \par Contributors: 71! ================== 72!> 73!> Weslley Pereira, University of Colorado Denver, USA 74! 75!> \ingroup single_blas_level1 76! 77!> \par Further Details: 78! ===================== 79!> 80!> \verbatim 81!> 82!> Anderson E. (2017) 83!> Algorithm 978: Safe Scaling in the Level 1 BLAS 84!> ACM Trans Math Softw 44:1--28 85!> https://doi.org/10.1145/3061665 86!> 87!> \endverbatim 88! 89! ===================================================================== 90subroutine CROTG( a, b, c, s ) 91 integer, parameter :: wp = kind(1.e0) 92! 93! -- Reference BLAS level1 routine -- 94! -- Reference BLAS is a software package provided by Univ. of Tennessee, -- 95! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 96! 97! .. Constants .. 98 real(wp), parameter :: zero = 0.0_wp 99 real(wp), parameter :: one = 1.0_wp 100 complex(wp), parameter :: czero = 0.0_wp 101! .. 102! .. Scaling constants .. 103 real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( & 104 minexponent(0._wp)-1, & 105 1-maxexponent(0._wp) & 106 ) 107 real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( & 108 1-minexponent(0._wp), & 109 maxexponent(0._wp)-1 & 110 ) 111 real(wp), parameter :: rtmin = sqrt( real(radix(0._wp),wp)**max( & 112 minexponent(0._wp)-1, & 113 1-maxexponent(0._wp) & 114 ) / epsilon(0._wp) ) 115 real(wp), parameter :: rtmax = sqrt( real(radix(0._wp),wp)**max( & 116 1-minexponent(0._wp), & 117 maxexponent(0._wp)-1 & 118 ) * epsilon(0._wp) ) 119! .. 120! .. Scalar Arguments .. 121 real(wp) :: c 122 complex(wp) :: a, b, s 123! .. 124! .. Local Scalars .. 125 real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w 126 complex(wp) :: f, fs, g, gs, r, t 127! .. 128! .. Intrinsic Functions .. 129 intrinsic :: abs, aimag, conjg, max, min, real, sqrt 130! .. 131! .. Statement Functions .. 132 real(wp) :: ABSSQ 133! .. 134! .. Statement Function definitions .. 135 ABSSQ( t ) = real( t )**2 + aimag( t )**2 136! .. 137! .. Executable Statements .. 138! 139 f = a 140 g = b 141 if( g == czero ) then 142 c = one 143 s = czero 144 r = f 145 else if( f == czero ) then 146 c = zero 147 g1 = max( abs(real(g)), abs(aimag(g)) ) 148 if( g1 > rtmin .and. g1 < rtmax ) then 149! 150! Use unscaled algorithm 151! 152 g2 = ABSSQ( g ) 153 d = sqrt( g2 ) 154 s = conjg( g ) / d 155 r = d 156 else 157! 158! Use scaled algorithm 159! 160 u = min( safmax, max( safmin, g1 ) ) 161 uu = one / u 162 gs = g*uu 163 g2 = ABSSQ( gs ) 164 d = sqrt( g2 ) 165 s = conjg( gs ) / d 166 r = d*u 167 end if 168 else 169 f1 = max( abs(real(f)), abs(aimag(f)) ) 170 g1 = max( abs(real(g)), abs(aimag(g)) ) 171 if( f1 > rtmin .and. f1 < rtmax .and. & 172 g1 > rtmin .and. g1 < rtmax ) then 173! 174! Use unscaled algorithm 175! 176 f2 = ABSSQ( f ) 177 g2 = ABSSQ( g ) 178 h2 = f2 + g2 179 if( f2 > rtmin .and. h2 < rtmax ) then 180 d = sqrt( f2*h2 ) 181 else 182 d = sqrt( f2 )*sqrt( h2 ) 183 end if 184 p = 1 / d 185 c = f2*p 186 s = conjg( g )*( f*p ) 187 r = f*( h2*p ) 188 else 189! 190! Use scaled algorithm 191! 192 u = min( safmax, max( safmin, f1, g1 ) ) 193 uu = one / u 194 gs = g*uu 195 g2 = ABSSQ( gs ) 196 if( f1*uu < rtmin ) then 197! 198! f is not well-scaled when scaled by g1. 199! Use a different scaling for f. 200! 201 v = min( safmax, max( safmin, f1 ) ) 202 vv = one / v 203 w = v * uu 204 fs = f*vv 205 f2 = ABSSQ( fs ) 206 h2 = f2*w**2 + g2 207 else 208! 209! Otherwise use the same scaling for f and g. 210! 211 w = one 212 fs = f*uu 213 f2 = ABSSQ( fs ) 214 h2 = f2 + g2 215 end if 216 if( f2 > rtmin .and. h2 < rtmax ) then 217 d = sqrt( f2*h2 ) 218 else 219 d = sqrt( f2 )*sqrt( h2 ) 220 end if 221 p = 1 / d 222 c = ( f2*p )*w 223 s = conjg( gs )*( fs*p ) 224 r = ( fs*( h2*p ) )*u 225 end if 226 end if 227 a = r 228 return 229end subroutine 230