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
  http://icl.cs.utk.edu/lapack-forum/viewtopic.php?t=924
Also lead to finding a related (but different) error in the LAWN and
(accepted, still in editing) SISC paper.

Signed-off-by: Jason Riedy <[email protected]>

Author: jason

SRC/clarfp.f

@@ -88,7 +88,7 @@ SUBROUTINE CLARFP( N, ALPHA, X, INCX, TAU )
       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 @@ SUBROUTINE CLARFP( N, ALPHA, X, INCX, TAU )
 !              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 @@ SUBROUTINE CLARFP( N, ALPHA, X, INCX, TAU )
             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

SRC/zlarfp.f

@@ -88,7 +88,7 @@ SUBROUTINE ZLARFP( N, ALPHA, X, INCX, TAU )
       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 @@ SUBROUTINE ZLARFP( N, ALPHA, X, INCX, TAU )
 !              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