Return to dlarfgp.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b DLARFGP generates an elementary reflector (Householder matrix) with non-negatibe beta.
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at 
    6: *            http://www.netlib.org/lapack/explore-html/ 
    7: *
    8: *> \htmlonly
    9: *> Download DLARFGP + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfgp.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfgp.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfgp.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
   22: * 
   23: *       .. Scalar Arguments ..
   24: *       INTEGER            INCX, N
   25: *       DOUBLE PRECISION   ALPHA, TAU
   26: *       ..
   27: *       .. Array Arguments ..
   28: *       DOUBLE PRECISION   X( * )
   29: *       ..
   30: *  
   31: *
   32: *> \par Purpose:
   33: *  =============
   34: *>
   35: *> \verbatim
   36: *>
   37: *> DLARFGP generates a real elementary reflector H of order n, such
   38: *> that
   39: *>
   40: *>       H * ( alpha ) = ( beta ),   H**T * H = I.
   41: *>           (   x   )   (   0  )
   42: *>
   43: *> where alpha and beta are scalars, beta is non-negative, and x is
   44: *> an (n-1)-element real vector.  H is represented in the form
   45: *>
   46: *>       H = I - tau * ( 1 ) * ( 1 v**T ) ,
   47: *>                     ( v )
   48: *>
   49: *> where tau is a real scalar and v is a real (n-1)-element
   50: *> vector.
   51: *>
   52: *> If the elements of x are all zero, then tau = 0 and H is taken to be
   53: *> the unit matrix.
   54: *> \endverbatim
   55: *
   56: *  Arguments:
   57: *  ==========
   58: *
   59: *> \param[in] N
   60: *> \verbatim
   61: *>          N is INTEGER
   62: *>          The order of the elementary reflector.
   63: *> \endverbatim
   64: *>
   65: *> \param[in,out] ALPHA
   66: *> \verbatim
   67: *>          ALPHA is DOUBLE PRECISION
   68: *>          On entry, the value alpha.
   69: *>          On exit, it is overwritten with the value beta.
   70: *> \endverbatim
   71: *>
   72: *> \param[in,out] X
   73: *> \verbatim
   74: *>          X is DOUBLE PRECISION array, dimension
   75: *>                         (1+(N-2)*abs(INCX))
   76: *>          On entry, the vector x.
   77: *>          On exit, it is overwritten with the vector v.
   78: *> \endverbatim
   79: *>
   80: *> \param[in] INCX
   81: *> \verbatim
   82: *>          INCX is INTEGER
   83: *>          The increment between elements of X. INCX > 0.
   84: *> \endverbatim
   85: *>
   86: *> \param[out] TAU
   87: *> \verbatim
   88: *>          TAU is DOUBLE PRECISION
   89: *>          The value tau.
   90: *> \endverbatim
   91: *
   92: *  Authors:
   93: *  ========
   94: *
   95: *> \author Univ. of Tennessee 
   96: *> \author Univ. of California Berkeley 
   97: *> \author Univ. of Colorado Denver 
   98: *> \author NAG Ltd. 
   99: *
  100: *> \date September 2012
  101: *
  102: *> \ingroup doubleOTHERauxiliary
  103: *
  104: *  =====================================================================
  105:       SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
  106: *
  107: *  -- LAPACK auxiliary routine (version 3.4.2) --
  108: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  109: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  110: *     September 2012
  111: *
  112: *     .. Scalar Arguments ..
  113:       INTEGER            INCX, N
  114:       DOUBLE PRECISION   ALPHA, TAU
  115: *     ..
  116: *     .. Array Arguments ..
  117:       DOUBLE PRECISION   X( * )
  118: *     ..
  119: *
  120: *  =====================================================================
  121: *
  122: *     .. Parameters ..
  123:       DOUBLE PRECISION   TWO, ONE, ZERO
  124:       PARAMETER          ( TWO = 2.0D+0, ONE = 1.0D+0, ZERO = 0.0D+0 )
  125: *     ..
  126: *     .. Local Scalars ..
  127:       INTEGER            J, KNT
  128:       DOUBLE PRECISION   BETA, BIGNUM, SAVEALPHA, SMLNUM, XNORM
  129: *     ..
  130: *     .. External Functions ..
  131:       DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
  132:       EXTERNAL           DLAMCH, DLAPY2, DNRM2
  133: *     ..
  134: *     .. Intrinsic Functions ..
  135:       INTRINSIC          ABS, SIGN
  136: *     ..
  137: *     .. External Subroutines ..
  138:       EXTERNAL           DSCAL
  139: *     ..
  140: *     .. Executable Statements ..
  141: *
  142:       IF( N.LE.0 ) THEN
  143:          TAU = ZERO
  144:          RETURN
  145:       END IF
  146: *
  147:       XNORM = DNRM2( N-1, X, INCX )
  148: *
  149:       IF( XNORM.EQ.ZERO ) THEN
  150: *
  151: *        H  =  [+/-1, 0; I], sign chosen so ALPHA >= 0
  152: *
  153:          IF( ALPHA.GE.ZERO ) THEN
  154: *           When TAU.eq.ZERO, the vector is special-cased to be
  155: *           all zeros in the application routines.  We do not need
  156: *           to clear it.
  157:             TAU = ZERO
  158:          ELSE
  159: *           However, the application routines rely on explicit
  160: *           zero checks when TAU.ne.ZERO, and we must clear X.
  161:             TAU = TWO
  162:             DO J = 1, N-1
  163:                X( 1 + (J-1)*INCX ) = 0
  164:             END DO
  165:             ALPHA = -ALPHA
  166:          END IF
  167:       ELSE
  168: *
  169: *        general case
  170: *
  171:          BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
  172:          SMLNUM = DLAMCH( 'S' ) / DLAMCH( 'E' )
  173:          KNT = 0
  174:          IF( ABS( BETA ).LT.SMLNUM ) THEN
  175: *
  176: *           XNORM, BETA may be inaccurate; scale X and recompute them
  177: *
  178:             BIGNUM = ONE / SMLNUM
  179:    10       CONTINUE
  180:             KNT = KNT + 1
  181:             CALL DSCAL( N-1, BIGNUM, X, INCX )
  182:             BETA = BETA*BIGNUM
  183:             ALPHA = ALPHA*BIGNUM
  184:             IF( ABS( BETA ).LT.SMLNUM )
  185:      $         GO TO 10
  186: *
  187: *           New BETA is at most 1, at least SMLNUM
  188: *
  189:             XNORM = DNRM2( N-1, X, INCX )
  190:             BETA = SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
  191:          END IF
  192:          SAVEALPHA = ALPHA
  193:          ALPHA = ALPHA + BETA
  194:          IF( BETA.LT.ZERO ) THEN
  195:             BETA = -BETA
  196:             TAU = -ALPHA / BETA
  197:          ELSE
  198:             ALPHA = XNORM * (XNORM/ALPHA)
  199:             TAU = ALPHA / BETA
  200:             ALPHA = -ALPHA
  201:          END IF
  202: *
  203:          IF ( ABS(TAU).LE.SMLNUM ) THEN
  204: *
  205: *           In the case where the computed TAU ends up being a denormalized number,
  206: *           it loses relative accuracy. This is a BIG problem. Solution: flush TAU 
  207: *           to ZERO. This explains the next IF statement.
  208: *
  209: *           (Bug report provided by Pat Quillen from MathWorks on Jul 29, 2009.)
  210: *           (Thanks Pat. Thanks MathWorks.)
  211: *
  212:             IF( SAVEALPHA.GE.ZERO ) THEN
  213:                TAU = ZERO
  214:             ELSE
  215:                TAU = TWO
  216:                DO J = 1, N-1
  217:                   X( 1 + (J-1)*INCX ) = 0
  218:                END DO
  219:                BETA = -SAVEALPHA
  220:             END IF
  221: *
  222:          ELSE 
  223: *
  224: *           This is the general case.
  225: *
  226:             CALL DSCAL( N-1, ONE / ALPHA, X, INCX )
  227: *
  228:          END IF
  229: *
  230: *        If BETA is subnormal, it may lose relative accuracy
  231: *
  232:          DO 20 J = 1, KNT
  233:             BETA = BETA*SMLNUM
  234:  20      CONTINUE
  235:          ALPHA = BETA
  236:       END IF
  237: *
  238:       RETURN
  239: *
  240: *     End of DLARFGP
  241: *
  242:       END