296 SUBROUTINE ddrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
297 $ ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE,
298 $ RSCALE, S, DTRU, DIF, DIFTRU, WORK, LWORK,
299 $ IWORK, LIWORK, RESULT, BWORK, INFO )
306 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
308 DOUBLE PRECISION THRESH
313 DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
314 $ alphar( * ), b( lda, * ), beta( * ),
315 $ bi( lda, * ), dif( * ), diftru( * ), dtru( * ),
316 $ lscale( * ), result( 4 ), rscale( * ), s( * ),
317 $ vl( lda, * ), vr( lda, * ), work( * )
323 DOUBLE PRECISION ZERO, ONE, TEN, TNTH, HALF
324 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
325 $ tnth = 1.0d-1, half = 0.5d+0 )
328 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
329 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
330 DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
334 DOUBLE PRECISION WEIGHT( 5 )
338 DOUBLE PRECISION DLAMCH, DLANGE
339 EXTERNAL ILAENV, DLAMCH, DLANGE
345 INTRINSIC abs, max, sqrt
355 IF( nsize.LT.0 )
THEN
357 ELSE IF( thresh.LT.zero )
THEN
359 ELSE IF( nin.LE.0 )
THEN
361 ELSE IF( nout.LE.0 )
THEN
363 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
365 ELSE IF( liwork.LT.nmax+6 )
THEN
377 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
378 minwrk = 2*nmax*nmax + 12*nmax + 16
379 maxwrk = 6*nmax + nmax*ilaenv( 1,
'DGEQRF',
' ', nmax, 1, nmax,
381 maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
385 IF( lwork.LT.minwrk )
389 CALL xerbla(
'DDRGVX', -info )
409 weight( 4 ) = one / weight( 2 )
410 weight( 5 ) = one / weight( 1 )
420 CALL dlatm6( iptype, 5, a, lda, b, vr, lda, vl,
421 $ lda, weight( iwa ), weight( iwb ),
422 $ weight( iwx ), weight( iwy ), dtru,
429 CALL dlacpy(
'F', n, n, a, lda, ai, lda )
430 CALL dlacpy(
'F', n, n, b, lda, bi, lda )
432 CALL dggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
433 $ lda, alphar, alphai, beta, vl, lda,
434 $ vr, lda, ilo, ihi, lscale, rscale,
435 $ anorm, bnorm, s, dif, work, lwork,
436 $ iwork, bwork, linfo )
437 IF( linfo.NE.0 )
THEN
439 WRITE( nout, fmt = 9999 )
'DGGEVX', linfo, n,
446 CALL dlacpy(
'Full', n, n, ai, lda, work, n )
447 CALL dlacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
449 abnorm = dlange(
'Fro', n, 2*n, work, n, work )
454 CALL dget52( .true., n, a, lda, b, lda, vl, lda,
455 $ alphar, alphai, beta, work,
457 IF( result( 2 ).GT.thresh )
THEN
458 WRITE( nout, fmt = 9998 )
'Left',
'DGGEVX',
459 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
463 CALL dget52( .false., n, a, lda, b, lda, vr, lda,
464 $ alphar, alphai, beta, work,
466 IF( result( 3 ).GT.thresh )
THEN
467 WRITE( nout, fmt = 9998 )
'Right',
'DGGEVX',
468 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
475 IF( s( i ).EQ.zero )
THEN
476 IF( dtru( i ).GT.abnorm*ulp )
477 $ result( 3 ) = ulpinv
478 ELSE IF( dtru( i ).EQ.zero )
THEN
479 IF( s( i ).GT.abnorm*ulp )
480 $ result( 3 ) = ulpinv
482 work( i ) = max( abs( dtru( i ) / s( i ) ),
483 $ abs( s( i ) / dtru( i ) ) )
484 result( 3 ) = max( result( 3 ), work( i ) )
491 IF( dif( 1 ).EQ.zero )
THEN
492 IF( diftru( 1 ).GT.abnorm*ulp )
493 $ result( 4 ) = ulpinv
494 ELSE IF( diftru( 1 ).EQ.zero )
THEN
495 IF( dif( 1 ).GT.abnorm*ulp )
496 $ result( 4 ) = ulpinv
497 ELSE IF( dif( 5 ).EQ.zero )
THEN
498 IF( diftru( 5 ).GT.abnorm*ulp )
499 $ result( 4 ) = ulpinv
500 ELSE IF( diftru( 5 ).EQ.zero )
THEN
501 IF( dif( 5 ).GT.abnorm*ulp )
502 $ result( 4 ) = ulpinv
504 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
505 $ abs( dif( 1 ) / diftru( 1 ) ) )
506 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
507 $ abs( dif( 5 ) / diftru( 5 ) ) )
508 result( 4 ) = max( ratio1, ratio2 )
516 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
517 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
523 IF( nerrs.EQ.0 )
THEN
524 WRITE( nout, fmt = 9997 )
'DXV'
530 WRITE( nout, fmt = 9995 )
531 WRITE( nout, fmt = 9994 )
532 WRITE( nout, fmt = 9993 )
536 WRITE( nout, fmt = 9992 )
'''',
541 IF( result( j ).LT.10000.0d0 )
THEN
542 WRITE( nout, fmt = 9991 )iptype, iwa,
543 $ iwb, iwx, iwy, j, result( j )
545 WRITE( nout, fmt = 9990 )iptype, iwa,
546 $ iwb, iwx, iwy, j, result( j )
566 READ( nin, fmt = *,
END = 150 )n
570 READ( nin, fmt = * )( a( i, j ), j = 1, n )
573 READ( nin, fmt = * )( b( i, j ), j = 1, n )
575 READ( nin, fmt = * )( dtru( i ), i = 1, n )
576 READ( nin, fmt = * )( diftru( i ), i = 1, n )
584 CALL dlacpy(
'F', n, n, a, lda, ai, lda )
585 CALL dlacpy(
'F', n, n, b, lda, bi, lda )
587 CALL dggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alphar,
588 $ alphai, beta, vl, lda, vr, lda, ilo, ihi, lscale,
589 $ rscale, anorm, bnorm, s, dif, work, lwork, iwork,
592 IF( linfo.NE.0 )
THEN
594 WRITE( nout, fmt = 9987 )
'DGGEVX', linfo, n, nptknt
600 CALL dlacpy(
'Full', n, n, ai, lda, work, n )
601 CALL dlacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
602 abnorm =
dlange(
'Fro', n, 2*n, work, n, work )
607 CALL dget52( .true., n, a, lda, b, lda, vl, lda, alphar, alphai,
608 $ beta, work, result( 1 ) )
609 IF( result( 2 ).GT.thresh )
THEN
610 WRITE( nout, fmt = 9986 )
'Left',
'DGGEVX', result( 2 ), n,
615 CALL dget52( .false., n, a, lda, b, lda, vr, lda, alphar, alphai,
616 $ beta, work, result( 2 ) )
617 IF( result( 3 ).GT.thresh )
THEN
618 WRITE( nout, fmt = 9986 )
'Right',
'DGGEVX', result( 3 ), n,
626 IF( s( i ).EQ.zero )
THEN
627 IF( dtru( i ).GT.abnorm*ulp )
628 $ result( 3 ) = ulpinv
629 ELSE IF( dtru( i ).EQ.zero )
THEN
630 IF( s( i ).GT.abnorm*ulp )
631 $ result( 3 ) = ulpinv
633 work( i ) = max( abs( dtru( i ) / s( i ) ),
634 $ abs( s( i ) / dtru( i ) ) )
635 result( 3 ) = max( result( 3 ), work( i ) )
642 IF( dif( 1 ).EQ.zero )
THEN
643 IF( diftru( 1 ).GT.abnorm*ulp )
644 $ result( 4 ) = ulpinv
645 ELSE IF( diftru( 1 ).EQ.zero )
THEN
646 IF( dif( 1 ).GT.abnorm*ulp )
647 $ result( 4 ) = ulpinv
648 ELSE IF( dif( 5 ).EQ.zero )
THEN
649 IF( diftru( 5 ).GT.abnorm*ulp )
650 $ result( 4 ) = ulpinv
651 ELSE IF( diftru( 5 ).EQ.zero )
THEN
652 IF( dif( 5 ).GT.abnorm*ulp )
653 $ result( 4 ) = ulpinv
655 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
656 $ abs( dif( 1 ) / diftru( 1 ) ) )
657 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
658 $ abs( dif( 5 ) / diftru( 5 ) ) )
659 result( 4 ) = max( ratio1, ratio2 )
667 IF( result( j ).GE.thrsh2 )
THEN
672 IF( nerrs.EQ.0 )
THEN
673 WRITE( nout, fmt = 9997 )
'DXV'
679 WRITE( nout, fmt = 9996 )
683 WRITE( nout, fmt = 9992 )
'''',
'transpose',
''''
687 IF( result( j ).LT.10000.0d0 )
THEN
688 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
690 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
702 CALL alasvm(
'DXV', nout, nerrs, ntestt, 0 )
708 9999
FORMAT(
' DDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
709 $ i6,
', JTYPE=', i6,
')' )
711 9998
FORMAT(
' DDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
712 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
713 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
714 $
', IWX=', i5,
', IWY=', i5 )
716 9997
FORMAT( / 1x, a3,
' -- Real Expert Eigenvalue/vector',
717 $
' problem driver' )
719 9996
FORMAT(
' Input Example' )
721 9995
FORMAT(
' Matrix types: ', / )
723 9994
FORMAT(
' TYPE 1: Da is diagonal, Db is identity, ',
724 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
725 $ /
' YH and X are left and right eigenvectors. ', / )
727 9993
FORMAT(
' TYPE 2: Da is quasi-diagonal, Db is identity, ',
728 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
729 $ /
' YH and X are left and right eigenvectors. ', / )
731 9992
FORMAT( /
' Tests performed: ', / 4x,
732 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
733 $
' r is a right eigenvector and ', a,
' means ', a,
'.',
734 $ /
' 1 = max | ( b A - a B )', a,
' l | / const.',
735 $ /
' 2 = max | ( b A - a B ) r | / const.',
736 $ /
' 3 = max ( Sest/Stru, Stru/Sest ) ',
737 $
' over all eigenvalues', /
738 $
' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
739 $
' over the 1st and 5th eigenvectors', / )
741 9991
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
742 $ i2,
', IWY=', i2,
', result ', i2,
' is', 0p, f8.2 )
743 9990
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
744 $ i2,
', IWY=', i2,
', result ', i2,
' is', 1p, d10.3 )
745 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
746 $
' result ', i2,
' is', 0p, f8.2 )
747 9988
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
748 $
' result ', i2,
' is', 1p, d10.3 )
749 9987
FORMAT(
' DDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
750 $ i6,
', Input example #', i2,
')' )
752 9986
FORMAT(
' DDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
753 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
754 $
'N=', i6,
', Input Example #', i2,
')' )
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine xerbla(srname, info)
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 dlatm6(type, n, a, lda, b, x, ldx, y, ldy, alpha, beta, wx, wy, s, dif)
DLATM6
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
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
double precision function dlange(norm, m, n, a, lda, work)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...