Skip to content

Commit 81d1029

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #5694 from martin-frbg/lapack1191
Update step length selection in ?LAED4 fallback (Reference-LAPACK PR 1191)
2 parents 4956446 + aa6a59a commit 81d1029

2 files changed

Lines changed: 75 additions & 22 deletions

File tree

‎lapack-netlib/SRC/dlaed4.f‎

Lines changed: 38 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download DLAED4 + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed4.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed4.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed4.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -132,7 +130,7 @@
132130
*> \author Univ. of Colorado Denver
133131
*> \author NAG Ltd.
134132
*
135-
*> \ingroup auxOTHERcomputational
133+
*> \ingroup laed4
136134
*
137135
*> \par Contributors:
138136
* ==================
@@ -142,6 +140,7 @@
142140
*>
143141
* =====================================================================
144142
SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
143+
IMPLICIT NONE
145144
*
146145
* -- LAPACK computational routine --
147146
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -169,8 +168,9 @@ SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
169168
LOGICAL ORGATI, SWTCH, SWTCH3
170169
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
171170
DOUBLE PRECISION A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
172-
$ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
173-
$ RHOINV, TAU, TEMP, TEMP1, W
171+
$ EPS, ERRETM, ETA, ETA1, ETA2, MIDPT, PHI,
172+
$ PREW, PSI, RHOINV, TAU, TEMP, TEMP1, W
173+
174174
* ..
175175
* .. Local Arrays ..
176176
DOUBLE PRECISION ZZ( 3 )
@@ -183,7 +183,7 @@ SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
183183
EXTERNAL DLAED5, DLAED6
184184
* ..
185185
* .. Intrinsic Functions ..
186-
INTRINSIC ABS, MAX, MIN, SQRT
186+
INTRINSIC ABS, MAX, MIN, SIGN, SQRT
187187
* ..
188188
* .. Executable Statements ..
189189
*
@@ -350,10 +350,23 @@ SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
350350
$ ETA = -W / ( DPSI+DPHI )
351351
TEMP = TAU + ETA
352352
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
353+
ETA1 = -W / ( DPSI+DPHI )
354+
TEMP = TAU + ETA1
353355
IF( W.LT.ZERO ) THEN
354-
ETA = ( DLTUB-TAU ) / TWO
356+
ETA2 = ( DLTUB-TAU ) / TWO
355357
ELSE
356-
ETA = ( DLTLB-TAU ) / TWO
358+
ETA2 = ( DLTLB-TAU ) / TWO
359+
END IF
360+
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
361+
*
362+
* If Newton step is within boundaries,
363+
* use the geometric mean of the distance
364+
* and keep the direction (sign).
365+
*
366+
ETA = SIGN(ONE, ETA1) *
367+
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
368+
ELSE
369+
ETA = ETA2
357370
END IF
358371
END IF
359372
DO 50 J = 1, N
@@ -832,7 +845,8 @@ SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
832845
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
833846
END IF
834847
END IF
835-
CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA,
848+
CALL DLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W,
849+
$ ETA,
836850
$ INFO )
837851
IF( INFO.NE.0 )
838852
$ GO TO 250
@@ -848,10 +862,23 @@ SUBROUTINE DLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
848862
$ ETA = -W / DW
849863
TEMP = TAU + ETA
850864
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
865+
ETA1 = -W / DW
866+
TEMP = TAU + ETA1
851867
IF( W.LT.ZERO ) THEN
852-
ETA = ( DLTUB-TAU ) / TWO
868+
ETA2 = ( DLTUB-TAU ) / TWO
853869
ELSE
854-
ETA = ( DLTLB-TAU ) / TWO
870+
ETA2 = ( DLTLB-TAU ) / TWO
871+
END IF
872+
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
873+
*
874+
* If Newton step is within boundaries,
875+
* use the geometric mean of the distance
876+
* and keep the direction (sign).
877+
*
878+
ETA = SIGN( ONE, ETA1 ) *
879+
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
880+
ELSE
881+
ETA = ETA2
855882
END IF
856883
END IF
857884
*

‎lapack-netlib/SRC/slaed4.f‎

