diff options
| author | julie <julielangou@users.noreply.github.com> | 2012-04-18 16:35:07 +0000 |
|---|---|---|
| committer | julie <julielangou@users.noreply.github.com> | 2012-04-18 16:35:07 +0000 |
| commit | 9769575ceb533886d0bfffc22d47ce55682727b5 (patch) | |
| tree | 66595c865b74ec4f835fc59813eba694b0510e5e | |
| parent | 90489aa716f68f31bfac3db4e62c9b41963d098b (diff) | |
| download | lapack-9769575ceb533886d0bfffc22d47ce55682727b5.tar.gz lapack-9769575ceb533886d0bfffc22d47ce55682727b5.tar.bz2 lapack-9769575ceb533886d0bfffc22d47ce55682727b5.zip | |
Correct bug sent by NAG people
Got a new bug for you! Mick Pont found this problem in DLAED6. The code in question is the 40 loop - when DSCALE(I)=TAU you get a divide by zero (rare in practice). This can cause some compilers to immediately stop, e.g. the Sun compiler.
Mick proposed solution is below:
DO 40 I = 1, 3
IF (DSCALE( I ).NE.TAU) THEN
TEMP = ONE / ( DSCALE( I )-TAU )
TEMP1 = ZSCALE( I )*TEMP
TEMP2 = TEMP1*TEMP
TEMP3 = TEMP2*TEMP
TEMP4 = TEMP1 / DSCALE( I )
FC = FC + TEMP4
ERRETM = ERRETM + ABS( TEMP4 )
DF = DF + TEMP2
DDF = DDF + TEMP3
ELSE
* On rare occasions dscale(i) can be exactly equal to
* tau, leading to division by zero. If no trap occurs,
* there is no problem; the quantities above all overflow
* and the test on abs(f) below sends you to the end
* with good results. If a trap occurs, though, the
* user program will stop. Avoid that happening by
* jumping directly out.
GO TO 60
END IF
40 CONTINUE
This seems to work OK in our testing.
| -rw-r--r-- | SRC/dlaed6.f | 28 | ||||
| -rw-r--r-- | SRC/slaed6.f | 28 |
2 files changed, 32 insertions, 24 deletions
diff --git a/SRC/dlaed6.f b/SRC/dlaed6.f index 1b55c1a3..e8e9709e 100644 --- a/SRC/dlaed6.f +++ b/SRC/dlaed6.f @@ -115,7 +115,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date November 2011 +*> \date April 2012 * *> \ingroup auxOTHERcomputational * @@ -140,10 +140,10 @@ * ===================================================================== SUBROUTINE DLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO ) * -* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK computational routine (version 3.4.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2011 +* April 2012 * * .. Scalar Arguments .. LOGICAL ORGATI @@ -370,15 +370,19 @@ DF = ZERO DDF = ZERO DO 40 I = 1, 3 - TEMP = ONE / ( DSCALE( I )-TAU ) - TEMP1 = ZSCALE( I )*TEMP - TEMP2 = TEMP1*TEMP - TEMP3 = TEMP2*TEMP - TEMP4 = TEMP1 / DSCALE( I ) - FC = FC + TEMP4 - ERRETM = ERRETM + ABS( TEMP4 ) - DF = DF + TEMP2 - DDF = DDF + TEMP3 + IF (DSCALE( I ).NE.TAU) THEN + TEMP = ONE / ( DSCALE( I )-TAU ) + TEMP1 = ZSCALE( I )*TEMP + TEMP2 = TEMP1*TEMP + TEMP3 = TEMP2*TEMP + TEMP4 = TEMP1 / DSCALE( I ) + FC = FC + TEMP4 + ERRETM = ERRETM + ABS( TEMP4 ) + DF = DF + TEMP2 + DDF = DDF + TEMP3 + ELSE + GO TO 60 + END IF 40 CONTINUE F = FINIT + TAU*FC ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + diff --git a/SRC/slaed6.f b/SRC/slaed6.f index 515d5f04..94194199 100644 --- a/SRC/slaed6.f +++ b/SRC/slaed6.f @@ -115,7 +115,7 @@ *> \author Univ. of Colorado Denver *> \author NAG Ltd. * -*> \date November 2011 +*> \date April 2012 * *> \ingroup auxOTHERcomputational * @@ -140,10 +140,10 @@ * ===================================================================== SUBROUTINE SLAED6( KNITER, ORGATI, RHO, D, Z, FINIT, TAU, INFO ) * -* -- LAPACK computational routine (version 3.4.0) -- +* -- LAPACK computational routine (version 3.4.1) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- -* November 2011 +* April 2012 * * .. Scalar Arguments .. LOGICAL ORGATI @@ -370,15 +370,19 @@ DF = ZERO DDF = ZERO DO 40 I = 1, 3 - TEMP = ONE / ( DSCALE( I )-TAU ) - TEMP1 = ZSCALE( I )*TEMP - TEMP2 = TEMP1*TEMP - TEMP3 = TEMP2*TEMP - TEMP4 = TEMP1 / DSCALE( I ) - FC = FC + TEMP4 - ERRETM = ERRETM + ABS( TEMP4 ) - DF = DF + TEMP2 - DDF = DDF + TEMP3 + IF (DSCALE( I ).NE.TAU) THEN + TEMP = ONE / ( DSCALE( I )-TAU ) + TEMP1 = ZSCALE( I )*TEMP + TEMP2 = TEMP1*TEMP + TEMP3 = TEMP2*TEMP + TEMP4 = TEMP1 / DSCALE( I ) + FC = FC + TEMP4 + ERRETM = ERRETM + ABS( TEMP4 ) + DF = DF + TEMP2 + DDF = DDF + TEMP3 + ELSE + GO TO 60 + END IF 40 CONTINUE F = FINIT + TAU*FC ERRETM = EIGHT*( ABS( FINIT )+ABS( TAU )*ERRETM ) + |