LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ zrotg()

 subroutine zrotg ( complex(wp) a, complex(wp) b, real(wp) c, complex(wp) s )

ZROTG generates a Givens rotation with real cosine and complex sine.

Purpose:
``` The computation uses the formulas
|x| = sqrt( Re(x)**2 + Im(x)**2 )
sgn(x) = x / |x|  if x /= 0
= 1        if x  = 0
c = |a| / sqrt(|a|**2 + |b|**2)
s = sgn(a) * conjg(b) / sqrt(|a|**2 + |b|**2)
r = sgn(a)*sqrt(|a|**2 + |b|**2)
When a and b are real and r /= 0, the formulas simplify to
c = a / r
s = b / r
the same as in DROTG when |a| > |b|.  When |b| >= |a|, the
sign of c and s will be different from those computed by DROTG
if the signs of a and b are not the same.```
Parameters
 [in,out] A ``` A is DOUBLE COMPLEX On entry, the scalar a. On exit, the scalar r.``` [in] B ``` B is DOUBLE COMPLEX The scalar b.``` [out] C ``` C is DOUBLE PRECISION The scalar c.``` [out] S ``` S is DOUBLE COMPLEX The scalar s.```
Date
December 2021
Further Details:
``` Based on the algorithm from

Anderson E. (2017)
Algorithm 978: Safe Scaling in the Level 1 BLAS
ACM Trans Math Softw 44:1--28
https://doi.org/10.1145/3061665```

Definition at line 89 of file zrotg.f90.

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
Here is the caller graph for this function: