452 SUBROUTINE ddrvsx( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
453 $ niunit, nounit, a, lda, h, ht, wr, wi, wrt,
454 $ wit, wrtmp, witmp, vs, ldvs, vs1, result, work,
455 $ lwork, iwork, bwork, info )
463 INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
465 DOUBLE PRECISION THRESH
468 LOGICAL BWORK( * ), DOTYPE( * )
469 INTEGER ISEED( 4 ), IWORK( * ), NN( * )
470 DOUBLE PRECISION A( lda, * ), H( lda, * ), HT( lda, * ),
471 $ result( 17 ), vs( ldvs, * ), vs1( ldvs, * ),
472 $ wi( * ), wit( * ), witmp( * ), work( * ),
473 $ wr( * ), wrt( * ), wrtmp( * )
479 DOUBLE PRECISION ZERO, ONE
480 parameter ( zero = 0.0d0, one = 1.0d0 )
482 parameter ( maxtyp = 21 )
487 INTEGER I, IINFO, IMODE, ITYPE, IWK, J, JCOL, JSIZE,
488 $ jtype, mtypes, n, nerrs, nfail, nmax, nnwork,
489 $ nslct, ntest, ntestf, ntestt
490 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
491 $ rtulp, rtulpi, ulp, ulpinv, unfl
494 CHARACTER ADUMMA( 1 )
495 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ISLCT( 20 ),
496 $ kconds( maxtyp ), kmagn( maxtyp ),
497 $ kmode( maxtyp ), ktype( maxtyp )
501 DOUBLE PRECISION SELWI( 20 ), SELWR( 20 )
504 INTEGER SELDIM, SELOPT
507 COMMON / sslct / selopt, seldim, selval, selwr, selwi
510 DOUBLE PRECISION DLAMCH
518 INTRINSIC abs, max, min, sqrt
521 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
522 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
524 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
525 $ 1, 5, 5, 5, 4, 3, 1 /
526 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
530 path( 1: 1 ) =
'Double precision'
548 nmax = max( nmax, nn( j ) )
555 IF( nsizes.LT.0 )
THEN
557 ELSE IF( badnn )
THEN
559 ELSE IF( ntypes.LT.0 )
THEN
561 ELSE IF( thresh.LT.zero )
THEN
563 ELSE IF( niunit.LE.0 )
THEN
565 ELSE IF( nounit.LE.0 )
THEN
567 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
569 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax )
THEN
571 ELSE IF( max( 3*nmax, 2*nmax**2 ).GT.lwork )
THEN
576 CALL xerbla(
'DDRVSX', -info )
582 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
587 unfl = dlamch(
'Safe minimum' )
590 ulp = dlamch(
'Precision' )
599 DO 140 jsize = 1, nsizes
601 IF( nsizes.NE.1 )
THEN
602 mtypes = min( maxtyp, ntypes )
604 mtypes = min( maxtyp+1, ntypes )
607 DO 130 jtype = 1, mtypes
608 IF( .NOT.dotype( jtype ) )
614 ioldsd( j ) = iseed( j )
633 IF( mtypes.GT.maxtyp )
636 itype = ktype( jtype )
637 imode = kmode( jtype )
641 GO TO ( 30, 40, 50 )kmagn( jtype )
657 CALL dlaset(
'Full', lda, n, zero, zero, a, lda )
665 IF( itype.EQ.1 )
THEN
668 ELSE IF( itype.EQ.2 )
THEN
673 a( jcol, jcol ) = anorm
676 ELSE IF( itype.EQ.3 )
THEN
681 a( jcol, jcol ) = anorm
683 $ a( jcol, jcol-1 ) = one
686 ELSE IF( itype.EQ.4 )
THEN
690 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
691 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
694 ELSE IF( itype.EQ.5 )
THEN
698 CALL dlatms( n, n,
'S', iseed,
'S', work, imode, cond,
699 $ anorm, n, n,
'N', a, lda, work( n+1 ),
702 ELSE IF( itype.EQ.6 )
THEN
706 IF( kconds( jtype ).EQ.1 )
THEN
708 ELSE IF( kconds( jtype ).EQ.2 )
THEN
715 CALL dlatme( n,
'S', iseed, work, imode, cond, one,
716 $ adumma,
'T',
'T',
'T', work( n+1 ), 4,
717 $ conds, n, n, anorm, a, lda, work( 2*n+1 ),
720 ELSE IF( itype.EQ.7 )
THEN
724 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
725 $
'T',
'N', work( n+1 ), 1, one,
726 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
727 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
729 ELSE IF( itype.EQ.8 )
THEN
733 CALL dlatmr( n, n,
'S', iseed,
'S', work, 6, one, one,
734 $
'T',
'N', work( n+1 ), 1, one,
735 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
736 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
738 ELSE IF( itype.EQ.9 )
THEN
742 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
743 $
'T',
'N', work( n+1 ), 1, one,
744 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
745 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
747 CALL dlaset(
'Full', 2, n, zero, zero, a, lda )
748 CALL dlaset(
'Full', n-3, 1, zero, zero, a( 3, 1 ),
750 CALL dlaset(
'Full', n-3, 2, zero, zero, a( 3, n-1 ),
752 CALL dlaset(
'Full', 1, n, zero, zero, a( n, 1 ),
756 ELSE IF( itype.EQ.10 )
THEN
760 CALL dlatmr( n, n,
'S', iseed,
'N', work, 6, one, one,
761 $
'T',
'N', work( n+1 ), 1, one,
762 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
763 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
770 IF( iinfo.NE.0 )
THEN
771 WRITE( nounit, fmt = 9991 )
'Generator', iinfo, n, jtype,
785 nnwork = max( 3*n, 2*n*n )
787 nnwork = max( nnwork, 1 )
789 CALL dget24( .false., jtype, thresh, ioldsd, nounit, n,
790 $ a, lda, h, ht, wr, wi, wrt, wit, wrtmp,
791 $ witmp, vs, ldvs, vs1, rcdein, rcdvin, nslct,
792 $ islct, result, work, nnwork, iwork, bwork,
800 IF( result( j ).GE.zero )
802 IF( result( j ).GE.thresh )
807 $ ntestf = ntestf + 1
808 IF( ntestf.EQ.1 )
THEN
809 WRITE( nounit, fmt = 9999 )path
810 WRITE( nounit, fmt = 9998 )
811 WRITE( nounit, fmt = 9997 )
812 WRITE( nounit, fmt = 9996 )
813 WRITE( nounit, fmt = 9995 )thresh
814 WRITE( nounit, fmt = 9994 )
819 IF( result( j ).GE.thresh )
THEN
820 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
825 nerrs = nerrs + nfail
826 ntestt = ntestt + ntest
839 READ( niunit, fmt = *, end = 200 )n, nslct
845 $
READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
847 READ( niunit, fmt = * )( a( i, j ), j = 1, n )
849 READ( niunit, fmt = * )rcdein, rcdvin
851 CALL dget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
852 $ wr, wi, wrt, wit, wrtmp, witmp, vs, ldvs, vs1,
853 $ rcdein, rcdvin, nslct, islct, result, work, lwork,
854 $ iwork, bwork, info )
861 IF( result( j ).GE.zero )
863 IF( result( j ).GE.thresh )
868 $ ntestf = ntestf + 1
869 IF( ntestf.EQ.1 )
THEN
870 WRITE( nounit, fmt = 9999 )path
871 WRITE( nounit, fmt = 9998 )
872 WRITE( nounit, fmt = 9997 )
873 WRITE( nounit, fmt = 9996 )
874 WRITE( nounit, fmt = 9995 )thresh
875 WRITE( nounit, fmt = 9994 )
879 IF( result( j ).GE.thresh )
THEN
880 WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
884 nerrs = nerrs + nfail
885 ntestt = ntestt + ntest
891 CALL dlasum( path, nounit, nerrs, ntestt )
893 9999
FORMAT( / 1x, a3,
' -- Real Schur Form Decomposition Expert ',
894 $
'Driver', /
' Matrix types (see DDRVSX for details):' )
896 9998
FORMAT( /
' Special Matrices:', /
' 1=Zero matrix. ',
897 $
' ',
' 5=Diagonal: geometr. spaced entries.',
898 $ /
' 2=Identity matrix. ',
' 6=Diagona',
899 $
'l: clustered entries.', /
' 3=Transposed Jordan block. ',
900 $
' ',
' 7=Diagonal: large, evenly spaced.', /
' ',
901 $
'4=Diagonal: evenly spaced entries. ',
' 8=Diagonal: s',
902 $
'mall, evenly spaced.' )
903 9997
FORMAT(
' Dense, Non-Symmetric Matrices:', /
' 9=Well-cond., ev',
904 $
'enly spaced eigenvals.',
' 14=Ill-cond., geomet. spaced e',
905 $
'igenals.', /
' 10=Well-cond., geom. spaced eigenvals. ',
906 $
' 15=Ill-conditioned, clustered e.vals.', /
' 11=Well-cond',
907 $
'itioned, clustered e.vals. ',
' 16=Ill-cond., random comp',
908 $
'lex ', /
' 12=Well-cond., random complex ',
' ',
909 $
' 17=Ill-cond., large rand. complx ', /
' 13=Ill-condi',
910 $
'tioned, evenly spaced. ',
' 18=Ill-cond., small rand.',
912 9996
FORMAT(
' 19=Matrix with random O(1) entries. ',
' 21=Matrix ',
913 $
'with small random entries.', /
' 20=Matrix with large ran',
914 $
'dom entries. ', / )
915 9995
FORMAT(
' Tests performed with test threshold =', f8.2,
916 $ /
' ( A denotes A on input and T denotes A on output)',
917 $ / /
' 1 = 0 if T in Schur form (no sort), ',
918 $
' 1/ulp otherwise', /
919 $
' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
920 $ /
' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ', /
921 $
' 4 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (no sort),',
922 $
' 1/ulp otherwise', /
923 $
' 5 = 0 if T same no matter if VS computed (no sort),',
924 $
' 1/ulp otherwise', /
925 $
' 6 = 0 if WR, WI same no matter if VS computed (no sort)',
926 $
', 1/ulp otherwise' )
927 9994
FORMAT(
' 7 = 0 if T in Schur form (sort), ',
' 1/ulp otherwise',
928 $ /
' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
929 $ /
' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
930 $ /
' 10 = 0 if WR+sqrt(-1)*WI are eigenvalues of T (sort),',
931 $
' 1/ulp otherwise', /
932 $
' 11 = 0 if T same no matter what else computed (sort),',
933 $
' 1/ulp otherwise', /
934 $
' 12 = 0 if WR, WI same no matter what else computed ',
935 $
'(sort), 1/ulp otherwise', /
936 $
' 13 = 0 if sorting successful, 1/ulp otherwise',
937 $ /
' 14 = 0 if RCONDE same no matter what else computed,',
938 $
' 1/ulp otherwise', /
939 $
' 15 = 0 if RCONDv same no matter what else computed,',
940 $
' 1/ulp otherwise', /
941 $
' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
942 $ /
' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
943 9993
FORMAT(
' N=', i5,
', IWK=', i2,
', seed=', 4( i4,
',' ),
944 $
' type ', i2,
', test(', i2,
')=', g10.3 )
945 9992
FORMAT(
' N=', i5,
', input example =', i3,
', test(', i2,
')=',
947 9991
FORMAT(
' DDRVSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
948 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
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 dlatme(N, DIST, ISEED, D, MODE, COND, DMAX, EI, RSIGN, UPPER, SIM, DS, MODES, CONDS, KL, KU, ANORM, A, LDA, WORK, INFO)
DLATME
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...
subroutine dget24(COMP, JTYPE, THRESH, ISEED, NOUNIT, N, A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1, RCDEIN, RCDVIN, NSLCT, ISLCT, RESULT, WORK, LWORK, IWORK, BWORK, INFO)
DGET24
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine ddrvsx(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NIUNIT, NOUNIT, A, LDA, H, HT, WR, WI, WRT, WIT, WRTMP, WITMP, VS, LDVS, VS1, RESULT, WORK, LWORK, IWORK, BWORK, INFO)
DDRVSX
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dlasum(TYPE, IOUNIT, IE, NRUN)
DLASUM
subroutine dlatms(M, N, DIST, ISEED, SYM, D, MODE, COND, DMAX, KL, KU, PACK, A, LDA, WORK, INFO)
DLATMS