Skip to content

Commit 37a1a1e

Browse files
weslleyspereiraweslleyspereira
committed
Fix all other Givens rotation routines
1 parent ac11f62 commit 37a1a1e

File tree

6 files changed

+76
-124
lines changed

6 files changed

+76
-124
lines changed

‎BLAS/SRC/crotg.f90‎

Lines changed: 15 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -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, u, v, w
126126
complex(wp) :: f, fs, g, gs, r, t
127127
! ..
128128
! .. Intrinsic Functions ..
@@ -149,19 +149,16 @@ 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
157156
!
158157
! Use scaled algorithm
159158
!
160159
u = min( safmax, max( safmin, g1 ) )
161-
uu = one / u
162-
gs = g*uu
163-
g2 = ABSSQ( gs )
164-
d = sqrt( g2 )
160+
gs = g / u
161+
d = abs( gs )
165162
s = conjg( gs ) / d
166163
r = d*u
167164
end if
@@ -181,35 +178,32 @@ subroutine CROTG( a, b, c, s )
181178
else
182179
d = sqrt( f2 )*sqrt( h2 )
183180
end if
184-
p = 1 / d
185-
c = f2*p
186-
s = conjg( g )*( f*p )
187-
r = f*( h2*p )
181+
c = f2 / d
182+
s = conjg( g )*( f / d )
183+
r = f*( h2 / d )
188184
else
189185
!
190186
! Use scaled algorithm
191187
!
192188
u = min( safmax, max( safmin, f1, g1 ) )
193-
uu = one / u
194-
gs = g*uu
189+
gs = g / u
195190
g2 = ABSSQ( gs )
196-
if( f1*uu < rtmin ) then
191+
if( f1 / u < rtmin ) then
197192
!
198193
! f is not well-scaled when scaled by g1.
199194
! Use a different scaling for f.
200195
!
201196
v = min( safmax, max( safmin, f1 ) )
202-
vv = one / v
203-
w = v * uu
204-
fs = f*vv
197+
w = v / u
198+
fs = f / v
205199
f2 = ABSSQ( fs )
206200
h2 = f2*w**2 + g2
207201
else
208202
!
209203
! Otherwise use the same scaling for f and g.
210204
!
211205
w = one
212-
fs = f*uu
206+
fs = f / u
213207
f2 = ABSSQ( fs )
214208
h2 = f2 + g2
215209
end if
@@ -218,10 +212,9 @@ subroutine CROTG( a, b, c, s )
218212
else
219213
d = sqrt( f2 )*sqrt( h2 )
220214
end if
221-
p = 1 / d
222-
c = ( f2*p )*w
223-
s = conjg( gs )*( fs*p )
224-
r = ( fs*( h2*p ) )*u
215+
c = ( f2 / d )*w
216+
s = conjg( gs )*( fs / d )
217+
r = ( fs*( h2 / d ) )*u
225218
end if
226219
end if
227220
a = r

‎BLAS/SRC/zrotg.f90‎

