1 SUBROUTINE DLARFGP( N, ALPHA, X, INCX, TAU )
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**T * 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**T ) ,
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
2 *
3 * -- LAPACK auxiliary routine (version 3.3.1) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * -- April 2011 --
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**T * 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**T ) ,
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