@@ -122,7 +122,7 @@ subroutine CROTG( a, b, c, s )
122122 complex (wp) :: a, b, s
123123! ..
124124! .. Local Scalars ..
125- real (wp) :: d, f1, f2, g1, g2, h2, p , u, uu, v, vv, w
125+ real (wp) :: d, f1, f2, g1, g2, h2, w2 , u, uu, v, vv, w
126126 complex (wp) :: f, fs, g, gs, r, t
127127! ..
128128! .. Intrinsic Functions ..
@@ -149,8 +149,7 @@ subroutine CROTG( a, b, c, s )
149149!
150150! Use unscaled algorithm
151151!
152- g2 = ABSSQ( g )
153- d = sqrt ( g2 )
152+ d = abs ( g )
154153 s = conjg ( g ) / d
155154 r = d
156155 else
@@ -160,8 +159,7 @@ subroutine CROTG( a, b, c, s )
160159 u = min ( safmax, max ( safmin, g1 ) )
161160 uu = one / u
162161 gs = g* uu
163- g2 = ABSSQ( gs )
164- d = sqrt ( g2 )
162+ d = abs ( g2 )
165163 s = conjg ( gs ) / d
166164 r = d* u
167165 end if
@@ -176,15 +174,10 @@ subroutine CROTG( a, b, c, s )
176174 f2 = ABSSQ( f )
177175 g2 = ABSSQ( g )
178176 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 )
177+ d = sqrt ( one + ( g2/ f2 ) )
178+ r = f* d
179+ c = one / d
180+ s = conjg ( g )* ( r / h2 )
188181 else
189182!
190183! Use scaled algorithm
@@ -201,27 +194,25 @@ subroutine CROTG( a, b, c, s )
201194 v = min ( safmax, max ( safmin, f1 ) )
202195 vv = one / v
203196 w = v * uu
197+ w2 = w** 2
204198 fs = f* vv
205199 f2 = ABSSQ( fs )
206- h2 = f2* w ** 2 + g2
200+ h2 = f2* w2 + g2
207201 else
208202!
209203! Otherwise use the same scaling for f and g.
210204!
211205 w = one
206+ w2 = one
212207 fs = f* uu
213208 f2 = ABSSQ( fs )
214209 h2 = f2 + g2
215210 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
211+ d = sqrt ( w2 + ( g2/ f2 ) )
212+ c = w / d
213+ r = fs* d
214+ s = conjg ( gs )* ( r / h2 )
215+ r = r* u
225216 end if
226217 end if
227218 a = r
0 commit comments