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

Merge pull request #1191 from angsch/issue-1166

Browse Source 
Update step length selection in LAED4 overshoot fallback
This commit is contained in:
langou
2026-03-05 22:36:31 -07:00
committed by GitHub
parent c0c64400da b6b6288080
commit 00cf62bc1b
2 changed files with 67 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SRC/dlaed4.f
+34 -7
View File
@@ -168,8 +168,9 @@
LOGICAL ORGATI, SWTCH, SWTCH3
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
$ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
$ RHOINV, TAU, TEMP, TEMP1, W
$ EPS, ERRETM, ETA, ETA1, ETA2, MIDPT, PHI,
$ PREW, PSI, RHOINV, TAU, TEMP, TEMP1, W
* ..
* .. Local Arrays ..
DOUBLE PRECISION ZZ( 3 )
@@ -182,7 +183,7 @@
EXTERNAL DLAED5, DLAED6
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
INTRINSIC ABS, MAX, MIN, SIGN, SQRT
* ..
* .. Executable Statements ..
*
@@ -349,10 +350,23 @@
$ ETA = -W / ( DPSI+DPHI )
TEMP = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
ETA1 = -W / ( DPSI+DPHI )
TEMP = TAU + ETA1
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ETA2 = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
ETA2 = ( DLTLB-TAU ) / TWO
END IF
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
*
* If Newton step is within boundaries,
* use the geometric mean of the distance
* and keep the direction (sign).
*
ETA = SIGN(ONE, ETA1) *
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
ELSE
ETA = ETA2
END IF
END IF
DO 50 J = 1, N
@@ -848,10 +862,23 @@
$ ETA = -W / DW
TEMP = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
ETA1 = -W / DW
TEMP = TAU + ETA1
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ETA2 = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
ETA2 = ( DLTLB-TAU ) / TWO
END IF
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
*
* If Newton step is within boundaries,
* use the geometric mean of the distance
* and keep the direction (sign).
*
ETA = SIGN( ONE, ETA1 ) *
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
ELSE
ETA = ETA2
END IF
END IF
*
SRC/slaed4.f
+33 -7
View File
@@ -168,8 +168,8 @@
LOGICAL ORGATI, SWTCH, SWTCH3
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
REAL A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
$ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
$ RHOINV, TAU, TEMP, TEMP1, W
$ EPS, ERRETM, ETA, ETA1, ETA2, MIDPT, PHI,
$ PREW, PSI, RHOINV, TAU, TEMP, TEMP1, W
* ..
* .. Local Arrays ..
REAL ZZ( 3 )
@@ -182,7 +182,7 @@
EXTERNAL SLAED5, SLAED6
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, MAX, MIN, SQRT
INTRINSIC ABS, MAX, MIN, SIGN, SQRT
* ..
* .. Executable Statements ..
*
@@ -349,10 +349,23 @@
$ ETA = -W / ( DPSI+DPHI )
TEMP = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
ETA1 = -W / ( DPSI+DPHI )
TEMP = TAU + ETA1
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ETA2 = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
ETA2 = ( DLTLB-TAU ) / TWO
END IF
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
*
* If Newton step is within boundaries,
* use the geometric mean of the distance
* and keep the direction (sign).
*
ETA = SIGN(ONE, ETA1) *
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
ELSE
ETA = ETA2
END IF
END IF
DO 50 J = 1, N
@@ -848,10 +861,23 @@
$ ETA = -W / DW
TEMP = TAU + ETA
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
ETA1 = -W / DW
TEMP = TAU + ETA1
IF( W.LT.ZERO ) THEN
ETA = ( DLTUB-TAU ) / TWO
ETA2 = ( DLTUB-TAU ) / TWO
ELSE
ETA = ( DLTLB-TAU ) / TWO
ETA2 = ( DLTLB-TAU ) / TWO
END IF
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
*
* If Newton step is within boundaries,
* use the geometric mean of the distance
* and keep the direction (sign).
*
ETA = SIGN( ONE, ETA1 ) *
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
ELSE
ETA = ETA2
END IF
END IF
*