2
3
4
5
6
7
8 INTEGER INFO, KL, KU, LDA, M, N
9
10
11 INTEGER ISEED( 4 )
12 REAL D( * )
13 COMPLEX A( LDA, * ), WORK( * )
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 COMPLEX ZERO, ONE
66 parameter( zero = ( 0.0e+0, 0.0e+0 ),
67 $ one = ( 1.0e+0, 0.0e+0 ) )
68
69
70 INTEGER I, J
71 REAL WN
72 COMPLEX TAU, WA, WB
73
74
75 EXTERNAL cgemv, cgerc, clacgv,
clarnv, cscal, xerbla
76
77
78 INTRINSIC abs,
max,
min, real
79
80
81 REAL SCNRM2
82 EXTERNAL scnrm2
83
84
85
86
87
88 info = 0
89 IF( m.LT.0 ) THEN
90 info = -1
91 ELSE IF( n.LT.0 ) THEN
92 info = -2
93 ELSE IF( kl.LT.0 .OR. kl.GT.m-1 ) THEN
94 info = -3
95 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
96 info = -4
97 ELSE IF( lda.LT.
max( 1, m ) )
THEN
98 info = -7
99 END IF
100 IF( info.LT.0 ) THEN
101 CALL xerbla( 'CLAGGE', -info )
102 RETURN
103 END IF
104
105
106
107 DO 20 j = 1, n
108 DO 10 i = 1, m
109 a( i, j ) = zero
110 10 CONTINUE
111 20 CONTINUE
112 DO 30 i = 1,
min( m, n )
113 a( i, i ) = d( i )
114 30 CONTINUE
115
116
117
118 DO 40 i =
min( m, n ), 1, -1
119 IF( i.LT.m ) THEN
120
121
122
123 CALL clarnv( 3, iseed, m-i+1, work )
124 wn = scnrm2( m-i+1, work, 1 )
125 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
126 IF( wn.EQ.zero ) THEN
127 tau = zero
128 ELSE
129 wb = work( 1 ) + wa
130 CALL cscal( m-i, one / wb, work( 2 ), 1 )
131 work( 1 ) = one
132 tau = real( wb / wa )
133 END IF
134
135
136
137 CALL cgemv( 'Conjugate transpose', m-i+1, n-i+1, one,
138 $ a( i, i ), lda, work, 1, zero, work( m+1 ), 1 )
139 CALL cgerc( m-i+1, n-i+1, -tau, work, 1, work( m+1 ), 1,
140 $ a( i, i ), lda )
141 END IF
142 IF( i.LT.n ) THEN
143
144
145
146 CALL clarnv( 3, iseed, n-i+1, work )
147 wn = scnrm2( n-i+1, work, 1 )
148 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
149 IF( wn.EQ.zero ) THEN
150 tau = zero
151 ELSE
152 wb = work( 1 ) + wa
153 CALL cscal( n-i, one / wb, work( 2 ), 1 )
154 work( 1 ) = one
155 tau = real( wb / wa )
156 END IF
157
158
159
160 CALL cgemv( 'No transpose', m-i+1, n-i+1, one, a( i, i ),
161 $ lda, work, 1, zero, work( n+1 ), 1 )
162 CALL cgerc( m-i+1, n-i+1, -tau, work( n+1 ), 1, work, 1,
163 $ a( i, i ), lda )
164 END IF
165 40 CONTINUE
166
167
168
169
170 DO 70 i = 1,
max( m-1-kl, n-1-ku )
171 IF( kl.LE.ku ) THEN
172
173
174
175 IF( i.LE.
min( m-1-kl, n ) )
THEN
176
177
178
179 wn = scnrm2( m-kl-i+1, a( kl+i, i ), 1 )
180 wa = ( wn / abs( a( kl+i, i ) ) )*a( kl+i, i )
181 IF( wn.EQ.zero ) THEN
182 tau = zero
183 ELSE
184 wb = a( kl+i, i ) + wa
185 CALL cscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
186 a( kl+i, i ) = one
187 tau = real( wb / wa )
188 END IF
189
190
191
192 CALL cgemv( 'Conjugate transpose', m-kl-i+1, n-i, one,
193 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
194 $ work, 1 )
195 CALL cgerc( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work,
196 $ 1, a( kl+i, i+1 ), lda )
197 a( kl+i, i ) = -wa
198 END IF
199
200 IF( i.LE.
min( n-1-ku, m ) )
THEN
201
202
203
204 wn = scnrm2( n-ku-i+1, a( i, ku+i ), lda )
205 wa = ( wn / abs( a( i, ku+i ) ) )*a( i, ku+i )
206 IF( wn.EQ.zero ) THEN
207 tau = zero
208 ELSE
209 wb = a( i, ku+i ) + wa
210 CALL cscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
211 a( i, ku+i ) = one
212 tau = real( wb / wa )
213 END IF
214
215
216
217 CALL clacgv( n-ku-i+1, a( i, ku+i ), lda )
218 CALL cgemv( 'No transpose', m-i, n-ku-i+1, one,
219 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
220 $ work, 1 )
221 CALL cgerc( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
222 $ lda, a( i+1, ku+i ), lda )
223 a( i, ku+i ) = -wa
224 END IF
225 ELSE
226
227
228
229
230 IF( i.LE.
min( n-1-ku, m ) )
THEN
231
232
233
234 wn = scnrm2( n-ku-i+1, a( i, ku+i ), lda )
235 wa = ( wn / abs( a( i, ku+i ) ) )*a( i, ku+i )
236 IF( wn.EQ.zero ) THEN
237 tau = zero
238 ELSE
239 wb = a( i, ku+i ) + wa
240 CALL cscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
241 a( i, ku+i ) = one
242 tau = real( wb / wa )
243 END IF
244
245
246
247 CALL clacgv( n-ku-i+1, a( i, ku+i ), lda )
248 CALL cgemv( 'No transpose', m-i, n-ku-i+1, one,
249 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
250 $ work, 1 )
251 CALL cgerc( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
252 $ lda, a( i+1, ku+i ), lda )
253 a( i, ku+i ) = -wa
254 END IF
255
256 IF( i.LE.
min( m-1-kl, n ) )
THEN
257
258
259
260 wn = scnrm2( m-kl-i+1, a( kl+i, i ), 1 )
261 wa = ( wn / abs( a( kl+i, i ) ) )*a( kl+i, i )
262 IF( wn.EQ.zero ) THEN
263 tau = zero
264 ELSE
265 wb = a( kl+i, i ) + wa
266 CALL cscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
267 a( kl+i, i ) = one
268 tau = real( wb / wa )
269 END IF
270
271
272
273 CALL cgemv( 'Conjugate transpose', m-kl-i+1, n-i, one,
274 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
275 $ work, 1 )
276 CALL cgerc( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work,
277 $ 1, a( kl+i, i+1 ), lda )
278 a( kl+i, i ) = -wa
279 END IF
280 END IF
281
282 DO 50 j = kl + i + 1, m
283 a( j, i ) = zero
284 50 CONTINUE
285
286 DO 60 j = ku + i + 1, n
287 a( i, j ) = zero
288 60 CONTINUE
289 70 CONTINUE
290 RETURN
291
292
293
subroutine clarnv(idist, iseed, n, x)