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 )
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 dget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
DGET52
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
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
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine dlatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
DLATM6
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 ...
double precision function dlamch(CMACH)
DLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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 ...