phynexis-deps/lapack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #1245 from nakatamaho/fix/lartg-stable-complex-s

Browse Source 
avoid large intermediates in complex LARTG
This commit is contained in:
langou
2026-07-15 03:22:29 -06:00
committed by GitHub
parent 31a1d7844b cb988d7863
commit 379b30d311
6 changed files with 20 additions and 24 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SRC/cbbcsd.f
+3 -4
View File
@@ -351,9 +351,9 @@
* .. Parameters ..
INTEGER MAXITR
PARAMETER ( MAXITR = 6 )
REAL HUNDRED, MEIGHTH, ONE, TEN, ZERO
REAL HUNDRED, MEIGHTH, ZERO, ONE, TEN
PARAMETER ( HUNDRED = 100.0E0, MEIGHTH = -0.125E0,
$ ONE = 1.0E0, TEN = 10.0E0, ZERO = 0.0E0 )
$ ZERO = 0.0E0, ONE = 1.0E0, TEN = 10.0E0 )
COMPLEX NEGONECOMPLEX
PARAMETER ( NEGONECOMPLEX = (-1.0E0,0.0E0) )
REAL PIOVER2
@@ -576,7 +576,7 @@
END IF
ELSE
NU = SIGMA21
MU = SQRT( 1.0 - NU**2 )
MU = SQRT( ONE - NU**2 )
IF( NU .LT. THRESH ) THEN
MU = ONE
NU = ZERO
@@ -1114,4 +1114,3 @@
* End of CBBCSD
*
END
SRC/clartg.f90
+4 -4
View File
@@ -204,7 +204,7 @@ subroutine CLARTG( f, g, c, s, r )
rtmax = rtmax * 2
if( f2 > rtmin .and. h2 < rtmax ) then
! safmin <= sqrt( f2*h2 ) <= safmax
s = conjg( g ) * ( f / sqrt( f2*h2 ) )
s = ( f / sqrt( f2 ) ) * ( conjg( g ) / sqrt( h2 ) )
else
s = conjg( g ) * ( r / h2 )
end if
@@ -224,7 +224,7 @@ subroutine CLARTG( f, g, c, s, r )
! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
r = f * ( h2 / d )
end if
s = conjg( g ) * ( f / d )
s = ( f / sqrt( f2 ) ) * ( conjg( g ) / sqrt( h2 ) )
end if
else
!
@@ -260,7 +260,7 @@ subroutine CLARTG( f, g, c, s, r )
rtmax = rtmax * 2
if( f2 > rtmin .and. h2 < rtmax ) then
! safmin <= sqrt( f2*h2 ) <= safmax
s = conjg( gs ) * ( fs / sqrt( f2*h2 ) )
s = ( fs / sqrt( f2 ) ) * ( conjg( gs ) / sqrt( h2 ) )
else
s = conjg( gs ) * ( r / h2 )
end if
@@ -280,7 +280,7 @@ subroutine CLARTG( f, g, c, s, r )
! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
r = fs * ( h2 / d )
end if
s = conjg( gs ) * ( fs / d )
s = ( fs / sqrt( f2 ) ) * ( conjg( gs ) / sqrt( h2 ) )
end if
! Rescale c and r
c = c * w
SRC/dbbcsd.f
+3 -4
View File
@@ -351,9 +351,9 @@
* .. Parameters ..
INTEGER MAXITR
PARAMETER ( MAXITR = 6 )
DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, TEN, ZERO
DOUBLE PRECISION HUNDRED, MEIGHTH, ZERO, ONE, TEN
PARAMETER ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0,
$ ONE = 1.0D0, TEN = 10.0D0, ZERO = 0.0D0 )
$ ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
DOUBLE PRECISION NEGONE
PARAMETER ( NEGONE = -1.0D0 )
DOUBLE PRECISION PIOVER2
@@ -576,7 +576,7 @@
END IF
ELSE
NU = SIGMA21
MU = SQRT( 1.0 - NU**2 )
MU = SQRT( ONE - NU**2 )
IF( NU .LT. THRESH ) THEN
MU = ONE
NU = ZERO
@@ -1108,4 +1108,3 @@
* End of DBBCSD
*
END
SRC/sbbcsd.f
+3 -4
View File
@@ -351,9 +351,9 @@
* .. Parameters ..
INTEGER MAXITR
PARAMETER ( MAXITR = 6 )
REAL HUNDRED, MEIGHTH, ONE, TEN, ZERO
REAL HUNDRED, MEIGHTH, ZERO, ONE, TEN
PARAMETER ( HUNDRED = 100.0E0, MEIGHTH = -0.125E0,
$ ONE = 1.0E0, TEN = 10.0E0, ZERO = 0.0E0 )
$ ZERO = 0.0E0, ONE = 1.0E0, TEN = 10.0E0 )
REAL NEGONE
PARAMETER ( NEGONE = -1.0E0 )
REAL PIOVER2
@@ -576,7 +576,7 @@
END IF
ELSE
NU = SIGMA21
MU = SQRT( 1.0 - NU**2 )
MU = SQRT( ONE - NU**2 )
IF( NU .LT. THRESH ) THEN
MU = ONE
NU = ZERO
@@ -1108,4 +1108,3 @@
* End of SBBCSD
*
END
SRC/zbbcsd.f
+3 -4
View File
@@ -351,9 +351,9 @@
* .. Parameters ..
INTEGER MAXITR
PARAMETER ( MAXITR = 6 )
DOUBLE PRECISION HUNDRED, MEIGHTH, ONE, TEN, ZERO
DOUBLE PRECISION HUNDRED, MEIGHTH, ZERO, ONE, TEN
PARAMETER ( HUNDRED = 100.0D0, MEIGHTH = -0.125D0,
$ ONE = 1.0D0, TEN = 10.0D0, ZERO = 0.0D0 )
$ ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
COMPLEX*16 NEGONECOMPLEX
PARAMETER ( NEGONECOMPLEX = (-1.0D0,0.0D0) )
DOUBLE PRECISION PIOVER2
@@ -575,7 +575,7 @@
END IF
ELSE
NU = SIGMA21
MU = SQRT( 1.0 - NU**2 )
MU = SQRT( ONE - NU**2 )
IF( NU .LT. THRESH ) THEN
MU = ONE
NU = ZERO
@@ -1113,4 +1113,3 @@
* End of ZBBCSD
*
END
SRC/zlartg.f90
+4 -4
View File
@@ -204,7 +204,7 @@ subroutine ZLARTG( f, g, c, s, r )
rtmax = rtmax * 2
if( f2 > rtmin .and. h2 < rtmax ) then
! safmin <= sqrt( f2*h2 ) <= safmax
s = conjg( g ) * ( f / sqrt( f2*h2 ) )
s = ( f / sqrt( f2 ) ) * ( conjg( g ) / sqrt( h2 ) )
else
s = conjg( g ) * ( r / h2 )
end if
@@ -224,7 +224,7 @@ subroutine ZLARTG( f, g, c, s, r )
! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
r = f * ( h2 / d )
end if
s = conjg( g ) * ( f / d )
s = ( f / sqrt( f2 ) ) * ( conjg( g ) / sqrt( h2 ) )
end if
else
!
@@ -260,7 +260,7 @@ subroutine ZLARTG( f, g, c, s, r )
rtmax = rtmax * 2
if( f2 > rtmin .and. h2 < rtmax ) then
! safmin <= sqrt( f2*h2 ) <= safmax
s = conjg( gs ) * ( fs / sqrt( f2*h2 ) )
s = ( fs / sqrt( f2 ) ) * ( conjg( gs ) / sqrt( h2 ) )
else
s = conjg( gs ) * ( r / h2 )
end if
@@ -280,7 +280,7 @@ subroutine ZLARTG( f, g, c, s, r )
! sqrt(safmin) <= f2 * sqrt(safmax) <= h2 / sqrt(f2 * h2) <= h2 * (safmin / f2) <= h2 <= safmax
r = fs * ( h2 / d )
end if
s = conjg( gs ) * ( fs / d )
s = ( fs / sqrt( f2 ) ) * ( conjg( gs ) / sqrt( h2 ) )
end if
! Rescale c and r
c = c * w