diff --git a/lapack-netlib/SRC/dlanv2.f b/lapack-netlib/SRC/dlanv2.f index d68481f7e..61b016f16 100644 --- a/lapack-netlib/SRC/dlanv2.f +++ b/lapack-netlib/SRC/dlanv2.f @@ -140,13 +140,16 @@ * * .. Parameters .. DOUBLE PRECISION ZERO, HALF, ONE - PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0 ) + PARAMETER ( ZERO = 0.0D+0, HALF = 0.5D+0, ONE = 1.0D+0, + $ TWO = 2.0D0 ) DOUBLE PRECISION MULTPL PARAMETER ( MULTPL = 4.0D+0 ) * .. * .. Local Scalars .. DOUBLE PRECISION AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB, - $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z + $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z, SAFMIN, + $ SAFMN2, SAFMX2 + INTEGER COUNT * .. * .. External Functions .. DOUBLE PRECISION DLAMCH, DLAPY2 @@ -157,7 +160,11 @@ * .. * .. Executable Statements .. * + SAFMIN = DLAMCH( 'S' ) EPS = DLAMCH( 'P' ) + SAFMN2 = DLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / + $ LOG( DLAMCH( 'B' ) ) / TWO ) + SAFMX2 = ONE / SAFMN2 IF( C.EQ.ZERO ) THEN CS = ONE SN = ZERO @@ -212,7 +219,24 @@ * Complex eigenvalues, or real (almost) equal eigenvalues. * Make diagonal elements equal. * + COUNT = 0 SIGMA = B + C + 10 CONTINUE + COUNT = COUNT + 1 + SCALE = MAX( ABS(TEMP), ABS(SIGMA) ) + IF( SCALE.GE.SAFMX2 ) THEN + SIGMA = SIGMA * SAFMN2 + TEMP = TEMP * SAFMN2 + IF (COUNT .LE. 20) + $ GOTO 10 + END IF + IF( SCALE.LE.SAFMN2 ) THEN + SIGMA = SIGMA * SAFMX2 + TEMP = TEMP * SAFMX2 + IF (COUNT .LE. 20) + $ GOTO 10 + END IF + P = HALF*TEMP TAU = DLAPY2( SIGMA, TEMP ) CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) ) SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA ) diff --git a/lapack-netlib/SRC/slanv2.f b/lapack-netlib/SRC/slanv2.f index 1163446fa..e678305f2 100644 --- a/lapack-netlib/SRC/slanv2.f +++ b/lapack-netlib/SRC/slanv2.f @@ -140,13 +140,16 @@ * * .. Parameters .. REAL ZERO, HALF, ONE - PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0 ) + PARAMETER ( ZERO = 0.0E+0, HALF = 0.5E+0, ONE = 1.0E+0, + $ TWO = 2.0E+0 ) REAL MULTPL PARAMETER ( MULTPL = 4.0E+0 ) * .. * .. Local Scalars .. REAL AA, BB, BCMAX, BCMIS, CC, CS1, DD, EPS, P, SAB, - $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z + $ SAC, SCALE, SIGMA, SN1, TAU, TEMP, Z, SAFMIN, + $ SAFMN2, SAFMX2 + INTEGER COUNT * .. * .. External Functions .. REAL SLAMCH, SLAPY2 @@ -157,7 +160,11 @@ * .. * .. Executable Statements .. * + SAFMIN = SLAMCH( 'S' ) EPS = SLAMCH( 'P' ) + SAFMN2 = SLAMCH( 'B' )**INT( LOG( SAFMIN / EPS ) / + $ LOG( SLAMCH( 'B' ) ) / TWO ) + SAFMX2 = ONE / SAFMN2 IF( C.EQ.ZERO ) THEN CS = ONE SN = ZERO @@ -212,7 +219,24 @@ * Complex eigenvalues, or real (almost) equal eigenvalues. * Make diagonal elements equal. * + COUNT = 0 SIGMA = B + C + 10 CONTINUE + COUNT = COUNT + 1 + SCALE = MAX( ABS(TEMP), ABS(SIGMA) ) + IF( SCALE.GE.SAFMX2 ) THEN + SIGMA = SIGMA * SAFMN2 + TEMP = TEMP * SAFMN2 + IF (COUNT .LE. 20) + $ GOTO 10 + END IF + IF( SCALE.LE.SAFMN2 ) THEN + SIGMA = SIGMA * SAFMX2 + TEMP = TEMP * SAFMX2 + IF (COUNT .LE. 20) + $ GOTO 10 + END IF + P = HALF*TEMP TAU = SLAPY2( SIGMA, TEMP ) CS = SQRT( HALF*( ONE+ABS( SIGMA ) / TAU ) ) SN = -( P / ( TAU*CS ) )*SIGN( ONE, SIGMA )