Return to zlarfgp.f CVS log

Up to [local] / rpl / lapack / lapack

Ajout des nouveaux fichiers pour Lapack 3.2.2.

    1:       SUBROUTINE ZLARFGP( N, ALPHA, X, INCX, TAU )
    2: *
    3: *  -- LAPACK auxiliary routine (version 3.2.2) --
    4: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
    5: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
    6: *     June 2010
    7: *
    8: *     .. Scalar Arguments ..
    9:       INTEGER            INCX, N
   10:       COMPLEX*16         ALPHA, TAU
   11: *     ..
   12: *     .. Array Arguments ..
   13:       COMPLEX*16         X( * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  ZLARFGP generates a complex elementary reflector H of order n, such
   20: *  that
   21: *
   22: *        H' * ( alpha ) = ( beta ),   H' * H = I.
   23: *             (   x   )   (   0  )
   24: *
   25: *  where alpha and beta are scalars, beta is real and non-negative, and
   26: *  x is an (n-1)-element complex vector.  H is represented in the form
   27: *
   28: *        H = I - tau * ( 1 ) * ( 1 v' ) ,
   29: *                      ( v )
   30: *
   31: *  where tau is a complex scalar and v is a complex (n-1)-element
   32: *  vector. Note that H is not hermitian.
   33: *
   34: *  If the elements of x are all zero and alpha is real, then tau = 0
   35: *  and H is taken to be the unit matrix.
   36: *
   37: *  Arguments
   38: *  =========
   39: *
   40: *  N       (input) INTEGER
   41: *          The order of the elementary reflector.
   42: *
   43: *  ALPHA   (input/output) COMPLEX*16
   44: *          On entry, the value alpha.
   45: *          On exit, it is overwritten with the value beta.
   46: *
   47: *  X       (input/output) COMPLEX*16 array, dimension
   48: *                         (1+(N-2)*abs(INCX))
   49: *          On entry, the vector x.
   50: *          On exit, it is overwritten with the vector v.
   51: *
   52: *  INCX    (input) INTEGER
   53: *          The increment between elements of X. INCX > 0.
   54: *
   55: *  TAU     (output) COMPLEX*16
   56: *          The value tau.
   57: *
   58: *  =====================================================================
   59: *
   60: *     .. Parameters ..
   61:       DOUBLE PRECISION   TWO, ONE, ZERO
   62:       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
   63: *     ..
   64: *     .. Local Scalars ..
   65:       INTEGER            J, KNT
   66:       DOUBLE PRECISION   ALPHI, ALPHR, BETA, BIGNUM, SMLNUM, XNORM
   67:       COMPLEX*16         SAVEALPHA
   68: *     ..
   69: *     .. External Functions ..
   70:       DOUBLE PRECISION   DLAMCH, DLAPY3, DLAPY2, DZNRM2
   71:       COMPLEX*16         ZLADIV
   72:       EXTERNAL           DLAMCH, DLAPY3, DLAPY2, DZNRM2, ZLADIV
   73: *     ..
   74: *     .. Intrinsic Functions ..
   75:       INTRINSIC          ABS, DBLE, DCMPLX, DIMAG, SIGN
   76: *     ..
   77: *     .. External Subroutines ..
   78:       EXTERNAL           ZDSCAL, ZSCAL
   79: *     ..
   80: *     .. Executable Statements ..
   81: *
   82:       IF( N.LE.0 ) THEN
   83:          TAU = ZERO
   84:          RETURN
   85:       END IF
   86: *
   87:       XNORM = DZNRM2( N-1, X, INCX )
   88:       ALPHR = DBLE( ALPHA )
   89:       ALPHI = DIMAG( ALPHA )
   90: *
   91:       IF( XNORM.EQ.ZERO ) THEN
   92: *
   93: *        H  =  [1-alpha/abs(alpha) 0; 0 I], sign chosen so ALPHA >= 0.
   94: *
   95:          IF( ALPHI.EQ.ZERO ) THEN
   96:             IF( ALPHR.GE.ZERO ) THEN
   97: *              When TAU.eq.ZERO, the vector is special-cased to be
   98: *              all zeros in the application routines.  We do not need
   99: *              to clear it.
  100:                TAU = ZERO
  101:             ELSE
  102: *              However, the application routines rely on explicit
  103: *              zero checks when TAU.ne.ZERO, and we must clear X.
  104:                TAU = TWO
  105:                DO J = 1, N-1
  106:                   X( 1 + (J-1)*INCX ) = ZERO
  107:                END DO
  108:                ALPHA = -ALPHA
  109:             END IF
  110:          ELSE
  111: *           Only "reflecting" the diagonal entry to be real and non-negative.
  112:             XNORM = DLAPY2( ALPHR, ALPHI )
  113:             TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
  114:             DO J = 1, N-1
  115:                X( 1 + (J-1)*INCX ) = ZERO
  116:             END DO
  117:             ALPHA = XNORM
  118:          END IF
  119:       ELSE
  120: *
  121: *        general case
  122: *
  123:          BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
  124:          SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
  125:          BIGNUM = ONE / SMLNUM
  126: *
  127:          KNT = 0
  128:          IF( ABS( BETA ).LT.SMLNUM ) THEN
  129: *
  130: *           XNORM, BETA may be inaccurate; scale X and recompute them
  131: *
  132:    10       CONTINUE
  133:             KNT = KNT + 1
  134:             CALL ZDSCAL( N-1, BIGNUM, X, INCX )
  135:             BETA = BETA*BIGNUM
  136:             ALPHI = ALPHI*BIGNUM
  137:             ALPHR = ALPHR*BIGNUM
  138:             IF( ABS( BETA ).LT.SMLNUM )
  139:      $         GO TO 10
  140: *
  141: *           New BETA is at most 1, at least SMLNUM
  142: *
  143:             XNORM = DZNRM2( N-1, X, INCX )
  144:             ALPHA = DCMPLX( ALPHR, ALPHI )
  145:             BETA = SIGN( DLAPY3( ALPHR, ALPHI, XNORM ), ALPHR )
  146:          END IF
  147:          SAVEALPHA = ALPHA
  148:          ALPHA = ALPHA + BETA
  149:          IF( BETA.LT.ZERO ) THEN
  150:             BETA = -BETA
  151:             TAU = -ALPHA / BETA
  152:          ELSE
  153:             ALPHR = ALPHI * (ALPHI/DBLE( ALPHA ))
  154:             ALPHR = ALPHR + XNORM * (XNORM/DBLE( ALPHA ))
  155:             TAU = DCMPLX( ALPHR/BETA, -ALPHI/BETA )
  156:             ALPHA = DCMPLX( -ALPHR, ALPHI )
  157:          END IF
  158:          ALPHA = ZLADIV( DCMPLX( ONE ), ALPHA )
  159: *
  160:          IF ( ABS(TAU).LE.SMLNUM ) THEN
  161: *
  162: *           In the case where the computed TAU ends up being a denormalized number,
  163: *           it loses relative accuracy. This is a BIG problem. Solution: flush TAU 
  164: *           to ZERO (or TWO or whatever makes a nonnegative real number for BETA).
  165: *
  166: *           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
  167: *           (Thanks Pat. Thanks MathWorks.)
  168: *
  169:             ALPHR = DBLE( SAVEALPHA )
  170:             ALPHI = DIMAG( SAVEALPHA )
  171:             IF( ALPHI.EQ.ZERO ) THEN
  172:                IF( ALPHR.GE.ZERO ) THEN
  173:                   TAU = ZERO
  174:                ELSE
  175:                   TAU = TWO
  176:                   DO J = 1, N-1
  177:                      X( 1 + (J-1)*INCX ) = ZERO
  178:                   END DO
  179:                   BETA = -SAVEALPHA
  180:                END IF
  181:             ELSE
  182:                XNORM = DLAPY2( ALPHR, ALPHI )
  183:                TAU = DCMPLX( ONE - ALPHR / XNORM, -ALPHI / XNORM )
  184:                DO J = 1, N-1
  185:                   X( 1 + (J-1)*INCX ) = ZERO
  186:                END DO
  187:                BETA = XNORM
  188:             END IF
  189: *
  190:          ELSE 
  191: *
  192: *           This is the general case.
  193: *
  194:             CALL ZSCAL( N-1, ALPHA, X, INCX )
  195: *
  196:          END IF
  197: *
  198: *        If BETA is subnormal, it may lose relative accuracy
  199: *
  200:          DO 20 J = 1, KNT
  201:             BETA = BETA*SMLNUM
  202:  20      CONTINUE
  203:          ALPHA = BETA
  204:       END IF
  205: *
  206:       RETURN
  207: *
  208: *     End of ZLARFGP
  209: *
  210:       END