452 SUBROUTINE sdrvsx( 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,
468 LOGICAL bwork( * ), dotype( * )
469 INTEGER iseed( 4 ), iwork( * ), nn( * )
470 REAL a( lda, * ), h( lda, * ), ht( lda, * ),
471 $ result( 17 ), vs( ldvs, * ), vs1( ldvs, * ),
472 $ wi( * ), wit( * ), witmp( * ), work( * ),
473 $ wr( * ), wrt( * ), wrtmp( * )
480 parameter( zero = 0.0e0, one = 1.0e0 )
482 parameter( maxtyp = 21 )
487 INTEGER i, iinfo, imode, itype, iwk, j, jcol, jsize,
488 $ jtype, mtypes, n, nerrs, nfail, nmax,
489 $ nnwork, nslct, ntest, ntestf, ntestt
490 REAL 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 REAL selwi( 20 ), selwr( 20 )
504 INTEGER seldim, selopt
507 common / sslct / selopt, seldim, selval, selwr, selwi
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 ) =
'Single 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(
'SDRVSX', -info )
582 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
587 unfl =
slamch(
'Safe minimum' )
590 ulp =
slamch(
'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
slaset(
'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
slatms( 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
slatms( 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
slatme( 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
slatmr( 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
slatmr( 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
slatmr( 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
slaset(
'Full', 2, n, zero, zero, a, lda )
748 CALL
slaset(
'Full', n-3, 1, zero, zero, a( 3, 1 ),
750 CALL
slaset(
'Full', n-3, 2, zero, zero, a( 3, n-1 ),
752 CALL
slaset(
'Full', 1, n, zero, zero, a( n, 1 ),
756 ELSE IF( itype.EQ.10 )
THEN
760 CALL
slatmr( 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
sget24( .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
sget24( .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
slasum( path, nounit, nerrs, ntestt )
893 9999 format( / 1x, a3,
' -- Real Schur Form Decomposition Expert ',
894 $
'Driver', /
' Matrix types (see SDRVSX 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 succesful, 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(
' SDRVSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
948 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )