KTH framework for Nek5000 toolboxes; testing version  0.0.1
dlarfg.f
Go to the documentation of this file.
1  SUBROUTINE dlarfg( N, ALPHA, X, INCX, TAU )
2 *
3 * -- LAPACK auxiliary routine (version 3.0) --
4 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5 * Courant Institute, Argonne National Lab, and Rice University
6 * September 30, 1994
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 * DLARFG 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, and x is an (n-1)-element real
26 * 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 * Otherwise 1 <= tau <= 2.
38 *
39 * Arguments
40 * =========
41 *
42 * N (input) INTEGER
43 * The order of the elementary reflector.
44 *
45 * ALPHA (input/output) DOUBLE PRECISION
46 * On entry, the value alpha.
47 * On exit, it is overwritten with the value beta.
48 *
49 * X (input/output) DOUBLE PRECISION array, dimension
50 * (1+(N-2)*abs(INCX))
51 * On entry, the vector x.
52 * On exit, it is overwritten with the vector v.
53 *
54 * INCX (input) INTEGER
55 * The increment between elements of X. INCX > 0.
56 *
57 * TAU (output) DOUBLE PRECISION
58 * The value tau.
59 *
60 * =====================================================================
61 *
62 * .. Parameters ..
63  DOUBLE PRECISION ONE, ZERO
64  parameter( one = 1.0d+0, zero = 0.0d+0 )
65 * ..
66 * .. Local Scalars ..
67  INTEGER J, KNT
68  DOUBLE PRECISION BETA, RSAFMN, SAFMIN, XNORM
69 * ..
70 * .. External Functions ..
71  DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
72  EXTERNAL dlamch, dlapy2, dnrm2
73 * ..
74 * .. Intrinsic Functions ..
75  INTRINSIC abs, sign
76 * ..
77 * .. External Subroutines ..
78  EXTERNAL dscal
79 * ..
80 * .. Executable Statements ..
81 *
82  IF( n.LE.1 ) THEN
83  tau = zero
84  RETURN
85  END IF
86 *
87  xnorm = dnrm2( n-1, x, incx )
88 *
89  IF( xnorm.EQ.zero ) THEN
90 *
91 * H = I
92 *
93  tau = zero
94  ELSE
95 *
96 * general case
97 *
98  beta = -sign( dlapy2( alpha, xnorm ), alpha )
99  safmin = dlamch( 'S' ) / dlamch( 'E' )
100  IF( abs( beta ).LT.safmin ) THEN
101 *
102 * XNORM, BETA may be inaccurate; scale X and recompute them
103 *
104  rsafmn = one / safmin
105  knt = 0
106  10 CONTINUE
107  knt = knt + 1
108  CALL dscal( n-1, rsafmn, x, incx )
109  beta = beta*rsafmn
110  alpha = alpha*rsafmn
111  IF( abs( beta ).LT.safmin )
112  $ GO TO 10
113 *
114 * New BETA is at most 1, at least SAFMIN
115 *
116  xnorm = dnrm2( n-1, x, incx )
117  beta = -sign( dlapy2( alpha, xnorm ), alpha )
118  tau = ( beta-alpha ) / beta
119  CALL dscal( n-1, one / ( alpha-beta ), x, incx )
120 *
121 * If ALPHA is subnormal, it may lose relative accuracy
122 *
123  alpha = beta
124  DO 20 j = 1, knt
125  alpha = alpha*safmin
126  20 CONTINUE
127  ELSE
128  tau = ( beta-alpha ) / beta
129  CALL dscal( n-1, one / ( alpha-beta ), x, incx )
130  alpha = beta
131  END IF
132  END IF
133 *
134  RETURN
135 *
136 * End of DLARFG
137 *
138  END
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
Definition: dlarfg.f:2
subroutine dscal(n, da, dx, incx)
Definition: dscal.f:2