github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/internal/testdata/netlib/dlarfg.f (about)

     1  *> \brief \b DLARFG generates an elementary reflector (Householder matrix).
     2  *
     3  *  =========== DOCUMENTATION ===========
     4  *
     5  * Online html documentation available at 
     6  *            http://www.netlib.org/lapack/explore-html/ 
     7  *
     8  *> \htmlonly
     9  *> Download DLARFG + dependencies 
    10  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlarfg.f"> 
    11  *> [TGZ]</a> 
    12  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlarfg.f"> 
    13  *> [ZIP]</a> 
    14  *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlarfg.f"> 
    15  *> [TXT]</a>
    16  *> \endhtmlonly 
    17  *
    18  *  Definition:
    19  *  ===========
    20  *
    21  *       SUBROUTINE DLARFG( 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  *> DLARFG 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, and x is an (n-1)-element real
    44  *> 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  *>
    55  *> Otherwise  1 <= tau <= 2.
    56  *> \endverbatim
    57  *
    58  *  Arguments:
    59  *  ==========
    60  *
    61  *> \param[in] N
    62  *> \verbatim
    63  *>          N is INTEGER
    64  *>          The order of the elementary reflector.
    65  *> \endverbatim
    66  *>
    67  *> \param[in,out] ALPHA
    68  *> \verbatim
    69  *>          ALPHA is DOUBLE PRECISION
    70  *>          On entry, the value alpha.
    71  *>          On exit, it is overwritten with the value beta.
    72  *> \endverbatim
    73  *>
    74  *> \param[in,out] X
    75  *> \verbatim
    76  *>          X is DOUBLE PRECISION array, dimension
    77  *>                         (1+(N-2)*abs(INCX))
    78  *>          On entry, the vector x.
    79  *>          On exit, it is overwritten with the vector v.
    80  *> \endverbatim
    81  *>
    82  *> \param[in] INCX
    83  *> \verbatim
    84  *>          INCX is INTEGER
    85  *>          The increment between elements of X. INCX > 0.
    86  *> \endverbatim
    87  *>
    88  *> \param[out] TAU
    89  *> \verbatim
    90  *>          TAU is DOUBLE PRECISION
    91  *>          The value tau.
    92  *> \endverbatim
    93  *
    94  *  Authors:
    95  *  ========
    96  *
    97  *> \author Univ. of Tennessee 
    98  *> \author Univ. of California Berkeley 
    99  *> \author Univ. of Colorado Denver 
   100  *> \author NAG Ltd. 
   101  *
   102  *> \date September 2012
   103  *
   104  *> \ingroup doubleOTHERauxiliary
   105  *
   106  *  =====================================================================
   107        SUBROUTINE DLARFG( N, ALPHA, X, INCX, TAU )
   108  *
   109  *  -- LAPACK auxiliary routine (version 3.4.2) --
   110  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
   111  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   112  *     September 2012
   113  *
   114  *     .. Scalar Arguments ..
   115        INTEGER            INCX, N
   116        DOUBLE PRECISION   ALPHA, TAU
   117  *     ..
   118  *     .. Array Arguments ..
   119        DOUBLE PRECISION   X( * )
   120  *     ..
   121  *
   122  *  =====================================================================
   123  *
   124  *     .. Parameters ..
   125        DOUBLE PRECISION   ONE, ZERO
   126        PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
   127  *     ..
   128  *     .. Local Scalars ..
   129        INTEGER            J, KNT
   130        DOUBLE PRECISION   BETA, RSAFMN, SAFMIN, XNORM
   131  *     ..
   132  *     .. External Functions ..
   133        DOUBLE PRECISION   DLAMCH, DLAPY2, DNRM2
   134        EXTERNAL           DLAMCH, DLAPY2, DNRM2
   135  *     ..
   136  *     .. Intrinsic Functions ..
   137        INTRINSIC          ABS, SIGN
   138  *     ..
   139  *     .. External Subroutines ..
   140        EXTERNAL           DSCAL
   141  *     ..
   142  *     .. Executable Statements ..
   143  *
   144        IF( N.LE.1 ) THEN
   145           TAU = ZERO
   146           RETURN
   147        END IF
   148  *
   149        XNORM = DNRM2( N-1, X, INCX )
   150  *
   151        IF( XNORM.EQ.ZERO ) THEN
   152  *
   153  *        H  =  I
   154  *
   155           TAU = ZERO
   156        ELSE
   157  *
   158  *        general case
   159  *
   160           BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
   161           SAFMIN = DLAMCH( 'S' ) / DLAMCH( 'E' )
   162           KNT = 0
   163           IF( ABS( BETA ).LT.SAFMIN ) THEN
   164  *
   165  *           XNORM, BETA may be inaccurate; scale X and recompute them
   166  *
   167              RSAFMN = ONE / SAFMIN
   168     10       CONTINUE
   169              KNT = KNT + 1
   170              CALL DSCAL( N-1, RSAFMN, X, INCX )
   171              BETA = BETA*RSAFMN
   172              ALPHA = ALPHA*RSAFMN
   173              IF( ABS( BETA ).LT.SAFMIN )
   174       $         GO TO 10
   175  *
   176  *           New BETA is at most 1, at least SAFMIN
   177  *
   178              XNORM = DNRM2( N-1, X, INCX )
   179              BETA = -SIGN( DLAPY2( ALPHA, XNORM ), ALPHA )
   180           END IF
   181           TAU = ( BETA-ALPHA ) / BETA
   182           CALL DSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
   183  *
   184  *        If ALPHA is subnormal, it may lose relative accuracy
   185  *
   186           DO 20 J = 1, KNT
   187              BETA = BETA*SAFMIN
   188   20      CONTINUE
   189           ALPHA = BETA
   190        END IF
   191  *
   192        RETURN
   193  *
   194  *     End of DLARFG
   195  *
   196        END