LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlartg.f90
Go to the documentation of this file.
1
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!
19! =============
58!
59! Arguments:
60! ==========
61!
91!
92! Authors:
93! ========
94!
96!
98!
100!
102! =====================
114!
115subroutine 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
290end subroutine
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
real(dp), parameter dtwo
real(dp), parameter dzero
real(dp), parameter dsafmin
integer, parameter dp
real(dp), parameter done
complex(sp), parameter czero
complex(dp), parameter zzero
real(dp), parameter dsafmax
LA_CONSTANTS is a module for the scaling constants for the compiled Fortran single and double precisi...