114
115
116
117
118
119
120 INTEGER INFO, KL, KU, LDA, M, N
121
122
123 INTEGER ISEED( 4 )
124 DOUBLE PRECISION D( * )
125 COMPLEX*16 A( LDA, * ), WORK( * )
126
127
128
129
130
131 COMPLEX*16 ZERO, ONE
132 parameter( zero = ( 0.0d+0, 0.0d+0 ),
133 $ one = ( 1.0d+0, 0.0d+0 ) )
134
135
136 INTEGER I, J
137 DOUBLE PRECISION WN
138 COMPLEX*16 TAU, WA, WB
139
140
142
143
144 INTRINSIC abs, dble, max, min
145
146
147 DOUBLE PRECISION DZNRM2
149
150
151
152
153
154 info = 0
155 IF( m.LT.0 ) THEN
156 info = -1
157 ELSE IF( n.LT.0 ) THEN
158 info = -2
159 ELSE IF( kl.LT.0 .OR. kl.GT.m-1 ) THEN
160 info = -3
161 ELSE IF( ku.LT.0 .OR. ku.GT.n-1 ) THEN
162 info = -4
163 ELSE IF( lda.LT.max( 1, m ) ) THEN
164 info = -7
165 END IF
166 IF( info.LT.0 ) THEN
167 CALL xerbla(
'ZLAGGE', -info )
168 RETURN
169 END IF
170
171
172
173 DO 20 j = 1, n
174 DO 10 i = 1, m
175 a( i, j ) = zero
176 10 CONTINUE
177 20 CONTINUE
178 DO 30 i = 1, min( m, n )
179 a( i, i ) = d( i )
180 30 CONTINUE
181
182
183
184 IF(( kl .EQ. 0 ).AND.( ku .EQ. 0)) RETURN
185
186
187
188 DO 40 i = min( m, n ), 1, -1
189 IF( i.LT.m ) THEN
190
191
192
193 CALL zlarnv( 3, iseed, m-i+1, work )
194 wn =
dznrm2( m-i+1, work, 1 )
195 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
196 IF( wn.EQ.zero ) THEN
197 tau = zero
198 ELSE
199 wb = work( 1 ) + wa
200 CALL zscal( m-i, one / wb, work( 2 ), 1 )
201 work( 1 ) = one
202 tau = dble( wb / wa )
203 END IF
204
205
206
207 CALL zgemv(
'Conjugate transpose', m-i+1, n-i+1, one,
208 $ a( i, i ), lda, work, 1, zero, work( m+1 ), 1 )
209 CALL zgerc( m-i+1, n-i+1, -tau, work, 1, work( m+1 ), 1,
210 $ a( i, i ), lda )
211 END IF
212 IF( i.LT.n ) THEN
213
214
215
216 CALL zlarnv( 3, iseed, n-i+1, work )
217 wn =
dznrm2( n-i+1, work, 1 )
218 wa = ( wn / abs( work( 1 ) ) )*work( 1 )
219 IF( wn.EQ.zero ) THEN
220 tau = zero
221 ELSE
222 wb = work( 1 ) + wa
223 CALL zscal( n-i, one / wb, work( 2 ), 1 )
224 work( 1 ) = one
225 tau = dble( wb / wa )
226 END IF
227
228
229
230 CALL zgemv(
'No transpose', m-i+1, n-i+1, one, a( i, i ),
231 $ lda, work, 1, zero, work( n+1 ), 1 )
232 CALL zgerc( m-i+1, n-i+1, -tau, work( n+1 ), 1, work, 1,
233 $ a( i, i ), lda )
234 END IF
235 40 CONTINUE
236
237
238
239
240 DO 70 i = 1, max( m-1-kl, n-1-ku )
241 IF( kl.LE.ku ) THEN
242
243
244
245 IF( i.LE.min( m-1-kl, n ) ) THEN
246
247
248
249 wn =
dznrm2( m-kl-i+1, a( kl+i, i ), 1 )
250 wa = ( wn / abs( a( kl+i, i ) ) )*a( kl+i, i )
251 IF( wn.EQ.zero ) THEN
252 tau = zero
253 ELSE
254 wb = a( kl+i, i ) + wa
255 CALL zscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
256 a( kl+i, i ) = one
257 tau = dble( wb / wa )
258 END IF
259
260
261
262 CALL zgemv(
'Conjugate transpose', m-kl-i+1, n-i, one,
263 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
264 $ work, 1 )
265 CALL zgerc( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work,
266 $ 1, a( kl+i, i+1 ), lda )
267 a( kl+i, i ) = -wa
268 END IF
269
270 IF( i.LE.min( n-1-ku, m ) ) THEN
271
272
273
274 wn =
dznrm2( n-ku-i+1, a( i, ku+i ), lda )
275 wa = ( wn / abs( a( i, ku+i ) ) )*a( i, ku+i )
276 IF( wn.EQ.zero ) THEN
277 tau = zero
278 ELSE
279 wb = a( i, ku+i ) + wa
280 CALL zscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
281 a( i, ku+i ) = one
282 tau = dble( wb / wa )
283 END IF
284
285
286
287 CALL zlacgv( n-ku-i+1, a( i, ku+i ), lda )
288 CALL zgemv(
'No transpose', m-i, n-ku-i+1, one,
289 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
290 $ work, 1 )
291 CALL zgerc( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
292 $ lda, a( i+1, ku+i ), lda )
293 a( i, ku+i ) = -wa
294 END IF
295 ELSE
296
297
298
299
300 IF( i.LE.min( n-1-ku, m ) ) THEN
301
302
303
304 wn =
dznrm2( n-ku-i+1, a( i, ku+i ), lda )
305 wa = ( wn / abs( a( i, ku+i ) ) )*a( i, ku+i )
306 IF( wn.EQ.zero ) THEN
307 tau = zero
308 ELSE
309 wb = a( i, ku+i ) + wa
310 CALL zscal( n-ku-i, one / wb, a( i, ku+i+1 ), lda )
311 a( i, ku+i ) = one
312 tau = dble( wb / wa )
313 END IF
314
315
316
317 CALL zlacgv( n-ku-i+1, a( i, ku+i ), lda )
318 CALL zgemv(
'No transpose', m-i, n-ku-i+1, one,
319 $ a( i+1, ku+i ), lda, a( i, ku+i ), lda, zero,
320 $ work, 1 )
321 CALL zgerc( m-i, n-ku-i+1, -tau, work, 1, a( i, ku+i ),
322 $ lda, a( i+1, ku+i ), lda )
323 a( i, ku+i ) = -wa
324 END IF
325
326 IF( i.LE.min( m-1-kl, n ) ) THEN
327
328
329
330 wn =
dznrm2( m-kl-i+1, a( kl+i, i ), 1 )
331 wa = ( wn / abs( a( kl+i, i ) ) )*a( kl+i, i )
332 IF( wn.EQ.zero ) THEN
333 tau = zero
334 ELSE
335 wb = a( kl+i, i ) + wa
336 CALL zscal( m-kl-i, one / wb, a( kl+i+1, i ), 1 )
337 a( kl+i, i ) = one
338 tau = dble( wb / wa )
339 END IF
340
341
342
343 CALL zgemv(
'Conjugate transpose', m-kl-i+1, n-i, one,
344 $ a( kl+i, i+1 ), lda, a( kl+i, i ), 1, zero,
345 $ work, 1 )
346 CALL zgerc( m-kl-i+1, n-i, -tau, a( kl+i, i ), 1, work,
347 $ 1, a( kl+i, i+1 ), lda )
348 a( kl+i, i ) = -wa
349 END IF
350 END IF
351
352 IF (i .LE. n) THEN
353 DO 50 j = kl + i + 1, m
354 a( j, i ) = zero
355 50 CONTINUE
356 END IF
357
358 IF (i .LE. m) THEN
359 DO 60 j = ku + i + 1, n
360 a( i, j ) = zero
361 60 CONTINUE
362 END IF
363 70 CONTINUE
364 RETURN
365
366
367
subroutine xerbla(srname, info)
subroutine zgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
ZGEMV
subroutine zgerc(m, n, alpha, x, incx, y, incy, a, lda)
ZGERC
subroutine zlacgv(n, x, incx)
ZLACGV conjugates a complex vector.
subroutine zlarnv(idist, iseed, n, x)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
real(wp) function dznrm2(n, x, incx)
DZNRM2
subroutine zscal(n, za, zx, incx)
ZSCAL