293 SUBROUTINE zdrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
294 $ ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
295 $ S, DTRU, DIF, DIFTRU, WORK, LWORK, RWORK,
296 $ IWORK, LIWORK, RESULT, BWORK, INFO )
303 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
305 DOUBLE PRECISION THRESH
310 DOUBLE PRECISION DIF( * ), DIFTRU( * ), DTRU( * ), LSCALE( * ),
311 $ result( 4 ), rscale( * ), rwork( * ), s( * )
312 COMPLEX*16 A( LDA, * ), AI( LDA, * ), ALPHA( * ),
313 $ b( lda, * ), beta( * ), bi( lda, * ),
314 $ vl( lda, * ), vr( lda, * ), work( * )
320 DOUBLE PRECISION ZERO, ONE, TEN, TNTH, HALF
321 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
322 $ tnth = 1.0d-1, half = 0.5d+0 )
325 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
326 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
327 DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
331 COMPLEX*16 WEIGHT( 5 )
335 DOUBLE PRECISION DLAMCH, ZLANGE
336 EXTERNAL ILAENV, DLAMCH, ZLANGE
342 INTRINSIC abs, dcmplx, max, sqrt
352 IF( nsize.LT.0 )
THEN
354 ELSE IF( thresh.LT.zero )
THEN
356 ELSE IF( nin.LE.0 )
THEN
358 ELSE IF( nout.LE.0 )
THEN
360 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
362 ELSE IF( liwork.LT.nmax+2 )
THEN
374 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
375 minwrk = 2*nmax*( nmax+1 )
376 maxwrk = nmax*( 1+ilaenv( 1,
'ZGEQRF',
' ', nmax, 1, nmax,
378 maxwrk = max( maxwrk, 2*nmax*( nmax+1 ) )
382 IF( lwork.LT.minwrk )
386 CALL xerbla(
'ZDRGVX', -info )
403 weight( 1 ) = dcmplx( tnth, zero )
404 weight( 2 ) = dcmplx( half, zero )
406 weight( 4 ) = one / weight( 2 )
407 weight( 5 ) = one / weight( 1 )
417 CALL zlatm6( iptype, 5, a, lda, b, vr, lda, vl,
418 $ lda, weight( iwa ), weight( iwb ),
419 $ weight( iwx ), weight( iwy ), dtru,
426 CALL zlacpy(
'F', n, n, a, lda, ai, lda )
427 CALL zlacpy(
'F', n, n, b, lda, bi, lda )
429 CALL zggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
430 $ lda, alpha, beta, vl, lda, vr, lda,
431 $ ilo, ihi, lscale, rscale, anorm,
432 $ bnorm, s, dif, work, lwork, rwork,
433 $ iwork, bwork, linfo )
434 IF( linfo.NE.0 )
THEN
435 WRITE( nout, fmt = 9999 )
'ZGGEVX', linfo, n,
436 $ iptype, iwa, iwb, iwx, iwy
442 CALL zlacpy(
'Full', n, n, ai, lda, work, n )
443 CALL zlacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
445 abnorm = zlange(
'Fro', n, 2*n, work, n, rwork )
450 CALL zget52( .true., n, a, lda, b, lda, vl, lda,
451 $ alpha, beta, work, rwork,
453 IF( result( 2 ).GT.thresh )
THEN
454 WRITE( nout, fmt = 9998 )
'Left',
'ZGGEVX',
455 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
459 CALL zget52( .false., n, a, lda, b, lda, vr, lda,
460 $ alpha, beta, work, rwork,
462 IF( result( 3 ).GT.thresh )
THEN
463 WRITE( nout, fmt = 9998 )
'Right',
'ZGGEVX',
464 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
471 IF( s( i ).EQ.zero )
THEN
472 IF( dtru( i ).GT.abnorm*ulp )
473 $ result( 3 ) = ulpinv
474 ELSE IF( dtru( i ).EQ.zero )
THEN
475 IF( s( i ).GT.abnorm*ulp )
476 $ result( 3 ) = ulpinv
478 rwork( i ) = max( abs( dtru( i ) / s( i ) ),
479 $ abs( s( i ) / dtru( i ) ) )
480 result( 3 ) = max( result( 3 ), rwork( i ) )
487 IF( dif( 1 ).EQ.zero )
THEN
488 IF( diftru( 1 ).GT.abnorm*ulp )
489 $ result( 4 ) = ulpinv
490 ELSE IF( diftru( 1 ).EQ.zero )
THEN
491 IF( dif( 1 ).GT.abnorm*ulp )
492 $ result( 4 ) = ulpinv
493 ELSE IF( dif( 5 ).EQ.zero )
THEN
494 IF( diftru( 5 ).GT.abnorm*ulp )
495 $ result( 4 ) = ulpinv
496 ELSE IF( diftru( 5 ).EQ.zero )
THEN
497 IF( dif( 5 ).GT.abnorm*ulp )
498 $ result( 4 ) = ulpinv
500 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
501 $ abs( dif( 1 ) / diftru( 1 ) ) )
502 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
503 $ abs( dif( 5 ) / diftru( 5 ) ) )
504 result( 4 ) = max( ratio1, ratio2 )
512 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
513 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
519 IF( nerrs.EQ.0 )
THEN
520 WRITE( nout, fmt = 9997 )
'ZXV'
526 WRITE( nout, fmt = 9995 )
527 WRITE( nout, fmt = 9994 )
528 WRITE( nout, fmt = 9993 )
532 WRITE( nout, fmt = 9992 )
'''',
537 IF( result( j ).LT.10000.0d0 )
THEN
538 WRITE( nout, fmt = 9991 )iptype, iwa,
539 $ iwb, iwx, iwy, j, result( j )
541 WRITE( nout, fmt = 9990 )iptype, iwa,
542 $ iwb, iwx, iwy, j, result( j )
562 READ( nin, fmt = *,
END = 150 )n
566 READ( nin, fmt = * )( a( i, j ), j = 1, n )
569 READ( nin, fmt = * )( b( i, j ), j = 1, n )
571 READ( nin, fmt = * )( dtru( i ), i = 1, n )
572 READ( nin, fmt = * )( diftru( i ), i = 1, n )
580 CALL zlacpy(
'F', n, n, a, lda, ai, lda )
581 CALL zlacpy(
'F', n, n, b, lda, bi, lda )
583 CALL zggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alpha, beta,
584 $ vl, lda, vr, lda, ilo, ihi, lscale, rscale, anorm,
585 $ bnorm, s, dif, work, lwork, rwork, iwork, bwork,
588 IF( linfo.NE.0 )
THEN
589 WRITE( nout, fmt = 9987 )
'ZGGEVX', linfo, n, nptknt
595 CALL zlacpy(
'Full', n, n, ai, lda, work, n )
596 CALL zlacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
597 abnorm =
zlange(
'Fro', n, 2*n, work, n, rwork )
602 CALL zget52( .true., n, a, lda, b, lda, vl, lda, alpha, beta,
603 $ work, rwork, result( 1 ) )
604 IF( result( 2 ).GT.thresh )
THEN
605 WRITE( nout, fmt = 9986 )
'Left',
'ZGGEVX', result( 2 ), n,
610 CALL zget52( .false., n, a, lda, b, lda, vr, lda, alpha, beta,
611 $ work, rwork, result( 2 ) )
612 IF( result( 3 ).GT.thresh )
THEN
613 WRITE( nout, fmt = 9986 )
'Right',
'ZGGEVX', result( 3 ), n,
621 IF( s( i ).EQ.zero )
THEN
622 IF( dtru( i ).GT.abnorm*ulp )
623 $ result( 3 ) = ulpinv
624 ELSE IF( dtru( i ).EQ.zero )
THEN
625 IF( s( i ).GT.abnorm*ulp )
626 $ result( 3 ) = ulpinv
628 rwork( i ) = max( abs( dtru( i ) / s( i ) ),
629 $ abs( s( i ) / dtru( i ) ) )
630 result( 3 ) = max( result( 3 ), rwork( i ) )
637 IF( dif( 1 ).EQ.zero )
THEN
638 IF( diftru( 1 ).GT.abnorm*ulp )
639 $ result( 4 ) = ulpinv
640 ELSE IF( diftru( 1 ).EQ.zero )
THEN
641 IF( dif( 1 ).GT.abnorm*ulp )
642 $ result( 4 ) = ulpinv
643 ELSE IF( dif( 5 ).EQ.zero )
THEN
644 IF( diftru( 5 ).GT.abnorm*ulp )
645 $ result( 4 ) = ulpinv
646 ELSE IF( diftru( 5 ).EQ.zero )
THEN
647 IF( dif( 5 ).GT.abnorm*ulp )
648 $ result( 4 ) = ulpinv
650 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
651 $ abs( dif( 1 ) / diftru( 1 ) ) )
652 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
653 $ abs( dif( 5 ) / diftru( 5 ) ) )
654 result( 4 ) = max( ratio1, ratio2 )
662 IF( result( j ).GE.thrsh2 )
THEN
667 IF( nerrs.EQ.0 )
THEN
668 WRITE( nout, fmt = 9997 )
'ZXV'
674 WRITE( nout, fmt = 9996 )
678 WRITE( nout, fmt = 9992 )
'''',
'transpose',
''''
682 IF( result( j ).LT.10000.0d0 )
THEN
683 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
685 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
697 CALL alasvm(
'ZXV', nout, nerrs, ntestt, 0 )
703 9999
FORMAT(
' ZDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
704 $ i6,
', JTYPE=', i6,
')' )
706 9998
FORMAT(
' ZDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
707 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
708 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
709 $
', IWX=', i5,
', IWY=', i5 )
711 9997
FORMAT( / 1x, a3,
' -- Complex Expert Eigenvalue/vector',
712 $
' problem driver' )
714 9996
FORMAT(
'Input Example' )
716 9995
FORMAT(
' Matrix types: ', / )
718 9994
FORMAT(
' TYPE 1: Da is diagonal, Db is identity, ',
719 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
720 $ /
' YH and X are left and right eigenvectors. ', / )
722 9993
FORMAT(
' TYPE 2: Da is quasi-diagonal, Db is identity, ',
723 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
724 $ /
' YH and X are left and right eigenvectors. ', / )
726 9992
FORMAT( /
' Tests performed: ', / 4x,
727 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
728 $
' r is a right eigenvector and ', a,
' means ', a,
'.',
729 $ /
' 1 = max | ( b A - a B )', a,
' l | / const.',
730 $ /
' 2 = max | ( b A - a B ) r | / const.',
731 $ /
' 3 = max ( Sest/Stru, Stru/Sest ) ',
732 $
' over all eigenvalues', /
733 $
' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
734 $
' over the 1st and 5th eigenvectors', / )
736 9991
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
737 $ i2,
', IWY=', i2,
', result ', i2,
' is', 0p, f8.2 )
739 9990
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
740 $ i2,
', IWY=', i2,
', result ', i2,
' is', 1p, d10.3 )
742 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
743 $
' result ', i2,
' is', 0p, f8.2 )
745 9988
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
746 $
' result ', i2,
' is', 1p, d10.3 )
748 9987
FORMAT(
' ZDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
749 $ i6,
', Input example #', i2,
')' )
751 9986
FORMAT(
' ZDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
752 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
753 $
'N=', i6,
', Input Example #', i2,
')' )