Return to dlarfgp.f CVS log

Up to [local] / rpl / lapack / lapack

Ajout des nouveaux fichiers pour Lapack 3.2.2.

    1:       SUBROUTINE DLARFGP( 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:       DOUBLE PRECISION   ALPHA, TAU
   11: *     ..
   12: *     .. Array Arguments ..
   13:       DOUBLE PRECISION   X( * )
   14: *     ..
   15: *
   16: *  Purpose
   17: *  =======
   18: *
   19: *  DLARFGP generates a real 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 non-negative, and x is
   26: *  an (n-1)-element real vector.  H is represented in the form
   27: *
   28: *        H = I - tau * ( 1 ) * ( 1 v' ) ,
   29: *                      ( v )
   30: *
   31: *  where tau is a real scalar and v is a real (n-1)-element
   32: *  vector.
   33: *
   34: *  If the elements of x are all zero, then tau = 0 and H is taken to be
   35: *  the unit matrix.
   36: *
   37: *  Arguments
   38: *  =========
   39: *
   40: *  N       (input) INTEGER
   41: *          The order of the elementary reflector.
   42: *
   43: *  ALPHA   (input/output) DOUBLE PRECISION
   44: *          On entry, the value alpha.
   45: *          On exit, it is overwritten with the value beta.
   46: *
   47: *  X       (input/output) DOUBLE PRECISION 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) DOUBLE PRECISION
   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   BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
   67: *     ..
   68: *     .. External Functions ..
   69:       DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
   70:       EXTERNAL           DLAMCH, DLAPY2, DNRM2
   71: *     ..
   72: *     .. Intrinsic Functions ..
   73:       INTRINSIC          ABS, SIGN
   74: *     ..
   75: *     .. External Subroutines ..
   76:       EXTERNAL           DSCAL
   77: *     ..
   78: *     .. Executable Statements ..
   79: *
   80:       IF( N.LE.0 ) THEN
   81:          TAU = ZERO
   82:          RETURN
   83:       END IF
   84: *
   85:       XNORM = DNRM2( N-1, X, INCX )
   86: *
   87:       IF( XNORM.EQ.ZERO ) THEN
   88: *
   89: *        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0
   90: *
   91:          IF( ALPHA.GE.ZERO ) THEN
   92: *           When TAU.eq.ZERO, the vector is special-cased to be
   93: *           all zeros in the application routines.  We do not need
   94: *           to clear it.
   95:             TAU = ZERO
   96:          ELSE
   97: *           However, the application routines rely on explicit
   98: *           zero checks when TAU.ne.ZERO, and we must clear X.
   99:             TAU = TWO
  100:             DO J = 1, N-1
  101:                X( 1 + (J-1)*INCX ) = 0
  102:             END DO
  103:             ALPHA = -ALPHA
  104:          END IF
  105:       ELSE
  106: *
  107: *        general case
  108: *
  109:          BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
  110:          SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
  111:          KNT = 0
  112:          IF( ABS( BETA ).LT.SMLNUM ) THEN
  113: *
  114: *           XNORM, BETA may be inaccurate; scale X and recompute them
  115: *
  116:             BIGNUM = ONE / SMLNUM
  117:    10       CONTINUE
  118:             KNT = KNT + 1
  119:             CALL DSCAL( N-1, BIGNUM, X, INCX )
  120:             BETA = BETA*BIGNUM
  121:             ALPHA = ALPHA*BIGNUM
  122:             IF( ABS( BETA ).LT.SMLNUM )
  123:      $         GO TO 10
  124: *
  125: *           New BETA is at most 1, at least SMLNUM
  126: *
  127:             XNORM = DNRM2( N-1, X, INCX )
  128:             BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
  129:          END IF
  130:          SAVEALPHA = ALPHA
  131:          ALPHA = ALPHA + BETA
  132:          IF( BETA.LT.ZERO ) THEN
  133:             BETA = -BETA
  134:             TAU = -ALPHA / BETA
  135:          ELSE
  136:             ALPHA = XNORM * (XNORM/ALPHA)
  137:             TAU = ALPHA / BETA
  138:             ALPHA = -ALPHA
  139:          END IF
  140: *
  141:          IF ( ABS(TAU).LE.SMLNUM ) THEN
  142: *
  143: *           In the case where the computed TAU ends up being a denormalized number,
  144: *           it loses relative accuracy. This is a BIG problem. Solution: flush TAU 
  145: *           to ZERO. This explains the next IF statement.
  146: *
  147: *           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
  148: *           (Thanks Pat. Thanks MathWorks.)
  149: *
  150:             IF( SAVEALPHA.GE.ZERO ) THEN
  151:                TAU = ZERO
  152:             ELSE
  153:                TAU = TWO
  154:                DO J = 1, N-1
  155:                   X( 1 + (J-1)*INCX ) = 0
  156:                END DO
  157:                BETA = -SAVEALPHA
  158:             END IF
  159: *
  160:          ELSE 
  161: *
  162: *           This is the general case.
  163: *
  164:             CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
  165: *
  166:          END IF
  167: *
  168: *        If BETA is subnormal, it may lose relative accuracy
  169: *
  170:          DO 20 J = 1, KNT
  171:             BETA = BETA*SMLNUM
  172:  20      CONTINUE
  173:          ALPHA = BETA
  174:       END IF
  175: *
  176:       RETURN
  177: *
  178: *     End of DLARFGP
  179: *
  180:       END