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 CALL dlabad( smlnum, bignum )
319 smlnum = sqrt( smlnum ) / eps
320 bignum = one / smlnum
321
322
323
324 anrm =
zlange(
'M', n, n, a, lda, dum )
325 scalea = .false.
326 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
327 scalea = .true.
328 cscale = smlnum
329 ELSE IF( anrm.GT.bignum ) THEN
330 scalea = .true.
331 cscale = bignum
332 END IF
333 IF( scalea )
334 $
CALL zlascl(
'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
335
336
337
338
339
340 ibal = 1
341 CALL zgebal(
'B', n, a, lda, ilo, ihi, rwork( ibal ), ierr )
342
343
344
345
346
347 itau = 1
348 iwrk = itau + n
349 CALL zgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
350 $ lwork-iwrk+1, ierr )
351
352 IF( wantvl ) THEN
353
354
355
356
357 side = 'L'
358 CALL zlacpy(
'L', n, n, a, lda, vl, ldvl )
359
360
361
362
363
364 CALL zunghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
365 $ lwork-iwrk+1, ierr )
366
367
368
369
370
371 iwrk = itau
372 CALL zhseqr(
'S',
'V', n, ilo, ihi, a, lda, w, vl, ldvl,
373 $ work( iwrk ), lwork-iwrk+1, info )
374
375 IF( wantvr ) THEN
376
377
378
379
380 side = 'B'
381 CALL zlacpy(
'F', n, n, vl, ldvl, vr, ldvr )
382 END IF
383
384 ELSE IF( wantvr ) THEN
385
386
387
388
389 side = 'R'
390 CALL zlacpy(
'L', n, n, a, lda, vr, ldvr )
391
392
393
394
395
396 CALL zunghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
397 $ lwork-iwrk+1, ierr )
398
399
400
401
402
403 iwrk = itau
404 CALL zhseqr(
'S',
'V', n, ilo, ihi, a, lda, w, vr, ldvr,
405 $ work( iwrk ), lwork-iwrk+1, info )
406
407 ELSE
408
409
410
411
412
413 iwrk = itau
414 CALL zhseqr(
'E',
'N', n, ilo, ihi, a, lda, w, vr, ldvr,
415 $ work( iwrk ), lwork-iwrk+1, info )
416 END IF
417
418
419
420 IF( info.NE.0 )
421 $ GO TO 50
422
423 IF( wantvl .OR. wantvr ) THEN
424
425
426
427
428
429 irwork = ibal + n
430 CALL ztrevc3( side,
'B',
SELECT, n, a, lda, vl, ldvl, vr, ldvr,
431 $ n, nout, work( iwrk ), lwork-iwrk+1,
432 $ rwork( irwork ), n, ierr )
433 END IF
434
435 IF( wantvl ) THEN
436
437
438
439
440
441 CALL zgebak(
'B',
'L', n, ilo, ihi, rwork( ibal ), n, vl, ldvl,
442 $ ierr )
443
444
445
446 DO 20 i = 1, n
447 scl = one /
dznrm2( n, vl( 1, i ), 1 )
448 CALL zdscal( n, scl, vl( 1, i ), 1 )
449 DO 10 k = 1, n
450 rwork( irwork+k-1 ) = dble( vl( k, i ) )**2 +
451 $ aimag( vl( k, i ) )**2
452 10 CONTINUE
453 k =
idamax( n, rwork( irwork ), 1 )
454 tmp = conjg( vl( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
455 CALL zscal( n, tmp, vl( 1, i ), 1 )
456 vl( k, i ) = dcmplx( dble( vl( k, i ) ), zero )
457 20 CONTINUE
458 END IF
459
460 IF( wantvr ) THEN
461
462
463
464
465
466 CALL zgebak(
'B',
'R', n, ilo, ihi, rwork( ibal ), n, vr, ldvr,
467 $ ierr )
468
469
470
471 DO 40 i = 1, n
472 scl = one /
dznrm2( n, vr( 1, i ), 1 )
473 CALL zdscal( n, scl, vr( 1, i ), 1 )
474 DO 30 k = 1, n
475 rwork( irwork+k-1 ) = dble( vr( k, i ) )**2 +
476 $ aimag( vr( k, i ) )**2
477 30 CONTINUE
478 k =
idamax( n, rwork( irwork ), 1 )
479 tmp = conjg( vr( k, i ) ) / sqrt( rwork( irwork+k-1 ) )
480 CALL zscal( n, tmp, vr( 1, i ), 1 )
481 vr( k, i ) = dcmplx( dble( vr( k, i ) ), zero )
482 40 CONTINUE
483 END IF
484
485
486
487 50 CONTINUE
488 IF( scalea ) THEN
489 CALL zlascl(
'G', 0, 0, cscale, anrm, n-info, 1, w( info+1 ),
490 $ max( n-info, 1 ), ierr )
491 IF( info.GT.0 ) THEN
492 CALL zlascl(
'G', 0, 0, cscale, anrm, ilo-1, 1, w, n, ierr )
493 END IF
494 END IF
495
496 work( 1 ) = maxwrk
497 RETURN
498
499
500
double precision function dlamch(CMACH)
DLAMCH
subroutine dlabad(SMALL, LARGE)
DLABAD
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
integer function idamax(N, DX, INCX)
IDAMAX
subroutine xerbla(SRNAME, INFO)
XERBLA
logical function lsame(CA, CB)
LSAME
subroutine zdscal(N, DA, ZX, INCX)
ZDSCAL
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
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 zgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
ZGEBAL
subroutine zgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
ZGEHRD
subroutine zgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
ZGEBAK
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.
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, W, Z, LDZ, WORK, LWORK, INFO)
ZHSEQR
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
real(wp) function dznrm2(n, x, incx)
DZNRM2