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,
')' )
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine xerbla(srname, info)
subroutine cdrgvx(nsize, thresh, nin, nout, a, lda, b, ai, bi, alpha, beta, vl, vr, ilo, ihi, lscale, rscale, s, stru, dif, diftru, work, lwork, rwork, iwork, liwork, result, bwork, info)
CDRGVX
subroutine cget52(left, n, a, lda, b, ldb, e, lde, alpha, beta, work, rwork, result)
CGET52
subroutine clatm6(type, n, a, lda, b, x, ldx, y, ldy, alpha, beta, wx, wy, s, dif)
CLATM6
subroutine cggevx(balanc, jobvl, jobvr, sense, n, a, lda, b, ldb, alpha, beta, vl, ldvl, vr, ldvr, ilo, ihi, lscale, rscale, abnrm, bbnrm, rconde, rcondv, work, lwork, rwork, iwork, bwork, info)
CGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...