Skip to content

Commit 2495f1c

Browse files
Solves a precision bug in clartg
1 parent 79aa0f2 commit 2495f1c

File tree

1 file changed

+10
-1
lines changed

1 file changed

+10
-1
lines changed

SRC/clartg.f90

Lines changed: 10 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -187,7 +187,11 @@ subroutine CLARTG( f, g, c, s, r )
187187
d = sqrt( f2 )*sqrt( h2 )
188188
end if
189189
p = 1 / d
190-
c = f2*p
190+
if( f2 > safmin * g2 ) then
191+
c = 1 / sqrt( one + g2/f2 )
192+
else
193+
c = f2*p
194+
end if
191195
s = conjg( g )*( f*p )
192196
r = f*( h2*p )
193197
else
@@ -224,6 +228,11 @@ subroutine CLARTG( f, g, c, s, r )
224228
d = sqrt( f2 )*sqrt( h2 )
225229
end if
226230
p = 1 / d
231+
if( f2 > safmin * g2 ) then
232+
c = (1 / sqrt( one + g2/f2 )) * w
233+
else
234+
c = ( f2*p )*w
235+
end if
227236
c = ( f2*p )*w
228237
s = conjg( gs )*( fs*p )
229238
r = ( fs*( h2*p ) )*u

0 commit comments

Comments
 (0)