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,
')' )
subroutine xerbla(srname, info)
subroutine ddrves(nsizes, nn, ntypes, dotype, iseed, thresh, nounit, a, lda, h, ht, wr, wi, wrt, wit, vs, ldvs, result, work, nwork, iwork, bwork, info)
DDRVES
subroutine dhst01(n, ilo, ihi, a, lda, h, ldh, q, ldq, work, lwork, result)
DHST01
subroutine dlasum(type, iounit, ie, nrun)
DLASUM
subroutine dlatme(n, dist, iseed, d, mode, cond, dmax, ei, rsign, upper, sim, ds, modes, conds, kl, ku, anorm, a, lda, work, info)
DLATME
subroutine dlatmr(m, n, dist, iseed, sym, d, mode, cond, dmax, rsign, grade, dl, model, condl, dr, moder, condr, pivtng, ipivot, kl, ku, sparse, anorm, pack, a, lda, iwork, info)
DLATMR
subroutine dlatms(m, n, dist, iseed, sym, d, mode, cond, dmax, kl, ku, pack, a, lda, work, info)
DLATMS
subroutine dgees(jobvs, sort, select, n, a, lda, sdim, wr, wi, vs, ldvs, work, lwork, bwork, info)
DGEES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE m...
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.