431 SUBROUTINE cdrvsx( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
432 $ NIUNIT, NOUNIT, A, LDA, H, HT, W, WT, WTMP, VS,
433 $ LDVS, VS1, RESULT, WORK, LWORK, RWORK, BWORK,
441 INTEGER INFO, LDA, LDVS, LWORK, NIUNIT, NOUNIT, NSIZES,
446 LOGICAL BWORK( * ), DOTYPE( * )
447 INTEGER ISEED( 4 ), NN( * )
448 REAL ( 17 ), RWORK( * )
449 COMPLEX A( , * ), H( LDA, * ), HT( LDA, * ),
450 $ vs( ldvs, * ), vs1( ldvs, * ), w( * ),
451 $ work( * ), wt( * ), wtmp( * )
458 PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ) )
460 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
462 parameter( zero = 0.0e+0, one = 1.0e+0 )
464 parameter( maxtyp = 21 )
469 INTEGER I, IINFO, IMODE, ISRT, ITYPE, IWK, J, JCOL,
470 $ jsize, jtype, mtypes, n, nerrs, nfail,
471 $ nmax, nnwork, nslct, ntest, ntestf, ntestt
472 REAL ANORM, COND, CONDS, OVFL, RCDEIN, RCDVIN,
473 $ RTULP, RTULPI, ULP, ULPINV, UNFL
476 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), ( 20 ),
477 $ KCONDS( MAXTYP ), KMAGN( MAXTYP ),
478 $ KMODE( MAXTYP ), KTYPE( MAXTYP )
482 REAL SELWI( 20 ), SELWR( 20 )
485 INTEGER SELDIM, SELOPT
488 COMMON / sslct / selopt, seldim, selval, selwr, selwi
499 INTRINSIC abs,
max,
min, sqrt
502 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
503 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
505 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
506 $ 1, 5, 5, 5, 4, 3, 1 /
507 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
511 path( 1: 1 ) =
'Complex precision'
529 nmax =
max( nmax, nn( j ) )
536 IF( nsizes.LT.0 )
THEN
538 ELSE IF( badnn )
THEN
540 ELSE IF( ntypes.LT.0 )
THEN
542 ELSE IF( thresh.LT.zero )
THEN
544 ELSE IF( niunit.LE.0 )
THEN
546 ELSE IF( nounit.LE.0 )
THEN
548 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
550 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax )
THEN
552 ELSE IF(
max( 3*nmax, 2*nmax**2 ).GT.lwork )
THEN
557 CALL xerbla(
'CDRVSX', -info )
563 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
568 unfl = slamch(
'Safe minimum' )
571 ulp = slamch(
'Precision' )
580 DO 140 jsize = 1, nsizes
582 IF( nsizes.NE.1 )
THEN
583 mtypes =
min( maxtyp, ntypes )
585 mtypes =
min( maxtyp+1, ntypes )
588 DO 130 jtype = 1, mtypes
589 IF( .NOT.dotype( jtype ) )
595 ioldsd( j ) = iseed( j )
614 IF( mtypes.GT.maxtyp )
617 itype = ktype( jtype )
618 imode = kmode( jtype )
622 GO TO ( 30, 40, 50 )kmagn( jtype )
638 CALL claset(
'Full', lda, n, czero, czero, a, lda )
644 IF( itype.EQ.1 )
THEN
650 ELSE IF( itype.EQ.2 )
THEN
655 a( jcol, jcol ) = anorm
658 ELSE IF( itype.EQ.3 )
THEN
663 a( jcol, jcol ) = anorm
665 $ a( jcol, jcol-1 ) = cone
668 ELSE IF( itype.EQ.4 )
THEN
672 CALL clatms( n, n,
'S', iseed,
'H', rwork, imode, cond,
673 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
676 ELSE IF( itype.EQ.5 )
THEN
680 CALL clatms( n, n,
'S', iseed,
'H', rwork, imode, cond,
681 $ anorm, n, n,
'N', a, lda, work( n+1 ),
684 ELSE IF( itype.EQ.6 )
THEN
688 IF( kconds( jtype ).EQ.1 )
THEN
690 ELSE IF( kconds( jtype ).EQ.2 )
THEN
696 CALL clatme( n,
'D', iseed, work, imode, cond, cone,
697 $
'T',
'T',
'T', rwork, 4, conds, n, n, anorm,
698 $ a, lda, work( 2*n+1 ), iinfo )
700 ELSE IF( itype.EQ.7 )
THEN
704 CALL clatmr( n, n,
'D', iseed,
'N', work, 6, one, cone,
705 $
'T',
'N', work( n+1 ), 1, one,
706 $ work( 2*n+1 ), 1, one,
'N', idumma
707 $ zero, anorm,
'NO', a, lda, idumma, iinfo )
709 ELSE IF( itype.EQ.8 )
THEN
713 CALL clatmr( n, n,
'D', iseed,
'H', work, 6, one, cone,
714 $
'T',
'N', work( n+1 ), 1, one,
715 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
716 $ zero, anorm,
'NO', a, lda, idumma, iinfo )
718 ELSE IF( itype.EQ.9 )
THEN
722 CALL clatmr( n, n,
'D', iseed,
'N', work, 6, one, cone,
723 $
'T',
'N', work( n+1 ), 1, one,
724 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
725 $ zero, anorm,
'NO', a, lda, idumma, iinfo )
727 CALL claset(
'Full', 2, n, czero, czero, a, lda )
728 CALL claset(
'Full', n-3, 1, czero, czero, a( 3, 1 ),
730 CALL claset(
'Full', n-3, 2, czero, czero,
732 CALL claset(
'Full', 1, n, czero, czero, a( n, 1 ),
736 ELSE IF( itype.EQ.10 )
THEN
740 CALL clatmr( n, n,
'D', iseed,
'N', work, 6, one, cone,
741 $
'T',
'N', work( n+1 ), 1, one,
742 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
743 $ zero, anorm,
'NO', a, lda, idumma, iinfo )
750 IF( iinfo.NE.0 )
THEN
751 WRITE( nounit, fmt = 9991 )
'Generator', iinfo, n, jtype,
765 nnwork =
max( 2*n, n*( n+1 ) / 2 )
767 nnwork =
max( nnwork, 1 )
769 CALL cget24( .false., jtype, thresh, ioldsd, nounit, n,
770 $ a, lda, h, ht, w, wt, wtmp, vs, ldvs, vs1,
771 $ rcdein, rcdvin, nslct, islct, 0, result,
772 $ work, nnwork, rwork, bwork, info )
779 IF( result( j ).GE.zero )
781 IF( result( j ).GE.thresh )
786 $ ntestf = ntestf + 1
787 IF( ntestf.EQ.1 )
THEN
788 WRITE( nounit, fmt = 9999 )path
789 WRITE( nounit, fmt = 9998 )
790 WRITE( nounit, fmt = 9997 )
791 WRITE( nounit, fmt = 9996 )
792 WRITE( nounit, fmt = 9995 )thresh
793 WRITE( nounit, fmt = 9994 )
798 IF( result( j ).GE.thresh )
THEN
799 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype
804 nerrs = nerrs + nfail
805 ntestt = ntestt + ntest
818 READ( niunit, fmt = *,
END = 200 )N, NSLCT, isrt
823 READ( niunit, fmt = * )( islct( i ), i = 1, nslct )
825 READ( niunit, fmt = * )( a( i, j ), j = 1, n )
827 READ( niunit, fmt = * )rcdein, rcdvin
829 CALL cget24( .true., 22, thresh, iseed, nounit, n, a, lda, h, ht,
830 $ w, wt, wtmp, vs, ldvs, vs1, rcdein, rcdvin, nslct
831 $ islct, isrt, result, work, lwork, rwork, bwork,
839 IF( result( j ).GE.zero )
841 IF( result( j ).GE.thresh )
846 $ ntestf = ntestf + 1
847 IF( ntestf.EQ.1 )
THEN
848 WRITE( nounit, fmt = 9999 )path
849 WRITE( nounit, fmt = 9998 )
850 WRITE( nounit, fmt = 9997 )
851 WRITE( nounit, fmt = 9996 )
852 WRITE( nounit, fmt = 9995 )thresh
853 WRITE( nounit, fmt = 9994 )
857 IF( result( j ).GE.thresh )
THEN
858 WRITE( nounit, fmt = 9992 )n, jtype, j, result( j )
862 nerrs = nerrs + nfail
863 ntestt = ntestt + ntest
869 CALL slasum( path, nounit, nerrs, ntestt )
871 9999
FORMAT( / 1x, a3,
' -- Complex Schur Form Decomposition Expert ',
872 $
'Driver', /
' Matrix types (see CDRVSX for details): ' )
874 9998
FORMAT( /
' Special Matrices:', / ' 1=zero matrix.
',
875 $ ' ', ' 5=diagonal: geometr. spaced entries.
',
876 $ / ' 2=identity matrix.
', ' 6=diagona
',
877 $ 'l: clustered entries.
', / ' 3=transposed jordan block.
',
878 $ ' ', ' 7=diagonal: large, evenly spaced.
', / ' ',
879 $ '4=diagonal: evenly spaced entries.
', ' 8=diagonal: s
',
880 $ 'mall, evenly spaced.
' )
881 9997 FORMAT( ' dense, non-symmetric matrices:
', / ' 9=well-cond., ev
',
882 $ 'enly spaced eigenvals.
', ' 14=ill-cond., geomet. spaced e
',
883 $ 'igenals.
', / ' 10=well-cond.,
geom. spaced eigenvals.
',
884 $ ' 15=ill-conditioned, clustered e.vals.
', / ' 11=well-cond
',
885 $ 'itioned, clustered e.vals.
', ' 16=ill-cond., random
comp',
886 $ 'lex
', / ' 12=well-cond., random
complex ', ' ',
887 $ ' 17=ill-cond., large rand. complx
', / ' 13=ill-condi',
888 $
'tioned, evenly spaced. ',
' 18=Ill-cond., small rand.',
890 9996
FORMAT(
' 19=Matrix with random O(1) entries. ',
' 21=Matrix ',
891 $
'with small random entries.', /
' 20=Matrix with large ran',
892 $
'dom entries. ', / )
893 9995
FORMAT(
' Tests performed with test threshold =', f8.2,
894 $ /
' ( A denotes A on input and T denotes A on output)',
895 $ / /
' 1 = 0 if T in Schur form (no sort), ',
896 $
' 1/ulp otherwise', /
897 $
' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
898 $ /
' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
899 $ /
' 4 = 0 if W are eigenvalues of T (no sort),',
900 $
' 1/ulp otherwise', /
901 $
' 5 = 0 if T same no matter if VS computed (no sort),',
902 $
' 1/ulp otherwise', /
903 $
' 6 = 0 if W same no matter if VS computed (no sort)',
904 $
', 1/ulp otherwise' )
905 9994
FORMAT(
' 7 = 0 if T in Schur form (sort), ',
' 1/ulp otherwise',
906 $ /
' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
907 $ /
' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
908 $ /
' 10 = 0 if W are eigenvalues of T (sort),',
909 $
' 1/ulp otherwise', /
910 $
' 11 = 0 if T same no matter what else computed (sort),',
911 $
' 1/ulp otherwise', /
912 $
' 12 = 0 if W same no matter what else computed ',
913 $
'(sort), 1/ulp otherwise', /
914 $
' 13 = 0 if sorting successful, 1/ulp otherwise',
915 $ /
' 14 = 0 if RCONDE same no matter what else computed,',
916 $
' 1/ulp otherwise', /
917 $
' 15 = 0 if RCONDv same no matter what else computed,',
918 $
' 1/ulp otherwise', /
919 $
' 16 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),',
920 $ /
' 17 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),' )
921 9993
FORMAT(
' N=', i5,
', IWK=', i2,
', seed=', 4( i4,
',' ),
922 $
' type ', i2,
', test(', i2,
')=', g10.3 )
923 9992
FORMAT(
' N=', i5,
', input example =', i3,
', test(', i2,
')=',
925 9991
FORMAT(
' CDRVSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
926 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )