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

## ◆ clartg()

 subroutine clartg ( complex(wp) f, complex(wp) g, real(wp) c, complex(wp) s, complex(wp) r )

CLARTG generates a plane rotation with real cosine and complex sine.

Purpose:
``` CLARTG generates a plane rotation so that

[  C         S  ] . [ F ]  =  [ R ]
[ -conjg(S)  C  ]   [ G ]     [ 0 ]

where C is real and C**2 + |S|**2 = 1.

The mathematical formulas used for C and S are

sgn(x) = {  x / |x|,   x != 0
{  1,         x  = 0

R = sgn(F) * sqrt(|F|**2 + |G|**2)

C = |F| / sqrt(|F|**2 + |G|**2)

S = sgn(F) * conjg(G) / sqrt(|F|**2 + |G|**2)

Special conditions:
If G=0, then C=1 and S=0.
If F=0, then C=0 and S is chosen so that R is real.

When F and G are real, the formulas simplify to C = F/R and
S = G/R, and the returned values of C, S, and R should be
identical to those returned by SLARTG.

The algorithm used to compute these quantities incorporates scaling
to avoid overflow or underflow in computing the square root of the
sum of squares.

This is the same routine CROTG fom BLAS1, except that
F and G are unchanged on return.

Below, wp=>sp stands for single precision from LA_CONSTANTS module.```
Parameters
 [in] F ``` F is COMPLEX(wp) The first component of vector to be rotated.``` [in] G ``` G is COMPLEX(wp) The second component of vector to be rotated.``` [out] C ``` C is REAL(wp) The cosine of the rotation.``` [out] S ``` S is COMPLEX(wp) The sine of the rotation.``` [out] R ``` R is COMPLEX(wp) The nonzero component of the rotated vector.```
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 115 of file clartg.f90.

116 use la_constants, &
117 only: wp=>sp, zero=>szero, one=>sone, two=>stwo, czero, &
118 safmin=>ssafmin, safmax=>ssafmax
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
real(sp), parameter sone
real(sp), parameter stwo
complex(sp), parameter czero
integer, parameter sp
real(sp), parameter ssafmin
real(sp), parameter ssafmax
real(sp), parameter szero
LA_CONSTANTS is a module for the scaling constants for the compiled Fortran single and double precisi...
Here is the caller graph for this function: