299 SUBROUTINE sdrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
300 $ alphar, alphai, beta, vl, vr, ilo, ihi, lscale,
301 $ rscale, s, stru, dif, diftru, work, lwork,
302 $ iwork, liwork, result, bwork, info )
310 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
317 REAL A( lda, * ), AI( lda, * ), ALPHAI( * ),
318 $ alphar( * ), b( lda, * ), beta( * ),
319 $ bi( lda, * ), dif( * ), diftru( * ),
320 $ lscale( * ), result( 4 ), rscale( * ), s( * ),
321 $ stru( * ), vl( lda, * ), vr( lda, * ),
328 REAL ZERO, ONE, TEN, TNTH, HALF
329 parameter ( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
330 $ tnth = 1.0e-1, half = 0.5d+0 )
333 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
334 $ maxwrk, minwrk, n, nerrs, nmax, nptknt, ntestt
335 REAL ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
344 EXTERNAL ilaenv, slamch, slange
350 INTRINSIC abs, max, sqrt
360 IF( nsize.LT.0 )
THEN
362 ELSE IF( thresh.LT.zero )
THEN
364 ELSE IF( nin.LE.0 )
THEN
366 ELSE IF( nout.LE.0 )
THEN
368 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
370 ELSE IF( liwork.LT.nmax+6 )
THEN
382 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
383 minwrk = 2*nmax*nmax + 12*nmax + 16
384 maxwrk = 6*nmax + nmax*ilaenv( 1,
'SGEQRF',
' ', nmax, 1, nmax,
386 maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
390 IF( lwork.LT.minwrk )
394 CALL xerbla(
'SDRGVX', -info )
414 weight( 4 ) = one / weight( 2 )
415 weight( 5 ) = one / weight( 1 )
425 CALL slatm6( iptype, 5, a, lda, b, vr, lda, vl,
426 $ lda, weight( iwa ), weight( iwb ),
427 $ weight( iwx ), weight( iwy ), stru,
434 CALL slacpy(
'F', n, n, a, lda, ai, lda )
435 CALL slacpy(
'F', n, n, b, lda, bi, lda )
437 CALL sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
438 $ lda, alphar, alphai, beta, vl, lda,
439 $ vr, lda, ilo, ihi, lscale, rscale,
440 $ anorm, bnorm, s, dif, work, lwork,
441 $ iwork, bwork, linfo )
442 IF( linfo.NE.0 )
THEN
444 WRITE( nout, fmt = 9999 )
'SGGEVX', linfo, n,
451 CALL slacpy(
'Full', n, n, ai, lda, work, n )
452 CALL slacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
454 abnorm = slange(
'Fro', n, 2*n, work, n, work )
459 CALL sget52( .true., n, a, lda, b, lda, vl, lda,
460 $ alphar, alphai, beta, work,
462 IF( result( 2 ).GT.thresh )
THEN
463 WRITE( nout, fmt = 9998 )
'Left',
'SGGEVX',
464 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
468 CALL sget52( .false., n, a, lda, b, lda, vr, lda,
469 $ alphar, alphai, beta, work,
471 IF( result( 3 ).GT.thresh )
THEN
472 WRITE( nout, fmt = 9998 )
'Right',
'SGGEVX',
473 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
480 IF( s( i ).EQ.zero )
THEN
481 IF( stru( i ).GT.abnorm*ulp )
482 $ result( 3 ) = ulpinv
483 ELSE IF( stru( i ).EQ.zero )
THEN
484 IF( s( i ).GT.abnorm*ulp )
485 $ result( 3 ) = ulpinv
487 work( i ) = max( abs( stru( i ) / s( i ) ),
488 $ abs( s( i ) / stru( i ) ) )
489 result( 3 ) = max( result( 3 ), work( i ) )
496 IF( dif( 1 ).EQ.zero )
THEN
497 IF( diftru( 1 ).GT.abnorm*ulp )
498 $ result( 4 ) = ulpinv
499 ELSE IF( diftru( 1 ).EQ.zero )
THEN
500 IF( dif( 1 ).GT.abnorm*ulp )
501 $ result( 4 ) = ulpinv
502 ELSE IF( dif( 5 ).EQ.zero )
THEN
503 IF( diftru( 5 ).GT.abnorm*ulp )
504 $ result( 4 ) = ulpinv
505 ELSE IF( diftru( 5 ).EQ.zero )
THEN
506 IF( dif( 5 ).GT.abnorm*ulp )
507 $ result( 4 ) = ulpinv
509 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
510 $ abs( dif( 1 ) / diftru( 1 ) ) )
511 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
512 $ abs( dif( 5 ) / diftru( 5 ) ) )
513 result( 4 ) = max( ratio1, ratio2 )
521 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
522 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
528 IF( nerrs.EQ.0 )
THEN
529 WRITE( nout, fmt = 9997 )
'SXV'
535 WRITE( nout, fmt = 9995 )
536 WRITE( nout, fmt = 9994 )
537 WRITE( nout, fmt = 9993 )
541 WRITE( nout, fmt = 9992 )
'''',
546 IF( result( j ).LT.10000.0 )
THEN
547 WRITE( nout, fmt = 9991 )iptype, iwa,
548 $ iwb, iwx, iwy, j, result( j )
550 WRITE( nout, fmt = 9990 )iptype, iwa,
551 $ iwb, iwx, iwy, j, result( j )
571 READ( nin, fmt = *, end = 150 )n
575 READ( nin, fmt = * )( a( i, j ), j = 1, n )
578 READ( nin, fmt = * )( b( i, j ), j = 1, n )
580 READ( nin, fmt = * )( stru( i ), i = 1, n )
581 READ( nin, fmt = * )( diftru( i ), i = 1, n )
589 CALL slacpy(
'F', n, n, a, lda, ai, lda )
590 CALL slacpy(
'F', n, n, b, lda, bi, lda )
592 CALL sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alphar,
593 $ alphai, beta, vl, lda, vr, lda, ilo, ihi, lscale,
594 $ rscale, anorm, bnorm, s, dif, work, lwork, iwork,
597 IF( linfo.NE.0 )
THEN
599 WRITE( nout, fmt = 9987 )
'SGGEVX', linfo, n, nptknt
605 CALL slacpy(
'Full', n, n, ai, lda, work, n )
606 CALL slacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
607 abnorm = slange(
'Fro', n, 2*n, work, n, work )
612 CALL sget52( .true., n, a, lda, b, lda, vl, lda, alphar, alphai,
613 $ beta, work, result( 1 ) )
614 IF( result( 2 ).GT.thresh )
THEN
615 WRITE( nout, fmt = 9986 )
'Left',
'SGGEVX', result( 2 ), n,
620 CALL sget52( .false., n, a, lda, b, lda, vr, lda, alphar, alphai,
621 $ beta, work, result( 2 ) )
622 IF( result( 3 ).GT.thresh )
THEN
623 WRITE( nout, fmt = 9986 )
'Right',
'SGGEVX', result( 3 ), n,
631 IF( s( i ).EQ.zero )
THEN
632 IF( stru( i ).GT.abnorm*ulp )
633 $ result( 3 ) = ulpinv
634 ELSE IF( stru( i ).EQ.zero )
THEN
635 IF( s( i ).GT.abnorm*ulp )
636 $ result( 3 ) = ulpinv
638 work( i ) = max( abs( stru( i ) / s( i ) ),
639 $ abs( s( i ) / stru( i ) ) )
640 result( 3 ) = max( result( 3 ), work( i ) )
647 IF( dif( 1 ).EQ.zero )
THEN
648 IF( diftru( 1 ).GT.abnorm*ulp )
649 $ result( 4 ) = ulpinv
650 ELSE IF( diftru( 1 ).EQ.zero )
THEN
651 IF( dif( 1 ).GT.abnorm*ulp )
652 $ result( 4 ) = ulpinv
653 ELSE IF( dif( 5 ).EQ.zero )
THEN
654 IF( diftru( 5 ).GT.abnorm*ulp )
655 $ result( 4 ) = ulpinv
656 ELSE IF( diftru( 5 ).EQ.zero )
THEN
657 IF( dif( 5 ).GT.abnorm*ulp )
658 $ result( 4 ) = ulpinv
660 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
661 $ abs( dif( 1 ) / diftru( 1 ) ) )
662 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
663 $ abs( dif( 5 ) / diftru( 5 ) ) )
664 result( 4 ) = max( ratio1, ratio2 )
672 IF( result( j ).GE.thrsh2 )
THEN
677 IF( nerrs.EQ.0 )
THEN
678 WRITE( nout, fmt = 9997 )
'SXV'
684 WRITE( nout, fmt = 9996 )
688 WRITE( nout, fmt = 9992 )
'''',
'transpose',
''''
692 IF( result( j ).LT.10000.0 )
THEN
693 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
695 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
707 CALL alasvm(
'SXV', nout, nerrs, ntestt, 0 )
713 9999
FORMAT(
' SDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
714 $ i6,
', JTYPE=', i6,
')' )
716 9998
FORMAT(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
717 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
718 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
719 $
', IWX=', i5,
', IWY=', i5 )
721 9997
FORMAT( / 1x, a3,
' -- Real Expert Eigenvalue/vector',
722 $
' problem driver' )
724 9996
FORMAT(
' Input Example' )
726 9995
FORMAT(
' Matrix types: ', / )
728 9994
FORMAT(
' TYPE 1: Da is diagonal, Db is identity, ',
729 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
730 $ /
' YH and X are left and right eigenvectors. ', / )
732 9993
FORMAT(
' TYPE 2: Da is quasi-diagonal, Db is identity, ',
733 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
734 $ /
' YH and X are left and right eigenvectors. ', / )
736 9992
FORMAT( /
' Tests performed: ', / 4x,
737 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
738 $
' r is a right eigenvector and ', a,
' means ', a,
'.',
739 $ /
' 1 = max | ( b A - a B )', a,
' l | / const.',
740 $ /
' 2 = max | ( b A - a B ) r | / const.',
741 $ /
' 3 = max ( Sest/Stru, Stru/Sest ) ',
742 $
' over all eigenvalues', /
743 $
' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
744 $
' over the 1st and 5th eigenvectors', / )
746 9991
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
747 $ i2,
', IWY=', i2,
', result ', i2,
' is', 0p, f8.2 )
748 9990
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
749 $ i2,
', IWY=', i2,
', result ', i2,
' is', 1p, e10.3 )
750 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
751 $
' result ', i2,
' is', 0p, f8.2 )
752 9988
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
753 $
' result ', i2,
' is', 1p, e10.3 )
754 9987
FORMAT(
' SDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
755 $ i6,
', Input example #', i2,
')' )
757 9986
FORMAT(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
758 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
759 $
'N=', i6,
', Input Example #', i2,
')' )
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine sggevx(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)
SGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
subroutine slatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
SLATM6
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sdrgvx(NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI, ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE, S, STRU, DIF, DIFTRU, WORK, LWORK, IWORK, LIWORK, RESULT, BWORK, INFO)
SDRGVX
subroutine sget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
SGET52