298 SUBROUTINE ddrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
299 $ alphar, alphai, beta, vl, vr, ilo, ihi, lscale,
300 $ rscale, s, dtru, dif, diftru, work, lwork,
301 $ iwork, liwork, result, bwork, info )
309 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
311 DOUBLE PRECISION THRESH
316 DOUBLE PRECISION A( lda, * ), AI( lda, * ), ALPHAI( * ),
317 $ alphar( * ), b( lda, * ), beta( * ),
318 $ bi( lda, * ), dif( * ), diftru( * ), dtru( * ),
319 $ lscale( * ), result( 4 ), rscale( * ), s( * ),
320 $ vl( lda, * ), vr( lda, * ), work( * )
326 DOUBLE PRECISION ZERO, ONE, TEN, TNTH, HALF
327 parameter ( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
328 $ tnth = 1.0d-1, half = 0.5d+0 )
331 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
332 $ maxwrk, minwrk, n, nerrs, nmax, nptknt, ntestt
333 DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
337 DOUBLE PRECISION WEIGHT( 5 )
341 DOUBLE PRECISION DLAMCH, DLANGE
342 EXTERNAL ilaenv, dlamch, dlange
348 INTRINSIC abs, max, sqrt
358 IF( nsize.LT.0 )
THEN
360 ELSE IF( thresh.LT.zero )
THEN
362 ELSE IF( nin.LE.0 )
THEN
364 ELSE IF( nout.LE.0 )
THEN
366 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
368 ELSE IF( liwork.LT.nmax+6 )
THEN
380 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
381 minwrk = 2*nmax*nmax + 12*nmax + 16
382 maxwrk = 6*nmax + nmax*ilaenv( 1,
'DGEQRF',
' ', nmax, 1, nmax,
384 maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
388 IF( lwork.LT.minwrk )
392 CALL xerbla(
'DDRGVX', -info )
412 weight( 4 ) = one / weight( 2 )
413 weight( 5 ) = one / weight( 1 )
423 CALL dlatm6( iptype, 5, a, lda, b, vr, lda, vl,
424 $ lda, weight( iwa ), weight( iwb ),
425 $ weight( iwx ), weight( iwy ), dtru,
432 CALL dlacpy(
'F', n, n, a, lda, ai, lda )
433 CALL dlacpy(
'F', n, n, b, lda, bi, lda )
435 CALL dggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
436 $ lda, alphar, alphai, beta, vl, lda,
437 $ vr, lda, ilo, ihi, lscale, rscale,
438 $ anorm, bnorm, s, dif, work, lwork,
439 $ iwork, bwork, linfo )
440 IF( linfo.NE.0 )
THEN
442 WRITE( nout, fmt = 9999 )
'DGGEVX', linfo, n,
449 CALL dlacpy(
'Full', n, n, ai, lda, work, n )
450 CALL dlacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
452 abnorm = dlange(
'Fro', n, 2*n, work, n, work )
457 CALL dget52( .true., n, a, lda, b, lda, vl, lda,
458 $ alphar, alphai, beta, work,
460 IF( result( 2 ).GT.thresh )
THEN
461 WRITE( nout, fmt = 9998 )
'Left',
'DGGEVX',
462 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
466 CALL dget52( .false., n, a, lda, b, lda, vr, lda,
467 $ alphar, alphai, beta, work,
469 IF( result( 3 ).GT.thresh )
THEN
470 WRITE( nout, fmt = 9998 )
'Right',
'DGGEVX',
471 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
478 IF( s( i ).EQ.zero )
THEN
479 IF( dtru( i ).GT.abnorm*ulp )
480 $ result( 3 ) = ulpinv
481 ELSE IF( dtru( i ).EQ.zero )
THEN
482 IF( s( i ).GT.abnorm*ulp )
483 $ result( 3 ) = ulpinv
485 work( i ) = max( abs( dtru( i ) / s( i ) ),
486 $ abs( s( i ) / dtru( i ) ) )
487 result( 3 ) = max( result( 3 ), work( i ) )
494 IF( dif( 1 ).EQ.zero )
THEN
495 IF( diftru( 1 ).GT.abnorm*ulp )
496 $ result( 4 ) = ulpinv
497 ELSE IF( diftru( 1 ).EQ.zero )
THEN
498 IF( dif( 1 ).GT.abnorm*ulp )
499 $ result( 4 ) = ulpinv
500 ELSE IF( dif( 5 ).EQ.zero )
THEN
501 IF( diftru( 5 ).GT.abnorm*ulp )
502 $ result( 4 ) = ulpinv
503 ELSE IF( diftru( 5 ).EQ.zero )
THEN
504 IF( dif( 5 ).GT.abnorm*ulp )
505 $ result( 4 ) = ulpinv
507 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
508 $ abs( dif( 1 ) / diftru( 1 ) ) )
509 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
510 $ abs( dif( 5 ) / diftru( 5 ) ) )
511 result( 4 ) = max( ratio1, ratio2 )
519 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
520 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
526 IF( nerrs.EQ.0 )
THEN
527 WRITE( nout, fmt = 9997 )
'DXV'
533 WRITE( nout, fmt = 9995 )
534 WRITE( nout, fmt = 9994 )
535 WRITE( nout, fmt = 9993 )
539 WRITE( nout, fmt = 9992 )
'''',
544 IF( result( j ).LT.10000.0d0 )
THEN
545 WRITE( nout, fmt = 9991 )iptype, iwa,
546 $ iwb, iwx, iwy, j, result( j )
548 WRITE( nout, fmt = 9990 )iptype, iwa,
549 $ iwb, iwx, iwy, j, result( j )
569 READ( nin, fmt = *, end = 150 )n
573 READ( nin, fmt = * )( a( i, j ), j = 1, n )
576 READ( nin, fmt = * )( b( i, j ), j = 1, n )
578 READ( nin, fmt = * )( dtru( i ), i = 1, n )
579 READ( nin, fmt = * )( diftru( i ), i = 1, n )
587 CALL dlacpy(
'F', n, n, a, lda, ai, lda )
588 CALL dlacpy(
'F', n, n, b, lda, bi, lda )
590 CALL dggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alphar,
591 $ alphai, beta, vl, lda, vr, lda, ilo, ihi, lscale,
592 $ rscale, anorm, bnorm, s, dif, work, lwork, iwork,
595 IF( linfo.NE.0 )
THEN
597 WRITE( nout, fmt = 9987 )
'DGGEVX', linfo, n, nptknt
603 CALL dlacpy(
'Full', n, n, ai, lda, work, n )
604 CALL dlacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
605 abnorm = dlange(
'Fro', n, 2*n, work, n, work )
610 CALL dget52( .true., n, a, lda, b, lda, vl, lda, alphar, alphai,
611 $ beta, work, result( 1 ) )
612 IF( result( 2 ).GT.thresh )
THEN
613 WRITE( nout, fmt = 9986 )
'Left',
'DGGEVX', result( 2 ), n,
618 CALL dget52( .false., n, a, lda, b, lda, vr, lda, alphar, alphai,
619 $ beta, work, result( 2 ) )
620 IF( result( 3 ).GT.thresh )
THEN
621 WRITE( nout, fmt = 9986 )
'Right',
'DGGEVX', result( 3 ), n,
629 IF( s( i ).EQ.zero )
THEN
630 IF( dtru( i ).GT.abnorm*ulp )
631 $ result( 3 ) = ulpinv
632 ELSE IF( dtru( i ).EQ.zero )
THEN
633 IF( s( i ).GT.abnorm*ulp )
634 $ result( 3 ) = ulpinv
636 work( i ) = max( abs( dtru( i ) / s( i ) ),
637 $ abs( s( i ) / dtru( i ) ) )
638 result( 3 ) = max( result( 3 ), work( i ) )
645 IF( dif( 1 ).EQ.zero )
THEN
646 IF( diftru( 1 ).GT.abnorm*ulp )
647 $ result( 4 ) = ulpinv
648 ELSE IF( diftru( 1 ).EQ.zero )
THEN
649 IF( dif( 1 ).GT.abnorm*ulp )
650 $ result( 4 ) = ulpinv
651 ELSE IF( dif( 5 ).EQ.zero )
THEN
652 IF( diftru( 5 ).GT.abnorm*ulp )
653 $ result( 4 ) = ulpinv
654 ELSE IF( diftru( 5 ).EQ.zero )
THEN
655 IF( dif( 5 ).GT.abnorm*ulp )
656 $ result( 4 ) = ulpinv
658 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
659 $ abs( dif( 1 ) / diftru( 1 ) ) )
660 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
661 $ abs( dif( 5 ) / diftru( 5 ) ) )
662 result( 4 ) = max( ratio1, ratio2 )
670 IF( result( j ).GE.thrsh2 )
THEN
675 IF( nerrs.EQ.0 )
THEN
676 WRITE( nout, fmt = 9997 )
'DXV'
682 WRITE( nout, fmt = 9996 )
686 WRITE( nout, fmt = 9992 )
'''',
'transpose',
''''
690 IF( result( j ).LT.10000.0d0 )
THEN
691 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
693 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
705 CALL alasvm(
'DXV', nout, nerrs, ntestt, 0 )
711 9999
FORMAT(
' DDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
712 $ i6,
', JTYPE=', i6,
')' )
714 9998
FORMAT(
' DDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
715 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
716 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
717 $
', IWX=', i5,
', IWY=', i5 )
719 9997
FORMAT( / 1x, a3,
' -- Real Expert Eigenvalue/vector',
720 $
' problem driver' )
722 9996
FORMAT(
' Input Example' )
724 9995
FORMAT(
' Matrix types: ', / )
726 9994
FORMAT(
' TYPE 1: Da is diagonal, Db is identity, ',
727 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
728 $ /
' YH and X are left and right eigenvectors. ', / )
730 9993
FORMAT(
' TYPE 2: Da is quasi-diagonal, Db is identity, ',
731 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
732 $ /
' YH and X are left and right eigenvectors. ', / )
734 9992
FORMAT( /
' Tests performed: ', / 4x,
735 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
736 $
' r is a right eigenvector and ', a,
' means ', a,
'.',
737 $ /
' 1 = max | ( b A - a B )', a,
' l | / const.',
738 $ /
' 2 = max | ( b A - a B ) r | / const.',
739 $ /
' 3 = max ( Sest/Stru, Stru/Sest ) ',
740 $
' over all eigenvalues', /
741 $
' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
742 $
' over the 1st and 5th eigenvectors', / )
744 9991
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
745 $ i2,
', IWY=', i2,
', result ', i2,
' is', 0p, f8.2 )
746 9990
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
747 $ i2,
', IWY=', i2,
', result ', i2,
' is', 1p, d10.3 )
748 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
749 $
' result ', i2,
' is', 0p, f8.2 )
750 9988
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
751 $
' result ', i2,
' is', 1p, d10.3 )
752 9987
FORMAT(
' DDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
753 $ i6,
', Input example #', i2,
')' )
755 9986
FORMAT(
' DDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
756 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
757 $
'N=', i6,
', Input Example #', i2,
')' )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine dlatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
DLATM6
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ddrgvx(NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI, ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE, S, DTRU, DIF, DIFTRU, WORK, LWORK, IWORK, LIWORK, RESULT, BWORK, INFO)
DDRGVX
subroutine dget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
DGET52
subroutine dggevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO, IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, BWORK, INFO)
DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...