1: !> \brief \b ZROTG generates a Givens 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: ! ZROTG constructs a plane rotation
12: ! [ c s ] [ a ] = [ r ]
13: ! [ -conjg(s) c ] [ b ] [ 0 ]
14: ! where c is real, s is 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: !> r = sgn(a)*sqrt(|a|**2 + |b|**2)
28: !> When a and b are real and r /= 0, the formulas simplify to
29: !> c = a / r
30: !> s = b / r
31: !> the same as in DROTG when |a| > |b|. When |b| >= |a|, the
32: !> sign of c and s will be different from those computed by DROTG
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 DOUBLE 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 DOUBLE COMPLEX
50: !> The scalar b.
51: !> \endverbatim
52: !>
53: !> \param[out] C
54: !> \verbatim
55: !> C is DOUBLE PRECISION
56: !> The scalar c.
57: !> \endverbatim
58: !>
59: !> \param[out] S
60: !> \verbatim
61: !> S is DOUBLE COMPLEX
62: !> The scalar s.
63: !> \endverbatim
64: !
65: ! Authors:
66: ! ========
67: !
68: !> \author Weslley Pereira, University of Colorado Denver, USA
69: !
70: !> \date December 2021
71: !
72: !> \ingroup single_blas_level1
73: !
74: !> \par Further Details:
75: ! =====================
76: !>
77: !> \verbatim
78: !>
79: !> Based on the algorithm from
80: !>
81: !> Anderson E. (2017)
82: !> Algorithm 978: Safe Scaling in the Level 1 BLAS
83: !> ACM Trans Math Softw 44:1--28
84: !> https://doi.org/10.1145/3061665
85: !>
86: !> \endverbatim
87: !
88: ! =====================================================================
89: subroutine ZROTG( a, b, c, s )
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
277: end subroutine
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zrotg.f90 CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / blas