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