Lines changed: 15 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,7 @@ subroutine ZROTG( 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, u, v, w
126126
complex(wp) :: f, fs, g, gs, r, t
127127
! ..
128128
! .. Intrinsic Functions ..
@@ -149,19 +149,16 @@ subroutine ZROTG( 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
157156
!
158157
! Use scaled algorithm
159158
!
160159
u = min( safmax, max( safmin, g1 ) )
161-
uu = one / u
162-
gs = g*uu
163-
g2 = ABSSQ( gs )
164-
d = sqrt( g2 )
160+
gs = g / u
161+
d = abs( gs )
165162
s = conjg( gs ) / d
166163
r = d*u
167164
end if
@@ -181,35 +178,32 @@ subroutine ZROTG( a, b, c, s )
181178
else
182179
d = sqrt( f2 )*sqrt( h2 )
183180
end if
184-
p = 1 / d
185-
c = f2*p
186-
s = conjg( g )*( f*p )
187-
r = f*( h2*p )
181+
c = f2 / d
182+
s = conjg( g )*( f / d )
183+
r = f*( h2 / d )
188184
else
189185
!
190186
! Use scaled algorithm
191187
!
192188
u = min( safmax, max( safmin, f1, g1 ) )
193-
uu = one / u
194-
gs = g*uu
189+
gs = g / u
195190
g2 = ABSSQ( gs )
196-
if( f1*uu < rtmin ) then
191+
if( f1 / u < rtmin ) then
197192
!
198193
! f is not well-scaled when scaled by g1.
199194
! Use a different scaling for f.
200195
!
201196
v = min( safmax, max( safmin, f1 ) )
202-
vv = one / v
203-
w = v * uu
204-
fs = f*vv
197+
w = v / u
198+
fs = f / v
205199
f2 = ABSSQ( fs )
206200
h2 = f2*w**2 + g2
207201
else
208202
!
209203
! Otherwise use the same scaling for f and g.
210204
!
211205
w = one
212-
fs = f*uu
206+
fs = f / u
213207
f2 = ABSSQ( fs )
214208
h2 = f2 + g2
215209
end if
@@ -218,10 +212,9 @@ subroutine ZROTG( a, b, c, s )
218212
else
219213
d = sqrt( f2 )*sqrt( h2 )
220214
end if
221-
p = 1 / d
222-
c = ( f2*p )*w
223-
s = conjg( gs )*( fs*p )
224-
r = ( fs*( h2*p ) )*u
215+
c = ( f2 / d )*w
216+
s = conjg( gs )*( fs / d )
217+
r = ( fs*( h2 / d ) )*u
225218
end if
226219
end if
227220
a = r

‎SRC/clartg.f90‎

Lines changed: 13 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -178,33 +178,22 @@ subroutine CLARTG( f, g, c, s, r )
178178
f2 = ABSSQ( f )
179179
g2 = ABSSQ( g )
180180
h2 = f2 + g2
181-
if( f2 > safmin * g2 ) then
182-
d = sqrt( one + g2/f2 )
183-
c = one / d
184-
if( f2 > rtmin .and. h2 < rtmax ) then
185-
s = conjg( g )*( f / sqrt( f2*h2 ) )
186-
else
187-
s = conjg( g )*( f /( f2*d ) )
188-
end if
189-
r = f * d
181+
if( f2 > rtmin .and. h2 < rtmax ) then
182+
d = sqrt( f2*h2 )
190183
else
191-
if( f2 > rtmin .and. h2 < rtmax ) then
192-
d = sqrt( f2*h2 )
193-
else
194-
d = sqrt( f2 )*sqrt( h2 )
195-
end if
196-
c = f2 / d
197-
s = conjg( g )*( f / d )
198-
r = f*( h2 / d )
184+
d = sqrt( f2 )*sqrt( h2 )
199185
end if
186+
c = f2 / d
187+
s = conjg( g )*( f / d )
188+
r = f*( h2 / d )
200189
else
201190
!
202191
! Use scaled algorithm
203192
!
204193
u = min( safmax, max( safmin, f1, g1 ) )
205194
gs = g / u
206195
g2 = ABSSQ( gs )
207-
if( f1 < rtmin * u ) then
196+
if( f1 / u < rtmin ) then
208197
!
209198
! f is not well-scaled when scaled by g1.
210199
! Use a different scaling for f.
@@ -223,26 +212,14 @@ subroutine CLARTG( f, g, c, s, r )
223212
f2 = ABSSQ( fs )
224213
h2 = f2 + g2
225214
end if
226-
if( f2 > safmin * g2 ) then
227-
! Use a precise algorithm
228-
d = sqrt( w**2 + g2/f2 )
229-
c = w / d
230-
if( f2 > rtmin .and. h2 < rtmax ) then
231-
s = conjg( gs )*( fs / sqrt( f2*h2 ) )
232-
else
233-
s = conjg( gs )*( fs / ( f2*d ) )
234-
end if
235-
r = ( fs * d ) * u
215+
if( f2 > rtmin .and. h2 < rtmax ) then
216+
d = sqrt( f2*h2 )
236217
else
237-
if( f2 > rtmin .and. h2 < rtmax ) then
238-
d = sqrt( f2*h2 )
239-
else
240-
d = sqrt( f2 )*sqrt( h2 )
241-
end if
242-
c = ( f2 / d )*w
243-
s = conjg( gs )*( fs / d )
244-
r = ( fs*( h2 / d ) )*u
218+
d = sqrt( f2 )*sqrt( h2 )
245219
end if
220+
c = ( f2 / d )*w
221+
s = conjg( gs )*( fs / d )
222+
r = ( fs*( h2 / d ) )*u
246223
end if
247224
end if
248225
return

‎SRC/dlartg.f90‎

Lines changed: 9 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -123,7 +123,7 @@ subroutine DLARTG( f, g, c, s, r )
123123
real(wp) :: c, f, g, r, s
124124
! ..
125125
! .. Local Scalars ..
126-
real(wp) :: d, f1, fs, g1, gs, p, u, uu
126+
real(wp) :: d, f1, fs, g1, gs, u
127127
! ..
128128
! .. Intrinsic Functions ..
129129
intrinsic :: abs, sign, sqrt
@@ -143,20 +143,18 @@ subroutine DLARTG( f, g, c, s, r )
143143
else if( f1 > rtmin .and. f1 < rtmax .and. &
144144
g1 > rtmin .and. g1 < rtmax ) then
145145
d = sqrt( f*f + g*g )
146-
p = one / d
147-
c = f1*p
148-
s = g*sign( p, f )
146+
c = f1 / d
149147
r = sign( d, f )
148+
s = g / r
150149
else
151150
u = min( safmax, max( safmin, f1, g1 ) )
152-
uu = one / u
153-
fs = f*uu
154-
gs = g*uu
151+
fs = f / u
152+
gs = g / u
155153
d = sqrt( fs*fs + gs*gs )
156-
p = one / d
157-
c = abs( fs )*p
158-
s = gs*sign( p, f )
159-
r = sign( d, f )*u
154+
c = abs( fs ) / d
155+
r = sign( d, f )
156+
s = gs / r
157+
r = r * u
160158
end if
161159
return
162160
end subroutine

‎SRC/slartg.f90‎

Lines changed: 9 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -123,7 +123,7 @@ subroutine SLARTG( f, g, c, s, r )
123123
real(wp) :: c, f, g, r, s
124124
! ..
125125
! .. Local Scalars ..
126-
real(wp) :: d, f1, fs, g1, gs, p, u, uu
126+
real(wp) :: d, f1, fs, g1, gs, u
127127
! ..
128128
! .. Intrinsic Functions ..
129129
intrinsic :: abs, sign, sqrt
@@ -143,20 +143,18 @@ subroutine SLARTG( f, g, c, s, r )
143143
else if( f1 > rtmin .and. f1 < rtmax .and. &
144144
g1 > rtmin .and. g1 < rtmax ) then
145145
d = sqrt( f*f + g*g )
146-
p = one / d
147-
c = f1*p
148-
s = g*sign( p, f )
146+
c = f1 / d
149147
r = sign( d, f )
148+
s = g / r
150149
else
151150
u = min( safmax, max( safmin, f1, g1 ) )
152-
uu = one / u
153-
fs = f*uu
154-
gs = g*uu
151+
fs = f / u
152+
gs = g / u
155153
d = sqrt( fs*fs + gs*gs )
156-
p = one / d
157-
c = abs( fs )*p
158-
s = gs*sign( p, f )
159-
r = sign( d, f )*u
154+
c = abs( fs ) / d
155+
r = sign( d, f )
156+
s = gs / r
157+
r = r*u
160158
end if
161159
return
162160
end subroutine

0 commit comments

Comments
 (0)