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, JCOL, 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 IDUMMA( 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' )
528 ulp = dlamch(
'Precision' )
537 DO 270 jsize = 1, nsizes
539 IF( nsizes.NE.1 )
THEN
540 mtypes = min( maxtyp, ntypes )
542 mtypes = min( maxtyp+1, ntypes )
545 DO 260 jtype = 1, mtypes
546 IF( .NOT.dotype( jtype ) )
552 ioldsd( j ) = iseed( j )
571 IF( mtypes.GT.maxtyp )
574 itype = ktype( jtype )
575 imode = kmode( jtype )
579 GO TO ( 30, 40, 50 )kmagn( jtype )
595 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
603 IF( itype.EQ.1 )
THEN
606 ELSE IF( itype.EQ.2 )
THEN
611 a( jcol, jcol ) = anorm
614 ELSE IF( itype.EQ.3 )
THEN
619 a( jcol, jcol ) = anorm
621 $ a( jcol, jcol-1 ) = one
624 ELSE IF( itype.EQ.4 )
THEN
628 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
629 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
632 ELSE IF( itype.EQ.5 )
THEN
636 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
637 $ anorm, n, n,
'N', a, lda, work( n+1 ),
640 ELSE IF( itype.EQ.6 )
THEN
644 IF( kconds( jtype ).EQ.1 )
THEN
646 ELSE IF( kconds( jtype ).EQ.2 )
THEN
653 CALL dlatme( n,
'S', iseed, work, imode, cond, one,
654 $ adumma,
'T',
'T',
'T', work( n+1 ), 4,
655 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
658 ELSE IF( itype.EQ.7 )
THEN
662 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
663 $
'T',
'N', work( n+1 ), 1, one,
664 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
665 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
667 ELSE IF( itype.EQ.8 )
THEN
671 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
672 $
'T',
'N', work( n+1 ), 1, one,
673 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
674 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
676 ELSE IF( itype.EQ.9 )
THEN
680 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
681 $
'T',
'N', work( n+1 ), 1, one,
682 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
683 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
685 CALL dlaset(
'Full', 2, n, zero, zero, a, lda )
686 CALL dlaset(
'Full', n-3, 1, zero, zero, a( 3, 1 ),
688 CALL dlaset(
'Full', n-3, 2, zero, zero, a( 3, n-1 ),
690 CALL dlaset(
'Full', 1, n, zero, zero, a( n, 1 ),
694 ELSE IF( itype.EQ.10 )
THEN
698 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
699 $
'T',
'N', work( n+1 ), 1, one,
700 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
701 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
708 IF( iinfo.NE.0 )
THEN
709 WRITE( nounit, fmt = 9993 )
'Generator', iinfo, n, jtype,
723 nnwork = 5*n + 2*n**2
725 nnwork = max( nnwork, 1 )
735 CALL dlacpy(
'F', n, n, a, lda, h, lda )
736 CALL dgeev(
'V',
'V', n, h, lda, wr, wi, vl, ldvl, vr,
737 $ ldvr, work, nnwork, iinfo )
738 IF( iinfo.NE.0 )
THEN
740 WRITE( nounit, fmt = 9993 )
'DGEEV1', iinfo, n, jtype,
748 CALL dget22(
'N',
'N',
'N', n, a, lda, vr, ldvr, wr, wi,
750 result( 1 ) = res( 1 )
754 CALL dget22(
'T',
'N',
'T', n, a, lda, vl, ldvl, wr, wi,
756 result( 2 ) = res( 1 )
762 IF( wi( j ).EQ.zero )
THEN
763 tnrm = dnrm2( n, vr( 1, j ), 1 )
764 ELSE IF( wi( j ).GT.zero )
THEN
765 tnrm = dlapy2( dnrm2( n, vr( 1, j ), 1 ),
766 $ dnrm2( n, vr( 1, j+1 ), 1 ) )
768 result( 3 ) = max( result( 3 ),
769 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
770 IF( wi( j ).GT.zero )
THEN
774 vtst = dlapy2( vr( jj, j ), vr( jj, j+1 ) )
777 IF( vr( jj, j+1 ).EQ.zero .AND.
778 $ abs( vr( jj, j ) ).GT.vrmx )
779 $ vrmx = abs( vr( jj, j ) )
781 IF( vrmx / vmx.LT.one-two*ulp )
782 $ result( 3 ) = ulpinv
790 IF( wi( j ).EQ.zero )
THEN
791 tnrm = dnrm2( n, vl( 1, j ), 1 )
792 ELSE IF( wi( j ).GT.zero )
THEN
793 tnrm = dlapy2( dnrm2( n, vl( 1, j ), 1 ),
794 $ dnrm2( n, vl( 1, j+1 ), 1 ) )
796 result( 4 ) = max( result( 4 ),
797 $ min( ulpinv, abs( tnrm-one ) / ulp ) )
798 IF( wi( j ).GT.zero )
THEN
802 vtst = dlapy2( vl( jj, j ), vl( jj, j+1 ) )
805 IF( vl( jj, j+1 ).EQ.zero .AND.
806 $ abs( vl( jj, j ) ).GT.vrmx )
807 $ vrmx = abs( vl( jj, j ) )
809 IF( vrmx / vmx.LT.one-two*ulp )
810 $ result( 4 ) = ulpinv
816 CALL dlacpy(
'F', n, n, a, lda, h, lda )
817 CALL dgeev(
'N',
'N', n, h, lda, wr1, wi1, dum, 1, dum,
818 $ 1, work, nnwork, iinfo )
819 IF( iinfo.NE.0 )
THEN
821 WRITE( nounit, fmt = 9993 )
'DGEEV2', iinfo, n, jtype,
830 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
831 $ result( 5 ) = ulpinv
836 CALL dlacpy(
'F', n, n, a, lda, h, lda )
837 CALL dgeev(
'N',
'V', n, h, lda, wr1, wi1, dum, 1, lre,
838 $ ldlre, work, nnwork, iinfo )
839 IF( iinfo.NE.0 )
THEN
841 WRITE( nounit, fmt = 9993 )
'DGEEV3', iinfo, n, jtype,
850 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
851 $ result( 5 ) = ulpinv
858 IF( vr( j, jj ).NE.lre( j, jj ) )
859 $ result( 6 ) = ulpinv
865 CALL dlacpy(
'F', n, n, a, lda, h, lda )
866 CALL dgeev(
'V',
'N', n, h, lda, wr1, wi1, lre, ldlre,
867 $ dum, 1, work, nnwork, iinfo )
868 IF( iinfo.NE.0 )
THEN
870 WRITE( nounit, fmt = 9993 )
'DGEEV4', iinfo, n, jtype,
879 IF( wr( j ).NE.wr1( j ) .OR. wi( j ).NE.wi1( j ) )
880 $ result( 5 ) = ulpinv
887 IF( vl( j, jj ).NE.lre( j, jj ) )
888 $ result( 7 ) = ulpinv
899 IF( result( j ).GE.zero )
901 IF( result( j ).GE.thresh )
906 $ ntestf = ntestf + 1
907 IF( ntestf.EQ.1 )
THEN
908 WRITE( nounit, fmt = 9999 )path
909 WRITE( nounit, fmt = 9998 )
910 WRITE( nounit, fmt = 9997 )
911 WRITE( nounit, fmt = 9996 )
912 WRITE( nounit, fmt = 9995 )thresh
917 IF( result( j ).GE.thresh )
THEN
918 WRITE( nounit, fmt = 9994 )n, iwk, ioldsd, jtype,
923 nerrs = nerrs + nfail
924 ntestt = ntestt + ntest
932 CALL dlasum( path, nounit, nerrs, ntestt )
934 9999
FORMAT( / 1x, a3,
' -- Real Eigenvalue-Eigenvector Decomposition',
935 $
' Driver', /
' Matrix types (see DDRVEV for details): ' )
937 9998
FORMAT( /
' Special Matrices:', /
' 1=Zero matrix. ',
938 $
' ',
' 5=Diagonal: geometr. spaced entries.',
939 $ /
' 2=Identity matrix. ',
' 6=Diagona',
940 $
'l: clustered entries.', /
' 3=Transposed Jordan block. ',
941 $
' ',
' 7=Diagonal: large, evenly spaced.', /
' ',
942 $
'4=Diagonal: evenly spaced entries. ',
' 8=Diagonal: s',
943 $
'mall, evenly spaced.' )
944 9997
FORMAT(
' Dense, Non-Symmetric Matrices:', /
' 9=Well-cond., ev',
945 $
'enly spaced eigenvals.',
' 14=Ill-cond., geomet. spaced e',
946 $
'igenals.', /
' 10=Well-cond., geom. spaced eigenvals. ',
947 $
' 15=Ill-conditioned, clustered e.vals.', /
' 11=Well-cond',
948 $
'itioned, clustered e.vals. ',
' 16=Ill-cond., random comp',
949 $
'lex ', /
' 12=Well-cond., random complex ', 6x,
' ',
950 $
' 17=Ill-cond., large rand. complx ', /
' 13=Ill-condi',
951 $
'tioned, evenly spaced. ',
' 18=Ill-cond., small rand.',
953 9996
FORMAT(
' 19=Matrix with random O(1) entries. ',
' 21=Matrix ',
954 $
'with small random entries.', /
' 20=Matrix with large ran',
955 $
'dom entries. ', / )
956 9995
FORMAT(
' Tests performed with test threshold =', f8.2,
957 $ / /
' 1 = | A VR - VR W | / ( n |A| ulp ) ',
958 $ /
' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ',
959 $ /
' 3 = | |VR(i)| - 1 | / ulp ',
960 $ /
' 4 = | |VL(i)| - 1 | / ulp ',
961 $ /
' 5 = 0 if W same no matter if VR or VL computed,',
962 $
' 1/ulp otherwise', /
963 $
' 6 = 0 if VR same no matter if VL computed,',
964 $
' 1/ulp otherwise', /
965 $
' 7 = 0 if VL same no matter if VR computed,',
966 $
' 1/ulp otherwise', / )
967 9994
FORMAT(
' N=', i5,
', IWK=', i2,
', seed=', 4( i4,
',' ),
968 $
' type ', i2,
', test(', i2,
')=', g10.3 )
969 9993
FORMAT(
' DDRVEV: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
970 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )