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,
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,
'.', / 9
x,
'N=',
714 $ i6,
', JTYPE=', i6,
')' )
716 9998
FORMAT(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
717 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9
x,
718 $
'N=', i6,
', JTYPE=', i6,
', IWA=', i5,
', IWB=', i5,
719 $
', IWX=', i5,
', IWY=', i5 )
721 9997
FORMAT( / 1
x, 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: ', / 4
x,
737 $
' a is alpha, b is beta, l is a left eigenvector, ', / 4
x,
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,
'.', / 9
x,
'N=',
755 $ i6,
', Input example #', i2,
')' )
757 9986
FORMAT(
' SDRGVX: ', a,
' Eigenvectors from ', a,
' incorrectly ',
758 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9
x,
759 $
'N=', i6,
', Input Example #', i2,
')' )
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
REAL function slamch(CMACH)
SLAMCH
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 sget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
SGET52
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
INTEGER function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
SLATM6