Skip to content

Fix accuracy issues in the Givens rotations #630

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
wants to merge 3 commits into from
Closed

Fix accuracy issues in the Givens rotations #630

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
bd11fe9
First attempt to fix the precision issue
weslleyspereira Oct 11, 2021
35685b9
Avoids one square operation
weslleyspereira Oct 12, 2021
72e100c
Updates all complex Givens rotations
weslleyspereira Oct 12, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Updates all complex Givens rotations
  • Loading branch information
@weslleyspereira
weslleyspereira committed Oct 12, 2021
commit 72e100c6c5e5c7146d657beb0ead4a1634aae9e3
39 changes: 15 additions & 24 deletions BLAS/SRC/crotg.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ subroutine CROTG( a, b, c, s )
complex(wp) :: a, b, s
! ..
! .. Local Scalars ..
real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
real(wp) :: d, f1, f2, g1, g2, h2, w2, u, uu, v, vv, w
complex(wp) :: f, fs, g, gs, r, t
! ..
! .. Intrinsic Functions ..
Expand All @@ -149,8 +149,7 @@ subroutine CROTG( a, b, c, s )
!
! Use unscaled algorithm
!
g2 = ABSSQ( g )
d = sqrt( g2 )
d = abs( g )
s = conjg( g ) / d
r = d
else
Expand All @@ -160,8 +159,7 @@ subroutine CROTG( a, b, c, s )
u = min( safmax, max( safmin, g1 ) )
uu = one / u
gs = g*uu
g2 = ABSSQ( gs )
d = sqrt( g2 )
d = abs( g2 )
s = conjg( gs ) / d
r = d*u
end if
Expand All @@ -176,15 +174,10 @@ subroutine CROTG( a, b, c, s )
f2 = ABSSQ( f )
g2 = ABSSQ( g )
h2 = f2 + g2
if( f2 > rtmin .and. h2 < rtmax ) then
d = sqrt( f2*h2 )
else
d = sqrt( f2 )*sqrt( h2 )
end if
p = 1 / d
c = f2*p
s = conjg( g )*( f*p )
r = f*( h2*p )
d = sqrt( one + ( g2/f2 ) )
r = f*d
c = one / d
s = conjg( g )*( r / h2 )
else
!
! Use scaled algorithm
Expand All @@ -201,27 +194,25 @@ subroutine CROTG( a, b, c, s )
v = min( safmax, max( safmin, f1 ) )
vv = one / v
w = v * uu
w2 = w**2
fs = f*vv
f2 = ABSSQ( fs )
h2 = f2*w**2 + g2
h2 = f2*w2 + g2
else
!
! Otherwise use the same scaling for f and g.
!
w = one
w2 = one
fs = f*uu
f2 = ABSSQ( fs )
h2 = f2 + g2
end if
if( f2 > rtmin .and. h2 < rtmax ) then
d = sqrt( f2*h2 )
else
d = sqrt( f2 )*sqrt( h2 )
end if
p = 1 / d
c = ( f2*p )*w
s = conjg( gs )*( fs*p )
r = ( fs*( h2*p ) )*u
d = sqrt( w2 + ( g2/f2 ) )
c = w / d
r = fs*d
s = conjg( gs )*( r / h2 )
r = r*u
end if
end if
a = r
Expand Down
39 changes: 15 additions & 24 deletions BLAS/SRC/zrotg.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ subroutine ZROTG( a, b, c, s )
complex(wp) :: a, b, s
! ..
! .. Local Scalars ..
real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
real(wp) :: d, f1, f2, g1, g2, h2, w2, u, uu, v, vv, w
complex(wp) :: f, fs, g, gs, r, t
! ..
! .. Intrinsic Functions ..
Expand All @@ -149,8 +149,7 @@ subroutine ZROTG( a, b, c, s )
!
! Use unscaled algorithm
!
g2 = ABSSQ( g )
d = sqrt( g2 )
d = abs( g )
s = conjg( g ) / d
r = d
else
Expand All @@ -160,8 +159,7 @@ subroutine ZROTG( a, b, c, s )
u = min( safmax, max( safmin, g1 ) )
uu = one / u
gs = g*uu
g2 = ABSSQ( gs )
d = sqrt( g2 )
d = abs( g2 )
s = conjg( gs ) / d
r = d*u
end if
Expand All @@ -176,15 +174,10 @@ subroutine ZROTG( a, b, c, s )
f2 = ABSSQ( f )
g2 = ABSSQ( g )
h2 = f2 + g2
if( f2 > rtmin .and. h2 < rtmax ) then
d = sqrt( f2*h2 )
else
d = sqrt( f2 )*sqrt( h2 )
end if
p = 1 / d
c = f2*p
s = conjg( g )*( f*p )
r = f*( h2*p )
d = sqrt( one + ( g2/f2 ) )
r = f*d
c = one / d
s = conjg( g )*( r / h2 )
else
!
! Use scaled algorithm
Expand All @@ -201,27 +194,25 @@ subroutine ZROTG( a, b, c, s )
v = min( safmax, max( safmin, f1 ) )
vv = one / v
w = v * uu
w2 = w**2
fs = f*vv
f2 = ABSSQ( fs )
h2 = f2*w**2 + g2
h2 = f2*w2 + g2
else
!
! Otherwise use the same scaling for f and g.
!
w = one
w2 = one
fs = f*uu
f2 = ABSSQ( fs )
h2 = f2 + g2
end if
if( f2 > rtmin .and. h2 < rtmax ) then
d = sqrt( f2*h2 )
else
d = sqrt( f2 )*sqrt( h2 )
end if
p = 1 / d
c = ( f2*p )*w
s = conjg( gs )*( fs*p )
r = ( fs*( h2*p ) )*u
d = sqrt( w2 + ( g2/f2 ) )
c = w / d
r = fs*d
s = conjg( gs )*( r / h2 )
r = r*u
end if
end if
a = r
Expand Down
39 changes: 15 additions & 24 deletions SRC/zlartg.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ subroutine ZLARTG( f, g, c, s, r )
complex(wp) f, g, r, s
! ..
! .. Local Scalars ..
real(wp) :: d, f1, f2, g1, g2, h2, p, u, uu, v, vv, w
real(wp) :: d, f1, f2, g1, g2, h2, w2, u, uu, v, vv, w
complex(wp) :: fs, gs, t
! ..
! .. Intrinsic Functions ..
Expand All @@ -154,8 +154,7 @@ subroutine ZLARTG( f, g, c, s, r )
!
! Use unscaled algorithm
!
g2 = ABSSQ( g )
d = sqrt( g2 )
d = abs( g )
s = conjg( g ) / d
r = d
else
Expand All @@ -165,8 +164,7 @@ subroutine ZLARTG( f, g, c, s, r )
u = min( safmax, max( safmin, g1 ) )
uu = one / u
gs = g*uu
g2 = ABSSQ( gs )
d = sqrt( g2 )
d = abs( g2 )
s = conjg( gs ) / d
r = d*u
end if
Expand All @@ -181,15 +179,10 @@ subroutine ZLARTG( f, g, c, s, r )
f2 = ABSSQ( f )
g2 = ABSSQ( g )
h2 = f2 + g2
if( f2 > rtmin .and. h2 < rtmax ) then
d = sqrt( f2*h2 )
else
d = sqrt( f2 )*sqrt( h2 )
end if
p = 1 / d
c = f2*p
s = conjg( g )*( f*p )
r = f*( h2*p )
d = sqrt( one + ( g2/f2 ) )
r = f*d
c = one / d
s = conjg( g )*( r / h2 )
else
!
! Use scaled algorithm
Expand All @@ -206,27 +199,25 @@ subroutine ZLARTG( f, g, c, s, r )
v = min( safmax, max( safmin, f1 ) )
vv = one / v
w = v * uu
w2 = w**2
fs = f*vv
f2 = ABSSQ( fs )
h2 = f2*w**2 + g2
h2 = f2*w2 + g2
else
!
! Otherwise use the same scaling for f and g.
!
w = one
w2 = one
fs = f*uu
f2 = ABSSQ( fs )
h2 = f2 + g2
end if
if( f2 > rtmin .and. h2 < rtmax ) then
d = sqrt( f2*h2 )
else
d = sqrt( f2 )*sqrt( h2 )
end if
p = 1 / d
c = ( f2*p )*w
s = conjg( gs )*( fs*p )
r = ( fs*( h2*p ) )*u
d = sqrt( w2 + ( g2/f2 ) )
c = w / d
r = fs*d
s = conjg( gs )*( r / h2 )
r = r*u
end if
end if
return
Expand Down