Lines changed: 37 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -5,15 +5,13 @@
55
* Online html documentation available at
66
* http://www.netlib.org/lapack/explore-html/
77
*
8-
*> \htmlonly
98
*> Download SLAED4 + dependencies
109
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed4.f">
1110
*> [TGZ]</a>
1211
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed4.f">
1312
*> [ZIP]</a>
1413
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed4.f">
1514
*> [TXT]</a>
16-
*> \endhtmlonly
1715
*
1816
* Definition:
1917
* ===========
@@ -132,7 +130,7 @@
132130
*> \author Univ. of Colorado Denver
133131
*> \author NAG Ltd.
134132
*
135-
*> \ingroup auxOTHERcomputational
133+
*> \ingroup laed4
136134
*
137135
*> \par Contributors:
138136
* ==================
@@ -142,6 +140,7 @@
142140
*>
143141
* =====================================================================
144142
SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
143+
IMPLICIT NONE
145144
*
146145
* -- LAPACK computational routine --
147146
* -- LAPACK is a software package provided by Univ. of Tennessee, --
@@ -169,8 +168,8 @@ SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
169168
LOGICAL ORGATI, SWTCH, SWTCH3
170169
INTEGER II, IIM1, IIP1, IP1, ITER, J, NITER
171170
REAL A, B, C, DEL, DLTLB, DLTUB, DPHI, DPSI, DW,
172-
$ EPS, ERRETM, ETA, MIDPT, PHI, PREW, PSI,
173-
$ RHOINV, TAU, TEMP, TEMP1, W
171+
$ EPS, ERRETM, ETA, ETA1, ETA2, MIDPT, PHI,
172+
$ PREW, PSI, RHOINV, TAU, TEMP, TEMP1, W
174173
* ..
175174
* .. Local Arrays ..
176175
REAL ZZ( 3 )
@@ -183,7 +182,7 @@ SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
183182
EXTERNAL SLAED5, SLAED6
184183
* ..
185184
* .. Intrinsic Functions ..
186-
INTRINSIC ABS, MAX, MIN, SQRT
185+
INTRINSIC ABS, MAX, MIN, SIGN, SQRT
187186
* ..
188187
* .. Executable Statements ..
189188
*
@@ -350,10 +349,23 @@ SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
350349
$ ETA = -W / ( DPSI+DPHI )
351350
TEMP = TAU + ETA
352351
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
352+
ETA1 = -W / ( DPSI+DPHI )
353+
TEMP = TAU + ETA1
353354
IF( W.LT.ZERO ) THEN
354-
ETA = ( DLTUB-TAU ) / TWO
355+
ETA2 = ( DLTUB-TAU ) / TWO
355356
ELSE
356-
ETA = ( DLTLB-TAU ) / TWO
357+
ETA2 = ( DLTLB-TAU ) / TWO
358+
END IF
359+
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
360+
*
361+
* If Newton step is within boundaries,
362+
* use the geometric mean of the distance
363+
* and keep the direction (sign).
364+
*
365+
ETA = SIGN(ONE, ETA1) *
366+
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
367+
ELSE
368+
ETA = ETA2
357369
END IF
358370
END IF
359371
DO 50 J = 1, N
@@ -832,7 +844,8 @@ SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
832844
ZZ( 3 ) = Z( IIP1 )*Z( IIP1 )
833845
END IF
834846
END IF
835-
CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W, ETA,
847+
CALL SLAED6( NITER, ORGATI, C, DELTA( IIM1 ), ZZ, W,
848+
$ ETA,
836849
$ INFO )
837850
IF( INFO.NE.0 )
838851
$ GO TO 250
@@ -848,10 +861,23 @@ SUBROUTINE SLAED4( N, I, D, Z, DELTA, RHO, DLAM, INFO )
848861
$ ETA = -W / DW
849862
TEMP = TAU + ETA
850863
IF( TEMP.GT.DLTUB .OR. TEMP.LT.DLTLB ) THEN
864+
ETA1 = -W / DW
865+
TEMP = TAU + ETA1
851866
IF( W.LT.ZERO ) THEN
852-
ETA = ( DLTUB-TAU ) / TWO
867+
ETA2 = ( DLTUB-TAU ) / TWO
853868
ELSE
854-
ETA = ( DLTLB-TAU ) / TWO
869+
ETA2 = ( DLTLB-TAU ) / TWO
870+
END IF
871+
IF ( DLTLB.LE.TEMP .AND. TEMP.LE.DLTUB ) THEN
872+
*
873+
* If Newton step is within boundaries,
874+
* use the geometric mean of the distance
875+
* and keep the direction (sign).
876+
*
877+
ETA = SIGN( ONE, ETA1 ) *
878+
$ SQRT( ABS( ETA1 ) ) * SQRT( ABS( ETA2 ) )
879+
ELSE
880+
ETA = ETA2
855881
END IF
856882
END IF
857883
*

0 commit comments

Comments
 (0)