Hello, I am posting the fix of Jason on this issue. Thanks a lot for reporting it to us. So Jason fixed this and commited it in the svn repository. Patch below. And we will release this soon. Best wishes, Julien.

=======================================================================================

Author: jason

Date: Sat Mar 14 16:06:19 2009

New Revision: 647

URL:

https://icl.cs.utk.edu/trac/lapack-dev/changeset/647Log: Fix ZLARFP and CLARFP optimizations when the vector is zero.

The x == zero branch ignored complex alphas. The code still functioned, but it scaled the entire zero vector. Now that I think of it, I should scan upwards for the scaling, too. That will be a separate enhancement patch; this is just the bug fix.

This patch also fixes an accidental precision shortening in ZLARFP, which used CMPLX in the dead branch.

Reported by Igor Zhuravlov on the LAPACK web goo thingy, currently at

viewtopic.php?t=924Also lead to finding a related (but different) error in the LAWN and (accepted, still in editing) SISC paper.

Signed-off-by: Jason Riedy <ejr@cs.berkeley.edu>

Modified:

lapack/trunk/SRC/clarfp.f

lapack/trunk/SRC/zlarfp.f

Modified: lapack/trunk/SRC/clarfp.f

URL:

https://icl.cs.utk.edu/trac/lapack-dev/ ... .f?rev=647==============================================================================

--- lapack/trunk/SRC/clarfp.f (original)

+++ lapack/trunk/SRC/clarfp.f Sat Mar 14 16:06:19 2009

@@ -88,7 +88,7 @@

ALPHR = REAL( ALPHA )

ALPHI = AIMAG( ALPHA )

*

- IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN

+ IF( XNORM.EQ.ZERO ) THEN

*

* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.

*

@@ -103,7 +103,7 @@

! zero checks when TAU.ne.ZERO, and we must clear X.

TAU = TWO

DO J = 1, N-1

- X( 1 + (J-1)*INCX ) = 0

+ X( 1 + (J-1)*INCX ) = ZERO

END DO

ALPHA = -ALPHA

END IF

@@ -112,7 +112,7 @@

XNORM = SLAPY2( ALPHR, ALPHI )

TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )

DO J = 1, N-1

- X( 1 + (J-1)*INCX ) = 0

+ X( 1 + (J-1)*INCX ) = ZERO

END DO

ALPHA = XNORM

END IF

Modified: lapack/trunk/SRC/zlarfp.f

URL:

https://icl.cs.utk.edu/trac/lapack-dev/ ... .f?rev=647==============================================================================

--- lapack/trunk/SRC/zlarfp.f (original)

+++ lapack/trunk/SRC/zlarfp.f Sat Mar 14 16:06:19 2009

@@ -88,7 +88,7 @@

ALPHR = DBLE( ALPHA )

ALPHI = DIMAG( ALPHA )

*

- IF( XNORM.EQ.ZERO .AND. ALPHI.EQ.ZERO ) THEN

+ IF( XNORM.EQ.ZERO ) THEN

*

* H = [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.

*

@@ -103,16 +103,16 @@

! zero checks when TAU.ne.ZERO, and we must clear X.

TAU = TWO

DO J = 1, N-1

- X( 1 + (J-1)*INCX ) = 0

+ X( 1 + (J-1)*INCX ) = ZERO

END DO

ALPHA = -ALPHA

END IF

ELSE

! Only "reflecting" the diagonal entry to be real and non-negative.

XNORM = DLAPY2( ALPHR, ALPHI )

- TAU = CMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )

+ TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )

DO J = 1, N-1

- X( 1 + (J-1)*INCX ) = 0

+ X( 1 + (J-1)*INCX ) = ZERO

END DO

ALPHA = XNORM

END IF