294 SUBROUTINE cdrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
295 $ ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
296 $ S, STRU, DIF, DIFTRU, WORK, LWORK, RWORK,
297 $ IWORK, LIWORK, RESULT, BWORK, INFO )
304 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
311 REAL DIF( * ), DIFTRU( * ), LSCALE( * ),
312 $ result( 4 ), rscale( * ), rwork( * ), s( * ),
314 COMPLEX A( LDA, * ), AI( LDA, * ), ALPHA( * ),
315 $ B( LDA, * ), BETA( * ), BI( LDA, * ),
316 $ VL( LDA, * ), VR( LDA, * ), WORK( * )
322 REAL ZERO, ONE, TEN, TNTH, HALF
323 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
324 $ tnth = 1.0e-1, half = 0.5e+0 )
327 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
328 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
329 REAL ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
338 EXTERNAL ILAENV, CLANGE, SLAMCH
344 INTRINSIC abs, cmplx, max, sqrt
354 IF( nsize.LT.0 )
THEN
356 ELSE IF( thresh.LT.zero )
THEN
358 ELSE IF( nin.LE.0 )
THEN
360 ELSE IF( nout.LE.0 )
THEN
362 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
364 ELSE IF( liwork.LT.nmax+2 )
THEN
376 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
377 minwrk = 2*nmax*( nmax+1 )
378 maxwrk = nmax*( 1+ilaenv( 1,
'CGEQRF',
' ', nmax, 1, nmax,
380 maxwrk = max( maxwrk, 2*nmax*( nmax+1 ) )
384 IF( lwork.LT.minwrk )
388 CALL xerbla(
'CDRGVX', -info )
405 weight( 1 ) = cmplx( tnth, zero )
406 weight( 2 ) = cmplx( half, zero )
408 weight( 4 ) = one / weight( 2 )
409 weight( 5 ) = one / weight( 1 )
419 CALL clatm6( iptype, 5, a, lda, b, vr, lda, vl,
420 $ lda, weight( iwa ), weight( iwb ),
421 $ weight( iwx ), weight( iwy ), stru,
428 CALL clacpy(
'F', n, n, a, lda, ai, lda )
429 CALL clacpy(
'F', n, n, b, lda, bi, lda )
431 CALL cggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
432 $ lda, alpha, beta, vl, lda, vr, lda,
433 $ ilo, ihi, lscale, rscale, anorm,
434 $ bnorm, s, dif, work, lwork, rwork,
435 $ iwork, bwork, linfo )
436 IF( linfo.NE.0 )
THEN
437 WRITE( nout, fmt = 9999 )
'CGGEVX', linfo, n,
438 $ iptype, iwa, iwb, iwx, iwy
444 CALL clacpy(
'Full', n, n, ai, lda, work, n )
445 CALL clacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
447 abnorm = clange(
'Fro', n, 2*n, work, n, rwork )
452 CALL cget52( .true., n, a, lda, b, lda, vl, lda,
453 $ alpha, beta, work, rwork,
455 IF( result( 2 ).GT.thresh )
THEN
456 WRITE( nout, fmt = 9998 )
'Left',
'CGGEVX',
457 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
461 CALL cget52( .false., n, a, lda, b, lda, vr, lda,
462 $ alpha, beta, work, rwork,
464 IF( result( 3 ).GT.thresh )
THEN
465 WRITE( nout, fmt = 9998 )
'Right',
'CGGEVX',
466 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
473 IF( s( i ).EQ.zero )
THEN
474 IF( stru( i ).GT.abnorm*ulp )
475 $ result( 3 ) = ulpinv
476 ELSE IF( stru( i ).EQ.zero )
THEN
477 IF( s( i ).GT.abnorm*ulp )
478 $ result( 3 ) = ulpinv
480 rwork( i ) = max( abs( stru( i ) / s( i ) ),
481 $ abs( s( i ) / stru( i ) ) )
482 result( 3 ) = max( result( 3 ), rwork( i ) )
489 IF( dif( 1 ).EQ.zero )
THEN
490 IF( diftru( 1 ).GT.abnorm*ulp )
491 $ result( 4 ) = ulpinv
492 ELSE IF( diftru( 1 ).EQ.zero )
THEN
493 IF( dif( 1 ).GT.abnorm*ulp )
494 $ result( 4 ) = ulpinv
495 ELSE IF( dif( 5 ).EQ.zero )
THEN
496 IF( diftru( 5 ).GT.abnorm*ulp )
497 $ result( 4 ) = ulpinv
498 ELSE IF( diftru( 5 ).EQ.zero )
THEN
499 IF( dif( 5 ).GT.abnorm*ulp )
500 $ result( 4 ) = ulpinv
502 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
503 $ abs( dif( 1 ) / diftru( 1 ) ) )
504 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
505 $ abs( dif( 5 ) / diftru( 5 ) ) )
506 result( 4 ) = max( ratio1, ratio2 )
514 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
515 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
521 IF( nerrs.EQ.0 )
THEN
522 WRITE( nout, fmt = 9997 )
'CXV'
528 WRITE( nout, fmt = 9995 )
529 WRITE( nout, fmt = 9994 )
530 WRITE( nout, fmt = 9993 )
534 WRITE( nout, fmt = 9992 )
'''',
539 IF( result( j ).LT.10000.0 )
THEN
540 WRITE( nout, fmt = 9991 )iptype, iwa,
541 $ iwb, iwx, iwy, j, result( j )
543 WRITE( nout, fmt = 9990 )iptype, iwa,
544 $ iwb, iwx, iwy, j, result( j )
564 READ( nin, fmt = *,
END = 150 )n
568 READ( nin, fmt = * )( a( i, j ), j = 1, n )
571 READ( nin, fmt = * )( b( i, j ), j = 1, n )
573 READ( nin, fmt = * )( stru( i ), i = 1, n )
574 READ( nin, fmt = * )( diftru( i ), i = 1, n )
582 CALL clacpy(
'F', n, n, a, lda, ai, lda )
583 CALL clacpy(
'F', n, n, b, lda, bi, lda )
585 CALL cggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alpha, beta,
586 $ vl, lda, vr, lda, ilo, ihi, lscale, rscale, anorm,
587 $ bnorm, s, dif, work, lwork, rwork, iwork, bwork,
590 IF( linfo.NE.0 )
THEN
591 WRITE( nout, fmt = 9987 )
'CGGEVX', linfo, n, nptknt
597 CALL clacpy(
'Full', n, n, ai, lda, work, n )
598 CALL clacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
599 abnorm =
clange(
'Fro', n, 2*n, work, n, rwork )
604 CALL cget52( .true., n, a, lda, b, lda, vl, lda, alpha, beta,
605 $ work, rwork, result( 1 ) )
606 IF( result( 2 ).GT.thresh )
THEN
607 WRITE( nout, fmt = 9986 )
'Left',
'CGGEVX', result( 2 ), n,
612 CALL cget52( .false., n, a, lda, b, lda, vr, lda, alpha, beta,
613 $ work, rwork, result( 2 ) )
614 IF( result( 3 ).GT.thresh )
THEN
615 WRITE( nout, fmt = 9986 )
'Right',
'CGGEVX', result( 3 ), n,
623 IF( s( i ).EQ.zero )
THEN
624 IF( stru( i ).GT.abnorm*ulp )
625 $ result( 3 ) = ulpinv
626 ELSE IF( stru( i ).EQ.zero )
THEN
627 IF( s( i ).GT.abnorm*ulp )
628 $ result( 3 ) = ulpinv
630 rwork( i ) = max( abs( stru( i ) / s( i ) ),
631 $ abs( s( i ) / stru( i ) ) )
632 result( 3 ) = max( result( 3 ), rwork( i ) )
639 IF( dif( 1 ).EQ.zero )
THEN
640 IF( diftru( 1 ).GT.abnorm*ulp )
641 $ result( 4 ) = ulpinv
642 ELSE IF( diftru( 1 ).EQ.zero )
THEN
643 IF( dif( 1 ).GT.abnorm*ulp )
644 $ result( 4 ) = ulpinv
645 ELSE IF( dif( 5 ).EQ.zero )
THEN
646 IF( diftru( 5 ).GT.abnorm*ulp )
647 $ result( 4 ) = ulpinv
648 ELSE IF( diftru( 5 ).EQ.zero )
THEN
649 IF( dif( 5 ).GT.abnorm*ulp )
650 $ result( 4 ) = ulpinv
652 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
653 $ abs( dif( 1 ) / diftru( 1 ) ) )
654 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
655 $ abs( dif( 5 ) / diftru( 5 ) ) )
656 result( 4 ) = max( ratio1, ratio2 )
664 IF( result( j ).GE.thrsh2 )
THEN
669 IF( nerrs.EQ.0 )
THEN
670 WRITE( nout, fmt = 9997 )
'CXV'
676 WRITE( nout, fmt = 9996 )
680 WRITE( nout, fmt = 9992 )
'''',
'transpose',
''''
684 IF( result( j ).LT.10000.0 )
THEN
685 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
687 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
699 CALL alasvm(
'CXV', nout, nerrs, ntestt, 0 )
705 9999
FORMAT(
' CDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
706 $ i6,
', JTYPE=', i6,
')' )
708 9998
FORMAT(
' CDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
709 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
710 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
711 $
', IWX=', i5,
', IWY=', i5 )
713 9997
FORMAT( / 1x, a3,
' -- Complex Expert Eigenvalue/vector',
714 $
' problem driver' )
716 9996
FORMAT(
'Input Example' )
718 9995
FORMAT(
' Matrix types: ', / )
720 9994
FORMAT(
' TYPE 1: Da is diagonal, Db is identity, ',
721 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
722 $ /
' YH and X are left and right eigenvectors. ', / )
724 9993
FORMAT(
' TYPE 2: Da is quasi-diagonal, Db is identity, ',
725 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
726 $ /
' YH and X are left and right eigenvectors. ', / )
728 9992
FORMAT( /
' Tests performed: ', / 4x,
729 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
730 $
' r is a right eigenvector and ', a,
' means ', a,
'.',
731 $ /
' 1 = max | ( b A - a B )', a,
' l | / const.',
732 $ /
' 2 = max | ( b A - a B ) r | / const.',
733 $ /
' 3 = max ( Sest/Stru, Stru/Sest ) ',
734 $
' over all eigenvalues', /
735 $
' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
736 $
' over the 1st and 5th eigenvectors', / )
738 9991
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
739 $ i2,
', IWY=', i2,
', result ', i2,
' is', 0p, f8.2 )
741 9990
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
742 $ i2,
', IWY=', i2,
', result ', i2,
' is', 1p, e10.3 )
744 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
745 $
' result ', i2,
' is', 0p, f8.2 )
747 9988
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
748 $
' result ', i2,
' is', 1p, e10.3 )
750 9987
FORMAT(
' CDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
751 $ i6,
', Input example #', i2,
')' )
753 9986
FORMAT(
' CDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
754 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
755 $
'N=', i6,
', Input Example #', i2,
')' )