LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

## ◆ zrotg()

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

ZROTG

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)
When a and b are real and r /= 0, the formulas simplify to
r = sgn(a)*sqrt(|a|**2 + |b|**2)
c = a / r
s = b / r
the same as in ZROTG when |a| > |b|.  When |b| >= |a|, the
sign of c and s will be different from those computed by ZROTG
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 PRECISION The scalar s.
Contributors:
Weslley Pereira, University of Colorado Denver, USA
Further Details:
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 90 of file zrotg.f90.

91  integer, parameter :: wp = kind(1.d0)
92 !
93 ! -- Reference BLAS level1 routine --
94 ! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
95 ! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
96 !
97 ! .. Constants ..
98  real(wp), parameter :: zero = 0.0_wp
99  real(wp), parameter :: one = 1.0_wp
100  complex(wp), parameter :: czero = 0.0_wp
101 ! ..
102 ! .. Scaling constants ..
103  real(wp), parameter :: safmin = real(radix(0._wp),wp)**max( &
104  minexponent(0._wp)-1, &
105  1-maxexponent(0._wp) &
106  )
107  real(wp), parameter :: safmax = real(radix(0._wp),wp)**max( &
108  1-minexponent(0._wp), &
109  maxexponent(0._wp)-1 &
110  )
111  real(wp), parameter :: rtmin = sqrt( real(radix(0._wp),wp)**max( &
112  minexponent(0._wp)-1, &
113  1-maxexponent(0._wp) &
114  ) / epsilon(0._wp) )
115  real(wp), parameter :: rtmax = sqrt( real(radix(0._wp),wp)**max( &
116  1-minexponent(0._wp), &
117  maxexponent(0._wp)-1 &
118  ) * epsilon(0._wp) )
119 ! ..
120 ! .. Scalar Arguments ..
121  real(wp) :: c
122  complex(wp) :: a, b, s
123 ! ..
124 ! .. Local Scalars ..
125  real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
126  complex(wp) :: f, fs, g, gs, r, t
127 ! ..
128 ! .. Intrinsic Functions ..
129  intrinsic :: abs, aimag, conjg, max, min, real, sqrt
130 ! ..
131 ! .. Statement Functions ..
132  real(wp) :: ABSSQ
133 ! ..
134 ! .. Statement Function definitions ..
135  abssq( t ) = real( t )**2 + aimag( t )**2
136 ! ..
137 ! .. Executable Statements ..
138 !
139  f = a
140  g = b
141  if( g == czero ) then
142  c = one
143  s = czero
144  r = f
145  else if( f == czero ) then
146  c = zero
147  g1 = max( abs(real(g)), abs(aimag(g)) )
148  if( g1 > rtmin .and. g1 < rtmax ) then
149 !
150 ! Use unscaled algorithm
151 !
152  g2 = abssq( g )
153  d = sqrt( g2 )
154  s = conjg( g ) / d
155  r = d
156  else
157 !
158 ! Use scaled algorithm
159 !
160  u = min( safmax, max( safmin, g1 ) )
161  uu = one / u
162  gs = g*uu
163  g2 = abssq( gs )
164  d = sqrt( g2 )
165  s = conjg( gs ) / d
166  r = d*u
167  end if
168  else
169  f1 = max( abs(real(f)), abs(aimag(f)) )
170  g1 = max( abs(real(g)), abs(aimag(g)) )
171  if( f1 > rtmin .and. f1 < rtmax .and. &
172  g1 > rtmin .and. g1 < rtmax ) then
173 !
174 ! Use unscaled algorithm
175 !
176  f2 = abssq( f )
177  g2 = abssq( g )
178  h2 = f2 + g2
179  if( f2 > rtmin .and. h2 < rtmax ) then
180  d = sqrt( f2*h2 )
181  else
182  d = sqrt( f2 )*sqrt( h2 )
183  end if
184  p = 1 / d
185  c = f2*p
186  s = conjg( g )*( f*p )
187  r = f*( h2*p )
188  else
189 !
190 ! Use scaled algorithm
191 !
192  u = min( safmax, max( safmin, f1, g1 ) )
193  uu = one / u
194  gs = g*uu
195  g2 = abssq( gs )
196  if( f1*uu < rtmin ) then
197 !
198 ! f is not well-scaled when scaled by g1.
199 ! Use a different scaling for f.
200 !
201  v = min( safmax, max( safmin, f1 ) )
202  vv = one / v
203  w = v * uu
204  fs = f*vv
205  f2 = abssq( fs )
206  h2 = f2*w**2 + g2
207  else
208 !
209 ! Otherwise use the same scaling for f and g.
210 !
211  w = one
212  fs = f*uu
213  f2 = abssq( fs )
214  h2 = f2 + g2
215  end if
216  if( f2 > rtmin .and. h2 < rtmax ) then
217  d = sqrt( f2*h2 )
218  else
219  d = sqrt( f2 )*sqrt( h2 )
220  end if
221  p = 1 / d
222  c = ( f2*p )*w
223  s = conjg( gs )*( fs*p )
224  r = ( fs*( h2*p ) )*u
225  end if
226  end if
227  a = r
228  return
Here is the caller graph for this function: