1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
|
SUBROUTINE SLARFG( N, ALPHA, X, INCX, TAU )
*
* -- LAPACK auxiliary routine (version 3.0) --
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
* Courant Institute, Argonne National Lab, and Rice University
* September 30, 1994
*
* .. Scalar Arguments ..
INTEGER INCX, N
REAL ALPHA, TAU
* ..
* .. Array Arguments ..
REAL X( * )
* ..
*
* Purpose
* =======
*
* SLARFG generates a real elementary reflector H of order n, such
* that
*
* H * ( alpha ) = ( beta ), H' * H = I.
* ( x ) ( 0 )
*
* where alpha and beta are scalars, and x is an (n-1)-element real
* vector. H is represented in the form
*
* H = I - tau * ( 1 ) * ( 1 v' ) ,
* ( v )
*
* where tau is a real scalar and v is a real (n-1)-element
* vector.
*
* If the elements of x are all zero, then tau = 0 and H is taken to be
* the unit matrix.
*
* Otherwise 1 <= tau <= 2.
*
* Arguments
* =========
*
* N (input) INTEGER
* The order of the elementary reflector.
*
* ALPHA (input/output) REAL
* On entry, the value alpha.
* On exit, it is overwritten with the value beta.
*
* X (input/output) REAL array, dimension
* (1+(N-2)*abs(INCX))
* On entry, the vector x.
* On exit, it is overwritten with the vector v.
*
* INCX (input) INTEGER
* The increment between elements of X. INCX > 0.
*
* TAU (output) REAL
* The value tau.
*
* =====================================================================
*
* .. Parameters ..
REAL ONE, ZERO
PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
* ..
* .. Local Scalars ..
INTEGER J, KNT
REAL BETA, RSAFMN, SAFMIN, XNORM
* ..
* .. External Functions ..
REAL SLAMCH, SLAPY2, SNRM2
EXTERNAL SLAMCH, SLAPY2, SNRM2
* ..
* .. Intrinsic Functions ..
INTRINSIC ABS, SIGN
* ..
* .. External Subroutines ..
EXTERNAL SSCAL
* ..
* .. Executable Statements ..
*
IF( N.LE.1 ) THEN
TAU = ZERO
RETURN
END IF
*
XNORM = SNRM2( N-1, X, INCX )
*
IF( XNORM.EQ.ZERO ) THEN
*
* H = I
*
TAU = ZERO
ELSE
*
* general case
*
BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
SAFMIN = SLAMCH( 'S' ) / SLAMCH( 'E' )
IF( ABS( BETA ).LT.SAFMIN ) THEN
*
* XNORM, BETA may be inaccurate; scale X and recompute them
*
RSAFMN = ONE / SAFMIN
KNT = 0
10 CONTINUE
KNT = KNT + 1
CALL SSCAL( N-1, RSAFMN, X, INCX )
BETA = BETA*RSAFMN
ALPHA = ALPHA*RSAFMN
IF( ABS( BETA ).LT.SAFMIN )
$ GO TO 10
*
* New BETA is at most 1, at least SAFMIN
*
XNORM = SNRM2( N-1, X, INCX )
BETA = -SIGN( SLAPY2( ALPHA, XNORM ), ALPHA )
TAU = ( BETA-ALPHA ) / BETA
CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
*
* If ALPHA is subnormal, it may lose relative accuracy
*
ALPHA = BETA
DO 20 J = 1, KNT
ALPHA = ALPHA*SAFMIN
20 CONTINUE
ELSE
TAU = ( BETA-ALPHA ) / BETA
CALL SSCAL( N-1, ONE / ( ALPHA-BETA ), X, INCX )
ALPHA = BETA
END IF
END IF
*
RETURN
*
* End of SLARFG
*
END
|
