402 SUBROUTINE ddrvev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
403 $ NOUNIT, A, LDA, H, WR, WI, WR1, WI1, VL, LDVL,
404 $ VR, LDVR, LRE, LDLRE, RESULT, WORK, NWORK,
412 INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NOUNIT, NSIZES,
414 DOUBLE PRECISION THRESH
418 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
419 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
420 $ result( 7 ), vl( ldvl, * ), vr( ldvr, * ),
421 $ wi( * ), wi1( * ), work( * ), wr( * ), wr1( * )
427 DOUBLE PRECISION ZERO, ONE
428 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
430 parameter( two = 2.0d0 )
432 parameter( maxtyp = 21 )
437 INTEGER IINFO, IMODE, ITYPE, IWK, J, , JJ, JSIZE,
438 $ jtype, mtypes, n, nerrs, nfail, nmax, nnwork,
439 $ ntest, ntestf, ntestt
440 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TNRM,
441 $ ULP, ULPINV, UNFL, VMX, VRMX, VTST
444 CHARACTER ADUMMA( 1 )
445 INTEGER ( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
446 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
448 DOUBLE PRECISION DUM( 1 ), RES( 2 )
451 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
452 EXTERNAL DLAMCH, DLAPY2, DNRM2
459 INTRINSIC abs,
max,
min, sqrt
462 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
463 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
465 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
466 $ 1, 5, 5, 5, 4, 3, 1 /
467 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
471 path( 1: 1 )
'Double precision'
485 nmax =
max( nmax, nn( j ) )
492 IF( nsizes.LT.0 )
THEN
494 ELSE IF( badnn )
THEN
496 ELSE IF( ntypes.LT.0 )
THEN
498 ELSE IF( thresh.LT.zero )
THEN
500 ELSE IF( nounit.LE.0 )
THEN
502 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
504 ELSE IF( ldvl.LT.1 .OR. ldvl.LT.nmax )
THEN
506 ELSE IF( ldvr.LT.1 .OR. ldvr.LT.nmax )
THEN
508 ELSE IF( ldlre.LT.1 .OR. ldlre.LT.nmax )
THEN
510 ELSE IF( 5*nmax+2*nmax**2.GT.nwork )
THEN
515 CALL xerbla(
'DDRVEV', -info )
521 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
526 unfl = dlamch(
'Safe minimum' )
529 ulp = dlamch(
'Precision' )
538 DO 270 jsize = 1, nsizes
540 IF( nsizes.NE.1 )
THEN
541 mtypes =
min( maxtyp, ntypes )
543 mtypes =
min( maxtyp+1, ntypes )
546 DO 260 jtype = 1, mtypes
547 IF( .NOT.dotype( jtype ) )
553 ioldsd( j ) = iseed( j )
572 IF( mtypes.GT.maxtyp )
575 itype = ktype( jtype )
576 imode = kmode( jtype )
580 GO TO ( 30, 40, 50 )kmagn( jtype )
596 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
604 IF( itype.EQ.1 )
THEN
607 ELSE IF( itype.EQ.2 )
THEN
612 a( jcol, jcol ) = anorm
615 ELSE IF( itype.EQ.3 )
THEN
620 a( jcol, jcol ) = anorm
622 $ a( jcol, jcol-1 ) = one
625 ELSE IF( itype.EQ.4 )
THEN
629 CALL dlatms( n, n,
'S', iseed,
'S', work, imode
630 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
633 ELSE IF( itype.EQ.5 )
THEN
637 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond
638 $ anorm, n, n,
'N', a, lda, work( n+1 ),
641 ELSE IF( itype.EQ.6 )
THEN
645 IF( kconds( jtype ).EQ.1 )
THEN
647 ELSE IF( kconds( jtype ).EQ.2 )
THEN
654 CALL dlatme( n,
'S', iseed, work, imode, cond, one,
655 $ adumma,
'T',
'T',
'T', work( n+1 ), 4,
656 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
659 ELSE IF( itype.EQ.7 )
THEN
663 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
664 $
'T',
'N', work( n+1 ), 1, one,
665 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
666 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
668 ELSE IF( itype.EQ.8 )
THEN
672 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
674 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
675 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
677 ELSE IF( itype.EQ.9 )
THEN
681 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
682 $
'T',
'N', work( n+1 ), 1, one,
683 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
684 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
686 CALL dlaset(
'Full', 2, n, zero, zero, a, lda )
687 CALL dlaset(
'Full', n-3, 1, zero, zero, a( 3, 1 ),
689 CALL dlaset(
'Full', n-3, 2, zero, zero, a( 3, n-1 ),
691 CALL dlaset(
'Full', 1, n, zero, zero, a( n, 1 ),
695 ELSE IF( itype.EQ.10 )
THEN
699 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
700 $
'T',
'N', work( n+1 ), 1, one,
701 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
702 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
709 IF( iinfo.NE.0 )
THEN
710 WRITE( nounit, fmt = 9993 )
'Generator', iinfo, n, jtype,
724 nnwork = 5*n + 2*n**2
726 nnwork =
max( nnwork, 1 )
736 CALL dlacpy(
'F', n, n, a, lda, h, lda )
737 CALL dgeev(
'V',
'V', n, h, lda, wr, wi, vl, ldvl, vr,
738 $ ldvr, work, nnwork, iinfo )
739 IF( iinfo.NE.0 )
THEN
741 WRITE( nounit, fmt = 9993 )
'DGEEV1', iinfo, n, jtype,
749 CALL dget22(
'N',
'N',
'N', n, a, lda, vr, ldvr, wr, wi,
751 result( 1 ) = res( 1 )
755 CALL dget22(
'T',
'N',
'T', n, a, lda, vl, ldvl, wr, wi,
757 result( 2 ) = res( 1 )
763 IF( wi( j ).EQ.zero )
THEN
764 tnrm = dnrm2( n, vr( 1, j ), 1 )
765 ELSE IF( wi( j ).GT.zero )
THEN
766 tnrm = dlapy2( dnrm2( n, vr( 1, j ), 1 ),
767 $ dnrm2( n, vr( 1, j+1 ), 1 ) )
769 result( 3 ) =
max( result( 3 ),
770 $
min( ulpinv, abs( tnrm-one ) / ulp ) )
771 IF( wi( j ).GT.zero )
THEN
775 vtst = dlapy2( vr( jj, j ), vr( jj, j+1 ) )
778 IF( vr( jj, j+1 ).EQ.zero .AND.
779 $ abs( vr( jj, j ) ).GT.vrmx )
780 $ vrmx = abs( vr( jj, j ) )
782 IF( vrmx / vmx.LT.one-two*ulp )
783 $ result( 3 ) = ulpinv
791 IF( wi( j ).EQ.zero )
THEN
792 tnrm = dnrm2( n, vl( 1, j ), 1 )
793 ELSE IF( wi( j ).GT.zero )
THEN
794 tnrm = dlapy2( dnrm2( n, vl( 1, j ), 1 ),
795 $ dnrm2( n, vl( 1, j+1 ), 1 ) )
797 result( 4 ) =
max( result( 4 ),
798 $
min( ulpinv, abs( tnrm-one ) / ulp ) )
799 IF( wi( j ).GT.zero )
THEN
803 vtst = dlapy2( vl( jj, j ), vl( jj, j+1 ) )
806 IF( vl( jj, j+1 ).EQ.zero .AND.
807 $ abs( vl( jj, j ) ).GT.vrmx )
808 $ vrmx = abs( vl( jj, j ) )
810 IF( vrmx / vmx.LT.one-two*ulp )
811 $ result( 4 ) = ulpinv
817 CALL dlacpy(
'F', n, n, a, lda, h, lda )
818 CALL dgeev(
'N',
'N', n, h, lda, wr1, wi1, dum, 1, dum,
819 $ 1, work, nnwork, iinfo )
820 IF( iinfo.NE.0 )
THEN
822 WRITE( nounit, fmt = 9993 )
'DGEEV2', iinfo, n, jtype,
831 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
832 $ result( 5 ) = ulpinv
837 CALL dlacpy(
'F', n, n, a, lda, h, lda )
838 CALL dgeev(
'N',
'V', n, h, lda, wr1, wi1, dum, 1, lre,
839 $ ldlre, work, nnwork, iinfo )
840 IF( iinfo.NE.0 )
THEN
842 WRITE( nounit, fmt = 9993 )
'DGEEV3', iinfo, n, jtype,
851 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
852 $ result( 5 ) = ulpinv
859 IF( vr( j, jj ).NE.lre( j, jj ) )
860 $ result( 6 ) = ulpinv
866 CALL dlacpy( 'f
', N, N, A, LDA, H, LDA )
867 CALL DGEEV( 'v
', 'n
', N, H, LDA, WR1, WI1, LRE, LDLRE,
868 $ DUM, 1, WORK, NNWORK, IINFO )
869.NE.
IF( IINFO0 ) THEN
871 WRITE( NOUNIT, FMT = 9993 )'dgeev4
', IINFO, N, JTYPE,
880.NE..OR..NE.
IF( WR( J )WR1( J ) WI( J )WI1( J ) )
881 $ RESULT( 5 ) = ULPINV
888.NE.
IF( VL( J, JJ )LRE( J, JJ ) )
889 $ RESULT( 7 ) = ULPINV
900.GE.
IF( RESULT( J )ZERO )
902.GE.
IF( RESULT( J )THRESH )
907 $ NTESTF = NTESTF + 1
908.EQ.
IF( NTESTF1 ) THEN
909 WRITE( NOUNIT, FMT = 9999 )PATH
910 WRITE( NOUNIT, FMT = 9998 )
911 WRITE( NOUNIT, FMT = 9997 )
912 WRITE( NOUNIT, FMT = 9996 )
913 WRITE( NOUNIT, FMT = 9995 )THRESH
918.GE.
IF( RESULT( J )THRESH ) THEN
919 WRITE( NOUNIT, FMT = 9994 )N, IWK, IOLDSD, JTYPE,
924 NERRS = NERRS + NFAIL
925 NTESTT = NTESTT + NTEST
933 CALL DLASUM( PATH, NOUNIT, NERRS, NTESTT )
935 9999 FORMAT( / 1X, A3, ' -- real eigenvalue-eigenvector decomposition
',
936 $ ' driver
', / ' matrix types
' )
938 9998 FORMAT( / ' special matrices:
', / ' 1=zero matrix.
',
939 $ ' ', ' 5=diagonal: geometr. spaced entries.
',
940 $ / ' 2=identity matrix
', ' 6=diagona
',
941 $ 'l: clustered entries.
', / ' 3=transposed jordan block.
',
943 $ '4=diagonal: evenly spaced entries
', ' 8=diagonal: s
',
944 $ 'mall, evenly spaced.
' )
945 9997 FORMAT( ' dense, non-symmetric matrices:
', / ' 9=well-cond., ev
',
946 $ 'enly spaced eigenvals.
', ' 14=ill-cond., geomet. spaced e
',
947 $ 'igenals.
', / ' 10=well-cond.,
geom. spaced eigenvals.
',
949 $ 'itioned, clustered e.vals.
', ' 16=ill-cond.,
',
950 $ 'lex
', / ' 12=well-cond., random
complex ', 6X, ' ',
951 $ ' 17=ill-cond., large rand. complx
', / ' 13=ill-condi
',
952 $ 'tioned, evenly spaced.
', ' 18=ill-cond., small rand.
',
954 9996 FORMAT( ' 19=matrix with random o(1) entries.
', ' 21=matrix
',
955 $ 'with small
', / ' 20=matrix with large ran
',
956 $ 'dom entries.
', / )
957 9995 FORMAT( ' tests performed with test threshold =
', F8.2,
958 $ / / ' 1 = | a vr - vr w | / ( n |a| ulp )
',
959 $ / ' 2 = | transpose(a) vl - vl w | / ( n |a| ulp )
',
960 $ / ' 3 = | |vr(i)| - 1 | / ulp
',
961 $ / ' 4 = | |vl(i)| - 1 | / ulp
',
962 $ / ' 5 = 0
if w same no matter
if vr or vl computed,
',
963 $ ' 1/ulp otherwise
', /
964 $ ' 6 = 0
if vr same no matter
if vl computed,
',
965 $ ' 1/ulp otherwise
', /
966 $ ' 7 = 0
if vl same no matter
if vr computed,
',
967 $ ' 1/ulp otherwise
', / )
968 9994 FORMAT( ' n=
', I5, '', I2, ',
seed=
', 4( I4, ',
' ),
969 $ ' type ', I2, ', test(
', I2, ')=
', G10.3 )
970 9993 FORMAT( ' ddrvev:
', A, ' returned info=
', I6, '.
', / 9X, 'n=
',
971 $ I6, ', jtype=
', I6, ', iseed=(
', 3( I5, ',
' ), I5, ')
' )