387 SUBROUTINE ddrves( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
388 $ nounit, a, lda, h, ht, wr, wi, wrt, wit, vs,
389 $ ldvs, result, work, nwork, iwork, bwork, info )
397 INTEGER info, lda, ldvs, nounit, nsizes, ntypes, nwork
398 DOUBLE PRECISION thresh
401 LOGICAL bwork( * ), dotype( * )
402 INTEGER iseed( 4 ), iwork( * ), nn( * )
403 DOUBLE PRECISION a( lda, * ), h( lda, * ), ht( lda, * ),
404 $ result( 13 ), vs( ldvs, * ), wi( * ), wit( * ),
405 $ work( * ), wr( * ), wrt( * )
411 DOUBLE PRECISION zero, one
412 parameter( zero = 0.0d0, one = 1.0d0 )
414 parameter( maxtyp = 21 )
420 INTEGER i, iinfo, imode, isort, itype, iwk, j, jcol,
421 $ jsize, jtype, knteig, lwork, mtypes, n, nerrs,
422 $ nfail, nmax, nnwork, ntest, ntestf, ntestt,
424 DOUBLE PRECISION anorm, cond, conds, ovfl, rtulp, rtulpi, tmp,
428 CHARACTER adumma( 1 )
429 INTEGER idumma( 1 ), ioldsd( 4 ), kconds( maxtyp ),
430 $ kmagn( maxtyp ), kmode( maxtyp ),
432 DOUBLE PRECISION res( 2 )
436 DOUBLE PRECISION selwi( 20 ), selwr( 20 )
439 INTEGER seldim, selopt
442 common / sslct / selopt, seldim, selval, selwr, selwi
454 INTRINSIC abs, max, sign, sqrt
457 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
458 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
460 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
461 $ 1, 5, 5, 5, 4, 3, 1 /
462 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
466 path( 1: 1 ) =
'Double precision'
481 nmax = max( nmax, nn( j ) )
488 IF( nsizes.LT.0 )
THEN
490 ELSE IF( badnn )
THEN
492 ELSE IF( ntypes.LT.0 )
THEN
494 ELSE IF( thresh.LT.zero )
THEN
496 ELSE IF( nounit.LE.0 )
THEN
498 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
500 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax )
THEN
502 ELSE IF( 5*nmax+2*nmax**2.GT.nwork )
THEN
507 CALL
xerbla(
'DDRVES', -info )
513 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
518 unfl =
dlamch(
'Safe minimum' )
521 ulp =
dlamch(
'Precision' )
530 DO 270 jsize = 1, nsizes
533 IF( nsizes.EQ.1 .AND. ntypes.EQ.maxtyp+1 )
534 $ mtypes = mtypes + 1
536 DO 260 jtype = 1, mtypes
537 IF( .NOT.dotype( jtype ) )
543 ioldsd( j ) = iseed( j )
562 IF( mtypes.GT.maxtyp )
565 itype = ktype( jtype )
566 imode = kmode( jtype )
570 go to( 30, 40, 50 )kmagn( jtype )
586 CALL
dlaset(
'Full', lda, n, zero, zero, a, lda )
594 IF( itype.EQ.1 )
THEN
597 ELSE IF( itype.EQ.2 )
THEN
602 a( jcol, jcol ) = anorm
605 ELSE IF( itype.EQ.3 )
THEN
610 a( jcol, jcol ) = anorm
612 $ a( jcol, jcol-1 ) = one
615 ELSE IF( itype.EQ.4 )
THEN
619 CALL
dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
620 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
623 ELSE IF( itype.EQ.5 )
THEN
627 CALL
dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
628 $ anorm, n, n,
'N', a, lda, work( n+1 ),
631 ELSE IF( itype.EQ.6 )
THEN
635 IF( kconds( jtype ).EQ.1 )
THEN
637 ELSE IF( kconds( jtype ).EQ.2 )
THEN
644 CALL
dlatme( n,
'S', iseed, work, imode, cond, one,
645 $ adumma,
'T',
'T',
'T', work( n+1 ), 4,
646 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
649 ELSE IF( itype.EQ.7 )
THEN
653 CALL
dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
654 $
'T',
'N', work( n+1 ), 1, one,
655 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
656 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
658 ELSE IF( itype.EQ.8 )
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, n, n,
665 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
667 ELSE IF( itype.EQ.9 )
THEN
671 CALL
dlatmr( n, n,
'S', iseed,
'N', 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 CALL
dlaset(
'Full', 2, n, zero, zero, a, lda )
677 CALL
dlaset(
'Full', n-3, 1, zero, zero, a( 3, 1 ),
679 CALL
dlaset(
'Full', n-3, 2, zero, zero, a( 3, n-1 ),
681 CALL
dlaset(
'Full', 1, n, zero, zero, a( n, 1 ),
685 ELSE IF( itype.EQ.10 )
THEN
689 CALL
dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
690 $
'T',
'N', work( n+1 ), 1, one,
691 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
692 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
699 IF( iinfo.NE.0 )
THEN
700 WRITE( nounit, fmt = 9992 )
'Generator', iinfo, n, jtype,
714 nnwork = 5*n + 2*n**2
716 nnwork = max( nnwork, 1 )
727 IF( isort.EQ.0 )
THEN
737 CALL
dlacpy(
'F', n, n, a, lda, h, lda )
738 CALL
dgees(
'V', sort,
dslect, n, h, lda, sdim, wr,
739 $ wi, vs, ldvs, work, nnwork, bwork, iinfo )
740 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
741 result( 1+rsub ) = ulpinv
742 WRITE( nounit, fmt = 9992 )
'DGEES1', iinfo, n,
750 result( 1+rsub ) = zero
753 IF( h( i, j ).NE.zero )
754 $ result( 1+rsub ) = ulpinv
758 IF( h( i+1, i ).NE.zero .AND. h( i+2, i+1 ).NE.
759 $ zero )result( 1+rsub ) = ulpinv
762 IF( h( i+1, i ).NE.zero )
THEN
763 IF( h( i, i ).NE.h( i+1, i+1 ) .OR.
764 $ h( i, i+1 ).EQ.zero .OR.
765 $ sign( one, h( i+1, i ) ).EQ.
766 $ sign( one, h( i, i+1 ) ) )result( 1+rsub )
773 lwork = max( 1, 2*n*n )
774 CALL
dhst01( n, 1, n, a, lda, h, lda, vs, ldvs, work,
776 result( 2+rsub ) = res( 1 )
777 result( 3+rsub ) = res( 2 )
781 result( 4+rsub ) = zero
783 IF( h( i, i ).NE.wr( i ) )
784 $ result( 4+rsub ) = ulpinv
787 IF( h( 2, 1 ).EQ.zero .AND. wi( 1 ).NE.zero )
788 $ result( 4+rsub ) = ulpinv
789 IF( h( n, n-1 ).EQ.zero .AND. wi( n ).NE.zero )
790 $ result( 4+rsub ) = ulpinv
793 IF( h( i+1, i ).NE.zero )
THEN
794 tmp = sqrt( abs( h( i+1, i ) ) )*
795 $ sqrt( abs( h( i, i+1 ) ) )
796 result( 4+rsub ) = max( result( 4+rsub ),
797 $ abs( wi( i )-tmp ) /
798 $ max( ulp*tmp, unfl ) )
799 result( 4+rsub ) = max( result( 4+rsub ),
800 $ abs( wi( i+1 )+tmp ) /
801 $ max( ulp*tmp, unfl ) )
802 ELSE IF( i.GT.1 )
THEN
803 IF( h( i+1, i ).EQ.zero .AND. h( i, i-1 ).EQ.
804 $ zero .AND. wi( i ).NE.zero )result( 4+rsub )
811 CALL
dlacpy(
'F', n, n, a, lda, ht, lda )
812 CALL
dgees(
'N', sort,
dslect, n, ht, lda, sdim, wrt,
813 $ wit, vs, ldvs, work, nnwork, bwork,
815 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
816 result( 5+rsub ) = ulpinv
817 WRITE( nounit, fmt = 9992 )
'DGEES2', iinfo, n,
823 result( 5+rsub ) = zero
826 IF( h( i, j ).NE.ht( i, j ) )
827 $ result( 5+rsub ) = ulpinv
833 result( 6+rsub ) = zero
835 IF( wr( i ).NE.wrt( i ) .OR. wi( i ).NE.wit( i ) )
836 $ result( 6+rsub ) = ulpinv
841 IF( isort.EQ.1 )
THEN
845 IF(
dslect( wr( i ), wi( i ) ) .OR.
846 $
dslect( wr( i ), -wi( i ) ) )
847 $ knteig = knteig + 1
850 $ wi( i+1 ) ) .OR.
dslect( wr( i+1 ),
851 $ -wi( i+1 ) ) ) .AND.
852 $ ( .NOT.(
dslect( wr( i ),
853 $ wi( i ) ) .OR.
dslect( wr( i ),
854 $ -wi( i ) ) ) ) .AND. iinfo.NE.n+2 )
855 $ result( 13 ) = ulpinv
858 IF( sdim.NE.knteig )
THEN
859 result( 13 ) = ulpinv
872 IF( result( j ).GE.zero )
874 IF( result( j ).GE.thresh )
879 $ ntestf = ntestf + 1
880 IF( ntestf.EQ.1 )
THEN
881 WRITE( nounit, fmt = 9999 )path
882 WRITE( nounit, fmt = 9998 )
883 WRITE( nounit, fmt = 9997 )
884 WRITE( nounit, fmt = 9996 )
885 WRITE( nounit, fmt = 9995 )thresh
886 WRITE( nounit, fmt = 9994 )
891 IF( result( j ).GE.thresh )
THEN
892 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
897 nerrs = nerrs + nfail
898 ntestt = ntestt + ntest
906 CALL
dlasum( path, nounit, nerrs, ntestt )
908 9999 format( / 1x, a3,
' -- Real Schur Form Decomposition Driver',
909 $ /
' Matrix types (see DDRVES for details): ' )
911 9998 format( /
' Special Matrices:', /
' 1=Zero matrix. ',
912 $
' ',
' 5=Diagonal: geometr. spaced entries.',
913 $ /
' 2=Identity matrix. ',
' 6=Diagona',
914 $
'l: clustered entries.', /
' 3=Transposed Jordan block. ',
915 $
' ',
' 7=Diagonal: large, evenly spaced.', /
' ',
916 $
'4=Diagonal: evenly spaced entries. ',
' 8=Diagonal: s',
917 $
'mall, evenly spaced.' )
918 9997 format(
' Dense, Non-Symmetric Matrices:', /
' 9=Well-cond., ev',
919 $
'enly spaced eigenvals.',
' 14=Ill-cond., geomet. spaced e',
920 $
'igenals.', /
' 10=Well-cond., geom. spaced eigenvals. ',
921 $
' 15=Ill-conditioned, clustered e.vals.', /
' 11=Well-cond',
922 $
'itioned, clustered e.vals. ',
' 16=Ill-cond., random comp',
923 $
'lex ', /
' 12=Well-cond., random complex ', 6x,
' ',
924 $
' 17=Ill-cond., large rand. complx ', /
' 13=Ill-condi',
925 $
'tioned, evenly spaced. ',
' 18=Ill-cond., small rand.',
927 9996 format(
' 19=Matrix with random O(1) entries. ',
' 21=Matrix ',
928 $
'with small random entries.', /
' 20=Matrix with large ran',
929 $
'dom entries. ', / )
930 9995 format(
' Tests performed with test threshold =', f8.2,
931 $ /
' ( A denotes A on input and T denotes A on output)',
932 $ / /
' 1 = 0 if T in Schur form (no sort), ',
933 $
' 1/ulp otherwise', /
934 $
' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
935 $ /
' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
936 $
' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
937 $
' 1/ulp otherwise', /
938 $
' 5 = 0 if T same no matter if VS computed (no sort),',
939 $
' 1/ulp otherwise', /
940 $
' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
941 $
', 1/ulp otherwise' )
942 9994 format(
' 7 = 0 if T in Schur form (sort), ',
' 1/ulp otherwise',
943 $ /
' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
944 $ /
' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
945 $ /
' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
946 $
' 1/ulp otherwise', /
947 $
' 11 = 0 if T same no matter if VS computed (sort),',
948 $
' 1/ulp otherwise', /
949 $
' 12 = 0 if WR, WI same no matter if VS computed (sort),',
950 $
' 1/ulp otherwise', /
951 $
' 13 = 0 if sorting succesful, 1/ulp otherwise', / )
952 9993 format(
' N=', i5,
', IWK=', i2,
', seed=', 4( i4,
',' ),
953 $
' type ', i2,
', test(', i2,
')=', g10.3 )
954 9992 format(
' DDRVES: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
955 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )