179 implicit none
180
181
182
183
184
185
186 CHARACTER JOBVL, JOBVR
187 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
188
189
190 DOUBLE PRECISION RWORK( * )
191 COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
192 $ W( * ), WORK( * )
193
194
195
196
197
198 DOUBLE PRECISION ZERO, ONE
199 parameter( zero = 0.0d0, one = 1.0d0 )
200
201
202 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
203 CHARACTER SIDE
204 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
205 $ IWRK, K, LWORK_TREVC, MAXWRK, MINWRK, NOUT
206 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
207 COMPLEX*16 TMP
208
209
210 LOGICAL SELECT( 1 )
211 DOUBLE PRECISION DUM( 1 )
212
213
217
218
219 LOGICAL LSAME
220 INTEGER IDAMAX, ILAENV
221 DOUBLE PRECISION DLAMCH, DZNRM2, ZLANGE
224
225
226 INTRINSIC dble, dcmplx, conjg, aimag, max, sqrt
227
228
229
230
231
232 info = 0
233 lquery = ( lwork.EQ.-1 )
234 wantvl =
lsame( jobvl,
'V' )
235 wantvr =
lsame( jobvr,
'V' )
236 IF( ( .NOT.wantvl ) .AND. ( .NOT.
lsame( jobvl,
'N' ) ) )
THEN
237 info = -1
238 ELSE IF( ( .NOT.wantvr ) .AND.
239 $ ( .NOT.
lsame( jobvr,
'N' ) ) )
THEN
240 info = -2
241 ELSE IF( n.LT.0 ) THEN
242 info = -3
243 ELSE IF( lda.LT.max( 1, n ) ) THEN
244 info = -5
245 ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
246 info = -8
247 ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
248 info = -10
249 END IF
250
251
252
253
254
255
256
257
258
259
260
261
262 IF( info.EQ.0 ) THEN
263 IF( n.EQ.0 ) THEN
264 minwrk = 1
265 maxwrk = 1
266 ELSE
267 maxwrk = n + n*
ilaenv( 1,
'ZGEHRD',
' ', n, 1, n, 0 )
268 minwrk = 2*n
269 IF( wantvl ) THEN
270 maxwrk = max( maxwrk, n + ( n - 1 )*
ilaenv( 1,
271 $ 'ZUNGHR',
272 $ ' ', n, 1, n, -1 ) )
273 CALL ztrevc3(
'L',
'B',
SELECT, n, a, lda,
274 $ vl, ldvl, vr, ldvr,
275 $ n, nout, work, -1, rwork, -1, ierr )
276 lwork_trevc = int( work(1) )
277 maxwrk = max( maxwrk, n + lwork_trevc )
278 CALL zhseqr(
'S',
'V', n, 1, n, a, lda, w, vl, ldvl,
279 $ work, -1, info )
280 ELSE IF( wantvr ) THEN
281 maxwrk = max( maxwrk, n + ( n - 1 )*
ilaenv( 1,
282 $ 'ZUNGHR',
283 $ ' ', n, 1, n, -1 ) )
284 CALL ztrevc3(
'R',
'B',
SELECT, n, a, lda,
285 $ vl, ldvl, vr, ldvr,
286 $ n, nout, work, -1, rwork, -1, ierr )
287 lwork_trevc = int( work(1) )
288 maxwrk = max( maxwrk, n + lwork_trevc )
289 CALL zhseqr(
'S',
'V', n, 1, n, a, lda, w, vr, ldvr,
290 $ work, -1, info )
291 ELSE
292 CALL zhseqr(
'E',
'N', n, 1, n, a, lda, w, vr, ldvr,
293 $ work, -1, info )
294 END IF
295 hswork = int( work(1) )
296 maxwrk = max( maxwrk, hswork, minwrk )
297 END IF
298 work( 1 ) = maxwrk
299
300 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
301 info = -12
302 END IF
303 END IF
304
305 IF( info.NE.0 ) THEN
306 CALL xerbla(
'ZGEEV ', -info )
307 RETURN
308 ELSE IF( lquery ) THEN
309 RETURN
310 END IF
311
312
313
314 IF( n.EQ.0 )
315 $ RETURN
316
317
318
321 bignum = one / smlnum
322 smlnum = sqrt( smlnum ) / eps
323 bignum = one / smlnum
324
325
326
327 anrm =
zlange(
'M', n, n, a, lda, dum )
328 scalea = .false.
329 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
330 scalea = .true.
331 cscale = smlnum
332 ELSE IF( anrm.GT.bignum ) THEN
333 scalea = .true.
334 cscale = bignum
335 END IF
336 IF( scalea )
337 $
CALL zlascl(
'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
338
339
340
341
342
343 ibal = 1
344 CALL zgebal(
'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
345
346
347
348
349
350 itau = 1
351 iwrk = itau + n
352 CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
353 $ lwork-iwrk+1, ierr )
354
355 IF( wantvl ) THEN
356
357
358
359
360 side = 'L'
361 CALL zlacpy(
'L', n, n, a, lda, vl, ldvl )
362
363
364
365
366
367 CALL zunghr( n, ilo, ihi, vl, ldvl, work( itau ),
368 $ work( iwrk ),
369 $ lwork-iwrk+1, ierr )
370
371
372
373
374
375 iwrk = itau
376 CALL zhseqr(
'S',
'V', n, ilo, ihi, a, lda, w, vl, ldvl,
377 $ work( iwrk ), lwork-iwrk+1, info )
378
379 IF( wantvr ) THEN
380
381
382
383
384 side = 'B'
385 CALL zlacpy(
'F', n, n, vl, ldvl, vr, ldvr )
386 END IF
387
388 ELSE IF( wantvr ) THEN
389
390
391
392
393 side = 'R'
394 CALL zlacpy(
'L', n, n, a, lda, vr, ldvr )
395
396
397
398
399
400 CALL zunghr( n, ilo, ihi, vr, ldvr, work( itau ),
401 $ work( iwrk ),
402 $ lwork-iwrk+1, ierr )
403
404
405
406
407
408 iwrk = itau
409 CALL zhseqr(
'S',
'V', n, ilo, ihi, a, lda, w, vr, ldvr,
410 $ work( iwrk ), lwork-iwrk+1, info )
411
412 ELSE
413
414
415
416
417
418 iwrk = itau
419 CALL zhseqr(
'E',
'N', n, ilo, ihi, a, lda, w, vr, ldvr,
420 $ work( iwrk ), lwork-iwrk+1, info )
421 END IF
422
423
424
425 IF( info.NE.0 )
426 $ GO TO 50
427
428 IF( wantvl .OR. wantvr ) THEN
429
430
431
432
433
434 irwork = ibal + n
435 CALL ztrevc3( side,
'B',
SELECT, n, a, lda, vl, ldvl, vr,
436 $ ldvr,
437 $ n, nout, work( iwrk ), lwork-iwrk+1,
438 $ rwork( irwork ), n, ierr )
439 END IF
440
441 IF( wantvl ) THEN
442
443
444
445
446
447 CALL zgebak(
'B',
'L', n, ilo, ihi, rwork( ibal ), n, vl,
448 $ ldvl,
449 $ ierr )
450
451
452
453 DO 20 i = 1, n
454 scl = one /
dznrm2( n, vl( 1, i ), 1 )
455 CALL zdscal( n, scl, vl( 1, i ), 1 )
456 DO 10 k = 1, n
457 rwork( irwork+k-1 ) = dble( vl( k, i ) )**2 +
458 $ aimag( vl( k, i ) )**2
459 10 CONTINUE
460 k =
idamax( n, rwork( irwork ), 1 )
461 tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
462 CALL zscal( n, tmp, vl( 1, i ), 1 )
463 vl( k, i ) = dcmplx( dble( vl( k, i ) ), zero )
464 20 CONTINUE
465 END IF
466
467 IF( wantvr ) THEN
468
469
470
471
472
473 CALL zgebak(
'B',
'R', n, ilo, ihi, rwork( ibal ), n, vr,
474 $ ldvr,
475 $ ierr )
476
477
478
479 DO 40 i = 1, n
480 scl = one /
dznrm2( n, vr( 1, i ), 1 )
481 CALL zdscal( n, scl, vr( 1, i ), 1 )
482 DO 30 k = 1, n
483 rwork( irwork+k-1 ) = dble( vr( k, i ) )**2 +
484 $ aimag( vr( k, i ) )**2
485 30 CONTINUE
486 k =
idamax( n, rwork( irwork ), 1 )
487 tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
488 CALL zscal( n, tmp, vr( 1, i ), 1 )
489 vr( k, i ) = dcmplx( dble( vr( k, i ) ), zero )
490 40 CONTINUE
491 END IF
492
493
494
495 50 CONTINUE
496 IF( scalea ) THEN
497 CALL zlascl(
'G', 0, 0, cscale, anrm, n-info, 1,
498 $ w( info+1 ),
499 $ max( n-info, 1 ), ierr )
500 IF( info.GT.0 ) THEN
501 CALL zlascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, w, n,
502 $ ierr )
503 END IF
504 END IF
505
506 work( 1 ) = maxwrk
507 RETURN
508
509
510
subroutine xerbla(srname, info)
subroutine zgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
ZGEBAK
subroutine zgebal(job, n, a, lda, ilo, ihi, scale, info)
ZGEBAL
subroutine zgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZGEHRD
subroutine zhseqr(job, compz, n, ilo, ihi, h, ldh, w, z, ldz, work, lwork, info)
ZHSEQR
integer function idamax(n, dx, incx)
IDAMAX
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
double precision function zlange(norm, m, n, a, lda, work)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine zlascl(type, kl, ku, cfrom, cto, m, n, a, lda, info)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
logical function lsame(ca, cb)
LSAME
real(wp) function dznrm2(n, x, incx)
DZNRM2
subroutine zdscal(n, da, zx, incx)
ZDSCAL
subroutine zscal(n, za, zx, incx)
ZSCAL
subroutine ztrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, rwork, lrwork, info)
ZTREVC3
subroutine zunghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
ZUNGHR