1: !> \brief \b ZLARTG 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 ZLARTG( 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: !> ZLARTG 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: !> Special conditions:
42: !> If G=0, then C=1 and S=0.
43: !> If F=0, then C=0 and S is chosen so that R is real.
44: !>
45: !> When F and G are real, the formulas simplify to C = F/R and
46: !> S = G/R, and the returned values of C, S, and R should be
47: !> identical to those returned by DLARTG.
48: !>
49: !> The algorithm used to compute these quantities incorporates scaling
50: !> to avoid overflow or underflow in computing the square root of the
51: !> sum of squares.
52: !>
53: !> This is the same routine ZROTG fom BLAS1, except that
54: !> F and G are unchanged on return.
55: !>
56: !> Below, wp=>dp stands for double precision from LA_CONSTANTS module.
57: !> \endverbatim
58: !
59: ! Arguments:
60: ! ==========
61: !
62: !> \param[in] F
63: !> \verbatim
64: !> F is COMPLEX(wp)
65: !> The first component of vector to be rotated.
66: !> \endverbatim
67: !>
68: !> \param[in] G
69: !> \verbatim
70: !> G is COMPLEX(wp)
71: !> The second component of vector to be rotated.
72: !> \endverbatim
73: !>
74: !> \param[out] C
75: !> \verbatim
76: !> C is REAL(wp)
77: !> The cosine of the rotation.
78: !> \endverbatim
79: !>
80: !> \param[out] S
81: !> \verbatim
82: !> S is COMPLEX(wp)
83: !> The sine of the rotation.
84: !> \endverbatim
85: !>
86: !> \param[out] R
87: !> \verbatim
88: !> R is COMPLEX(wp)
89: !> The nonzero component of the rotated vector.
90: !> \endverbatim
91: !
92: ! Authors:
93: ! ========
94: !
95: !> \author Weslley Pereira, University of Colorado Denver, USA
96: !
97: !> \date December 2021
98: !
99: !> \ingroup OTHERauxiliary
100: !
101: !> \par Further Details:
102: ! =====================
103: !>
104: !> \verbatim
105: !>
106: !> Based on the algorithm from
107: !>
108: !> Anderson E. (2017)
109: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
110: !> ACM Trans Math Softw 44:1--28
111: !> https://doi.org/10.1145/3061665
112: !>
113: !> \endverbatim
114: !
115: subroutine 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
290: end subroutine
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlartg.f90 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack