LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
crotg.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 ! CROTG constructs a plane rotation
12 ! [ c s ] [ a ] = [ r ]
13 ! [ -conjg(s) c ] [ b ] [ 0 ]
14 ! where c is real, s ic complex, and c**2 + conjg(s)*s = 1.
15 !
17 ! =============
36 !
37 ! Arguments:
38 ! ==========
39 !
64 !
65 ! Authors:
66 ! ========
67 !
69 !
71 ! ==================
74 !
76 !
78 ! =====================
88 !
89 ! =====================================================================
90 subroutine crotg( a, b, c, s )
91  integer, parameter :: wp = kind(1.e0)
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
229 end subroutine
subroutine crotg(a, b, c, s)
CROTG
Definition: crotg.f90:91