Rescale input vector more often to minimize relative error (Reference-LAPACK PR 981)

This commit is contained in:
Martin Kroeker 2024-02-05 15:35:24 +01:00 committed by GitHub
parent a4fde2c5ac
commit 479e4af089
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
2 changed files with 20 additions and 40 deletions

View File

@ -148,33 +148,23 @@
ALPHR = REAL( ALPHA )
ALPHI = AIMAG( ALPHA )
*
IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
IF( XNORM.LE.EPS*ABS(ALPHA) .AND. ALPHI.EQ.ZERO ) THEN
*
* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
*
IF( ALPHI.EQ.ZERO ) THEN
IF( ALPHR.GE.ZERO ) THEN
* When TAU.eq.ZERO, the vector is special-cased to be
* all zeros in the application routines. We do not need
* to clear it.
TAU = ZERO
ELSE
* However, the application routines rely on explicit
* zero checks when TAU.ne.ZERO, and we must clear X.
TAU = TWO
DO J = 1, N-1
X( 1 + (J-1)*INCX ) = ZERO
END DO
ALPHA = -ALPHA
END IF
IF( ALPHR.GE.ZERO ) THEN
* When TAU.eq.ZERO, the vector is special-cased to be
* all zeros in the application routines. We do not need
* to clear it.
TAU = ZERO
ELSE
* Only "reflecting" the diagonal entry to be real and non-negative.
XNORM = SLAPY2( ALPHR, ALPHI )
TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
* However, the application routines rely on explicit
* zero checks when TAU.ne.ZERO, and we must clear X.
TAU = TWO
DO J = 1, N-1
X( 1 + (J-1)*INCX ) = ZERO
END DO
ALPHA = XNORM
ALPHA = -ALPHA
END IF
ELSE
*

View File

@ -148,33 +148,23 @@
ALPHR = DBLE( ALPHA )
ALPHI = DIMAG( ALPHA )
*
IF( XNORM.LE.EPS*ABS(ALPHA) ) THEN
IF( XNORM.LE.EPS*ABS(ALPHA) .AND. ALPHI.EQ.ZERO ) THEN
*
* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
*
IF( ALPHI.EQ.ZERO ) THEN
IF( ALPHR.GE.ZERO ) THEN
* When TAU.eq.ZERO, the vector is special-cased to be
* all zeros in the application routines. We do not need
* to clear it.
TAU = ZERO
ELSE
* However, the application routines rely on explicit
* zero checks when TAU.ne.ZERO, and we must clear X.
TAU = TWO
DO J = 1, N-1
X( 1 + (J-1)*INCX ) = ZERO
END DO
ALPHA = -ALPHA
END IF
IF( ALPHR.GE.ZERO ) THEN
* When TAU.eq.ZERO, the vector is special-cased to be
* all zeros in the application routines. We do not need
* to clear it.
TAU = ZERO
ELSE
* Only "reflecting" the diagonal entry to be real and non-negative.
XNORM = DLAPY2( ALPHR, ALPHI )
TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
* However, the application routines rely on explicit
* zero checks when TAU.ne.ZERO, and we must clear X.
TAU = TWO
DO J = 1, N-1
X( 1 + (J-1)*INCX ) = ZERO
END DO
ALPHA = XNORM
ALPHA = -ALPHA
END IF
ELSE
*