385 SUBROUTINE ddrves( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
386 $ NOUNIT, A, LDA, H, HT, WR, WI, WRT, WIT, VS,
387 $ LDVS, RESULT, WORK, NWORK, IWORK, BWORK, INFO )
394 INTEGER INFO, LDA, LDVS, NOUNIT, NSIZES, NTYPES, NWORK
395 DOUBLE PRECISION THRESH
398 LOGICAL BWORK( * ), DOTYPE( * )
399 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
400 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), HT( LDA, * ),
401 $ result( 13 ), vs( ldvs, * ), wi( * ), wit( * ),
402 $ work( * ), wr( * ), wrt( * )
408 DOUBLE PRECISION ZERO, ONE
409 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
411 parameter( maxtyp = 21 )
417 INTEGER I, IINFO, IMODE, ISORT, ITYPE, IWK, J, JCOL,
418 $ jsize, jtype, knteig, lwork, mtypes, n, nerrs,
419 $ nfail, nmax, nnwork, ntest, ntestf, ntestt,
421 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RTULP, RTULPI, TMP,
425 CHARACTER ADUMMA( 1 )
426 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
427 $ kmagn( maxtyp ), kmode( maxtyp ),
429 DOUBLE PRECISION RES( 2 )
433 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
436 INTEGER SELDIM, SELOPT
439 COMMON / sslct / selopt, seldim, selval, selwr, selwi
443 DOUBLE PRECISION DLAMCH
444 EXTERNAL dslect, dlamch
451 INTRINSIC abs, max, sign, sqrt
454 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
455 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
457 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
458 $ 1, 5, 5, 5, 4, 3, 1 /
459 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
463 path( 1: 1 ) =
'Double precision'
478 nmax = max( nmax, nn( j ) )
485 IF( nsizes.LT.0 )
THEN
487 ELSE IF( badnn )
THEN
489 ELSE IF( ntypes.LT.0 )
THEN
491 ELSE IF( thresh.LT.zero )
THEN
493 ELSE IF( nounit.LE.0 )
THEN
495 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
497 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax )
THEN
499 ELSE IF( 5*nmax+2*nmax**2.GT.nwork )
THEN
504 CALL xerbla(
'DDRVES', -info )
510 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
515 unfl = dlamch(
'Safe minimum' )
517 ulp = dlamch(
'Precision' )
526 DO 270 jsize = 1, nsizes
529 IF( nsizes.EQ.1 .AND. ntypes.EQ.maxtyp+1 )
530 $ mtypes = mtypes + 1
532 DO 260 jtype = 1, mtypes
533 IF( .NOT.dotype( jtype ) )
539 ioldsd( j ) = iseed( j )
558 IF( mtypes.GT.maxtyp )
561 itype = ktype( jtype )
562 imode = kmode( jtype )
566 GO TO ( 30, 40, 50 )kmagn( jtype )
582 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
590 IF( itype.EQ.1 )
THEN
593 ELSE IF( itype.EQ.2 )
THEN
598 a( jcol, jcol ) = anorm
601 ELSE IF( itype.EQ.3 )
THEN
606 a( jcol, jcol ) = anorm
608 $ a( jcol, jcol-1 ) = one
611 ELSE IF( itype.EQ.4 )
THEN
615 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
616 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
619 ELSE IF( itype.EQ.5 )
THEN
623 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
624 $ anorm, n, n,
'N', a, lda, work( n+1 ),
627 ELSE IF( itype.EQ.6 )
THEN
631 IF( kconds( jtype ).EQ.1 )
THEN
633 ELSE IF( kconds( jtype ).EQ.2 )
THEN
640 CALL dlatme( n,
'S', iseed, work, imode, cond, one,
641 $ adumma,
'T',
'T',
'T', work( n+1 ), 4,
642 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
645 ELSE IF( itype.EQ.7 )
THEN
649 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
650 $
'T',
'N', work( n+1 ), 1, one,
651 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
652 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
654 ELSE IF( itype.EQ.8 )
THEN
658 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
659 $
'T',
'N', work( n+1 ), 1, one,
660 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
661 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
663 ELSE IF( itype.EQ.9 )
THEN
667 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
668 $
'T',
'N', work( n+1 ), 1, one,
669 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
670 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
672 CALL dlaset(
'Full', 2, n, zero, zero, a, lda )
673 CALL dlaset(
'Full', n-3, 1, zero, zero, a( 3, 1 ),
675 CALL dlaset(
'Full', n-3, 2, zero, zero, a( 3, n-1 ),
677 CALL dlaset(
'Full', 1, n, zero, zero, a( n, 1 ),
681 ELSE IF( itype.EQ.10 )
THEN
685 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
686 $
'T',
'N', work( n+1 ), 1, one,
687 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
688 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
695 IF( iinfo.NE.0 )
THEN
696 WRITE( nounit, fmt = 9992 )
'Generator', iinfo, n, jtype,
710 nnwork = 5*n + 2*n**2
712 nnwork = max( nnwork, 1 )
723 IF( isort.EQ.0 )
THEN
733 CALL dlacpy(
'F', n, n, a, lda, h, lda )
734 CALL dgees(
'V', sort, dslect, n, h, lda, sdim, wr,
735 $ wi, vs, ldvs, work, nnwork, bwork, iinfo )
736 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
737 result( 1+rsub ) = ulpinv
738 WRITE( nounit, fmt = 9992 )
'DGEES1', iinfo, n,
746 result( 1+rsub ) = zero
749 IF( h( i, j ).NE.zero )
750 $ result( 1+rsub ) = ulpinv
754 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.
755 $ zero )result( 1+rsub ) = ulpinv
758 IF( h( i+1, i ).NE.zero )
THEN
759 IF( h( i, i ).NE.h( i+1, i+1 ) .OR.
760 $ h( i, i+1 ).EQ.zero .OR.
761 $ sign( one, h( i+1, i ) ).EQ.
762 $ sign( one, h( i, i+1 ) ) )result( 1+rsub )
769 lwork = max( 1, 2*n*n )
770 CALL dhst01( n, 1, n, a, lda, h, lda, vs, ldvs, work,
772 result( 2+rsub ) = res( 1 )
773 result( 3+rsub ) = res( 2 )
777 result( 4+rsub ) = zero
779 IF( h( i, i ).NE.wr( i ) )
780 $ result( 4+rsub ) = ulpinv
783 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
784 $ result( 4+rsub ) = ulpinv
785 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
786 $ result( 4+rsub ) = ulpinv
789 IF( h( i+1, i ).NE.zero )
THEN
790 tmp = sqrt( abs( h( i+1, i ) ) )*
791 $ sqrt( abs( h( i, i+1 ) ) )
792 result( 4+rsub ) = max( result( 4+rsub ),
793 $ abs( wi( i )-tmp ) /
794 $ max( ulp*tmp, unfl ) )
795 result( 4+rsub ) = max( result( 4+rsub ),
796 $ abs( wi( i+1 )+tmp ) /
797 $ max( ulp*tmp, unfl ) )
798 ELSE IF( i.GT.1 )
THEN
799 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.
800 $ zero .AND. wi( i ).NE.zero )result( 4+rsub )
807 CALL dlacpy(
'F', n, n, a, lda, ht, lda )
808 CALL dgees(
'N', sort, dslect, n, ht, lda, sdim, wrt,
809 $ wit, vs, ldvs, work, nnwork, bwork,
811 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
812 result( 5+rsub ) = ulpinv
813 WRITE( nounit, fmt = 9992 )
'DGEES2', iinfo, n,
819 result( 5+rsub ) = zero
822 IF( h( i, j ).NE.ht( i, j ) )
823 $ result( 5+rsub ) = ulpinv
829 result( 6+rsub ) = zero
831 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
832 $ result( 6+rsub ) = ulpinv
837 IF( isort.EQ.1 )
THEN
841 IF( dslect( wr( i ), wi( i ) ) .OR.
842 $ dslect( wr( i ), -wi( i ) ) )
843 $ knteig = knteig + 1
845 IF( ( dslect( wr( i+1 ),
846 $ wi( i+1 ) ) .OR. dslect( wr( i+1 ),
847 $ -wi( i+1 ) ) ) .AND.
848 $ ( .NOT.( dslect( wr( i ),
849 $ wi( i ) ) .OR. dslect( wr( i ),
850 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )
851 $ result( 13 ) = ulpinv
854 IF( sdim.NE.knteig )
THEN
855 result( 13 ) = ulpinv
868 IF( result( j ).GE.zero )
870 IF( result( j ).GE.thresh )
875 $ ntestf = ntestf + 1
876 IF( ntestf.EQ.1 )
THEN
877 WRITE( nounit, fmt = 9999 )path
878 WRITE( nounit, fmt = 9998 )
879 WRITE( nounit, fmt = 9997 )
880 WRITE( nounit, fmt = 9996 )
881 WRITE( nounit, fmt = 9995 )thresh
882 WRITE( nounit, fmt = 9994 )
887 IF( result( j ).GE.thresh )
THEN
888 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
893 nerrs = nerrs + nfail
894 ntestt = ntestt + ntest
902 CALL dlasum( path, nounit, nerrs, ntestt )
904 9999
FORMAT( / 1x, a3,
' -- Real Schur Form Decomposition Driver',
905 $ /
' Matrix types (see DDRVES for details): ' )
907 9998
FORMAT( /
' Special Matrices:', /
' 1=Zero matrix. ',
908 $
' ',
' 5=Diagonal: geometr. spaced entries.',
909 $ /
' 2=Identity matrix. ',
' 6=Diagona',
910 $
'l: clustered entries.', /
' 3=Transposed Jordan block. ',
911 $
' ',
' 7=Diagonal: large, evenly spaced.', /
' ',
912 $
'4=Diagonal: evenly spaced entries. ',
' 8=Diagonal: s',
913 $
'mall, evenly spaced.' )
914 9997
FORMAT(
' Dense, Non-Symmetric Matrices:', /
' 9=Well-cond., ev',
915 $
'enly spaced eigenvals.',
' 14=Ill-cond., geomet. spaced e',
916 $
'igenals.', /
' 10=Well-cond., geom. spaced eigenvals. ',
917 $
' 15=Ill-conditioned, clustered e.vals.', /
' 11=Well-cond',
918 $
'itioned, clustered e.vals. ',
' 16=Ill-cond., random comp',
919 $
'lex ', /
' 12=Well-cond., random complex ', 6x,
' ',
920 $
' 17=Ill-cond., large rand. complx ', /
' 13=Ill-condi',
921 $
'tioned, evenly spaced. ',
' 18=Ill-cond., small rand.',
923 9996
FORMAT(
' 19=Matrix with random O(1) entries. ',
' 21=Matrix ',
924 $
'with small random entries.', /
' 20=Matrix with large ran',
925 $
'dom entries. ', / )
926 9995
FORMAT(
' Tests performed with test threshold =', f8.2,
927 $ /
' ( A denotes A on input and T denotes A on output)',
928 $ / /
' 1 = 0 if T in Schur form (no sort), ',
929 $
' 1/ulp otherwise', /
930 $
' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
931 $ /
' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
932 $
' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
933 $
' 1/ulp otherwise', /
934 $
' 5 = 0 if T same no matter if VS computed (no sort),',
935 $
' 1/ulp otherwise', /
936 $
' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
937 $
', 1/ulp otherwise' )
938 9994
FORMAT(
' 7 = 0 if T in Schur form (sort), ',
' 1/ulp otherwise',
939 $ /
' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
940 $ /
' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
941 $ /
' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
942 $
' 1/ulp otherwise', /
943 $
' 11 = 0 if T same no matter if VS computed (sort),',
944 $
' 1/ulp otherwise', /
945 $
' 12 = 0 if WR, WI same no matter if VS computed (sort),',
946 $
' 1/ulp otherwise', /
947 $
' 13 = 0 if sorting successful, 1/ulp otherwise', / )
948 9993
FORMAT(
' N=', i5,
', IWK=', i2,
', seed=', 4( i4,
',' ),
949 $
' type ', i2,
', test(', i2,
')=', g10.3 )
950 9992
FORMAT(
' DDRVES: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
951 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )