297 SUBROUTINE sdrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
298 $ ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE,
299 $ RSCALE, S, STRU, DIF, DIFTRU, WORK, LWORK,
300 $ IWORK, LIWORK, RESULT, BWORK, INFO )
307 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
314 REAL A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
315 $ alphar( * ), b( lda, * ), beta( * ),
316 $ bi( lda, * ), dif( * ), diftru( * ),
317 $ lscale( * ), result( 4 ), rscale( * ), s( * ),
318 $ stru( * ), vl( lda, * ), vr( lda, * ),
325 REAL ZERO, ONE, TEN, TNTH, HALF
326 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0, ten = 1.0e+1,
327 $ tnth = 1.0e-1, half = 0.5d+0 )
330 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
331 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
332 REAL ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
341 EXTERNAL ILAENV, SLAMCH, SLANGE
347 INTRINSIC abs, max, sqrt
357 IF( nsize.LT.0 )
THEN
359 ELSE IF( thresh.LT.zero )
THEN
361 ELSE IF( nin.LE.0 )
THEN
363 ELSE IF( nout.LE.0 )
THEN
365 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
367 ELSE IF( liwork.LT.nmax+6 )
THEN
379 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
380 minwrk = 2*nmax*nmax + 12*nmax + 16
381 maxwrk = 6*nmax + nmax*ilaenv( 1,
'SGEQRF',
' ', nmax, 1, nmax,
383 maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
387 IF( lwork.LT.minwrk )
391 CALL xerbla(
'SDRGVX', -info )
411 weight( 4 ) = one / weight( 2 )
412 weight( 5 ) = one / weight( 1 )
422 CALL slatm6( iptype, 5, a, lda, b, vr, lda, vl,
423 $ lda, weight( iwa ), weight( iwb ),
424 $ weight( iwx ), weight( iwy ), stru,
431 CALL slacpy(
'F', n, n, a, lda, ai, lda )
432 CALL slacpy(
'F', n, n, b, lda, bi, lda )
434 CALL sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi,
435 $ lda, alphar, alphai, beta, vl, lda,
436 $ vr, lda, ilo, ihi, lscale, rscale,
437 $ anorm, bnorm, s, dif, work, lwork,
438 $ iwork, bwork, linfo )
439 IF( linfo.NE.0 )
THEN
441 WRITE( nout, fmt = 9999 )
'SGGEVX', linfo, n,
448 CALL slacpy(
'Full', n, n, ai, lda, work, n )
449 CALL slacpy(
'Full', n, n, bi, lda, work( n*n+1 ),
451 abnorm = slange(
'Fro', n, 2*n, work, n, work )
456 CALL sget52( .true., n, a, lda, b, lda, vl, lda,
457 $ alphar, alphai, beta, work,
459 IF( result( 2 ).GT.thresh )
THEN
460 WRITE( nout, fmt = 9998 )
'Left',
'SGGEVX',
461 $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
465 CALL sget52( .false., n, a, lda, b, lda, vr, lda,
466 $ alphar, alphai, beta, work,
468 IF( result( 3 ).GT.thresh )
THEN
469 WRITE( nout, fmt = 9998 )
'Right',
'SGGEVX',
470 $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
477 IF( s( i ).EQ.zero )
THEN
478 IF( stru( i ).GT.abnorm*ulp )
479 $ result( 3 ) = ulpinv
480 ELSE IF( stru( i ).EQ.zero )
THEN
481 IF( s( i ).GT.abnorm*ulp )
482 $ result( 3 ) = ulpinv
484 work( i ) = max( abs( stru( i ) / s( i ) ),
485 $ abs( s( i ) / stru( i ) ) )
486 result( 3 ) = max( result( 3 ), work( i ) )
493 IF( dif( 1 ).EQ.zero )
THEN
494 IF( diftru( 1 ).GT.abnorm*ulp )
495 $ result( 4 ) = ulpinv
496 ELSE IF( diftru( 1 ).EQ.zero )
THEN
497 IF( dif( 1 ).GT.abnorm*ulp )
498 $ result( 4 ) = ulpinv
499 ELSE IF( dif( 5 ).EQ.zero )
THEN
500 IF( diftru( 5 ).GT.abnorm*ulp )
501 $ result( 4 ) = ulpinv
502 ELSE IF( diftru( 5 ).EQ.zero )
THEN
503 IF( dif( 5 ).GT.abnorm*ulp )
504 $ result( 4 ) = ulpinv
506 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
507 $ abs( dif( 1 ) / diftru( 1 ) ) )
508 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
509 $ abs( dif( 5 ) / diftru( 5 ) ) )
510 result( 4 ) = max( ratio1, ratio2 )
518 IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
519 $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
525 IF( nerrs.EQ.0 )
THEN
526 WRITE( nout, fmt = 9997 )
'SXV'
532 WRITE( nout, fmt = 9995 )
533 WRITE( nout, fmt = 9994 )
534 WRITE( nout, fmt = 9993 )
538 WRITE( nout, fmt = 9992 )
'''',
543 IF( result( j ).LT.10000.0 )
THEN
544 WRITE( nout, fmt = 9991 )iptype, iwa,
545 $ iwb, iwx, iwy, j, result( j )
547 WRITE( nout, fmt = 9990 )iptype, iwa,
548 $ iwb, iwx, iwy, j, result( j )
568 READ( nin, fmt = *,
END = 150 )n
572 READ( nin, fmt = * )( a( i, j ), j = 1, n )
575 READ( nin, fmt = * )( b( i, j ), j = 1, n )
577 READ( nin, fmt = * )( stru( i ), i = 1, n )
578 READ( nin, fmt = * )( diftru( i ), i = 1, n )
586 CALL slacpy(
'F', n, n, a, lda, ai, lda )
587 CALL slacpy(
'F', n, n, b, lda, bi, lda )
589 CALL sggevx(
'N',
'V',
'V',
'B', n, ai, lda, bi, lda, alphar,
590 $ alphai, beta, vl, lda, vr, lda, ilo, ihi, lscale,
591 $ rscale, anorm, bnorm, s, dif, work, lwork, iwork,
594 IF( linfo.NE.0 )
THEN
596 WRITE( nout, fmt = 9987 )
'SGGEVX', linfo, n, nptknt
602 CALL slacpy(
'Full', n, n, ai, lda, work, n )
603 CALL slacpy(
'Full', n, n, bi, lda, work( n*n+1 ), n )
604 abnorm =
slange(
'Fro', n, 2*n, work, n, work )
609 CALL sget52( .true., n, a, lda, b, lda, vl, lda, alphar, alphai,
610 $ beta, work, result( 1 ) )
611 IF( result( 2 ).GT.thresh )
THEN
612 WRITE( nout, fmt = 9986 )
'Left',
'SGGEVX', result( 2 ), n,
617 CALL sget52( .false., n, a, lda, b, lda, vr, lda, alphar, alphai,
618 $ beta, work, result( 2 ) )
619 IF( result( 3 ).GT.thresh )
THEN
620 WRITE( nout, fmt = 9986 )
'Right',
'SGGEVX', result( 3 ), n,
628 IF( s( i ).EQ.zero )
THEN
629 IF( stru( i ).GT.abnorm*ulp )
630 $ result( 3 ) = ulpinv
631 ELSE IF( stru( i ).EQ.zero )
THEN
632 IF( s( i ).GT.abnorm*ulp )
633 $ result( 3 ) = ulpinv
635 work( i ) = max( abs( stru( i ) / s( i ) ),
636 $ abs( s( i ) / stru( i ) ) )
637 result( 3 ) = max( result( 3 ), work( i ) )
644 IF( dif( 1 ).EQ.zero )
THEN
645 IF( diftru( 1 ).GT.abnorm*ulp )
646 $ result( 4 ) = ulpinv
647 ELSE IF( diftru( 1 ).EQ.zero )
THEN
648 IF( dif( 1 ).GT.abnorm*ulp )
649 $ result( 4 ) = ulpinv
650 ELSE IF( dif( 5 ).EQ.zero )
THEN
651 IF( diftru( 5 ).GT.abnorm*ulp )
652 $ result( 4 ) = ulpinv
653 ELSE IF( diftru( 5 ).EQ.zero )
THEN
654 IF( dif( 5 ).GT.abnorm*ulp )
655 $ result( 4 ) = ulpinv
657 ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
658 $ abs( dif( 1 ) / diftru( 1 ) ) )
659 ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
660 $ abs( dif( 5 ) / diftru( 5 ) ) )
661 result( 4 ) = max( ratio1, ratio2 )
669 IF( result( j ).GE.thrsh2 )
THEN
674 IF( nerrs.EQ.0 )
THEN
675 WRITE( nout, fmt = 9997 )
'SXV'
681 WRITE( nout, fmt = 9996 )
685 WRITE( nout, fmt = 9992 )
'''',
'transpose',
''''
689 IF( result( j ).LT.10000.0 )
THEN
690 WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
692 WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
704 CALL alasvm(
'SXV', nout, nerrs, ntestt, 0 )
710 9999
FORMAT(
' SDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
711 $ i6,
', JTYPE=', i6,
')' )
713 9998
FORMAT(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
714 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
715 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
716 $
', IWX=', i5,
', IWY=', i5 )
718 9997
FORMAT( / 1x, a3,
' -- Real Expert Eigenvalue/vector',
719 $
' problem driver' )
721 9996
FORMAT(
' Input Example' )
723 9995
FORMAT(
' Matrix types: ', / )
725 9994
FORMAT(
' TYPE 1: Da is diagonal, Db is identity, ',
726 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
727 $ /
' YH and X are left and right eigenvectors. ', / )
729 9993
FORMAT(
' TYPE 2: Da is quasi-diagonal, Db is identity, ',
730 $ /
' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
731 $ /
' YH and X are left and right eigenvectors. ', / )
733 9992
FORMAT( /
' Tests performed: ', / 4x,
734 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
735 $
' r is a right eigenvector and ', a,
' means ', a,
'.',
736 $ /
' 1 = max | ( b A - a B )', a,
' l | / const.',
737 $ /
' 2 = max | ( b A - a B ) r | / const.',
738 $ /
' 3 = max ( Sest/Stru, Stru/Sest ) ',
739 $
' over all eigenvalues', /
740 $
' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
741 $
' over the 1st and 5th eigenvectors', / )
743 9991
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
744 $ i2,
', IWY=', i2,
', result ', i2,
' is', 0p, f8.2 )
745 9990
FORMAT(
' Type=', i2,
',',
' IWA=', i2,
', IWB=', i2,
', IWX=',
746 $ i2,
', IWY=', i2,
', result ', i2,
' is', 1p, e10.3 )
747 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
748 $
' result ', i2,
' is', 0p, f8.2 )
749 9988
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
750 $
' result ', i2,
' is', 1p, e10.3 )
751 9987
FORMAT(
' SDRGVX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
752 $ i6,
', Input example #', i2,
')' )
754 9986
FORMAT(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
755 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
756 $
'N=', i6,
', Input Example #', i2,
')' )
subroutine alasvm(type, nout, nfail, nrun, nerrs)
ALASVM
subroutine xerbla(srname, info)
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 slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
real function slange(norm, m, n, a, lda, work)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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
subroutine slatm6(type, n, a, lda, b, x, ldx, y, ldy, alpha, beta, wx, wy, s, dif)
SLATM6