377 SUBROUTINE cdrves( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
378 $ nounit, a, lda, h, ht, w, wt, vs, ldvs, result,
379 $ work, nwork, rwork, iwork, bwork, info )
387 INTEGER info, lda, ldvs, nounit, nsizes, ntypes, nwork
391 LOGICAL bwork( * ), dotype( * )
392 INTEGER iseed( 4 ), iwork( * ), nn( * )
393 REAL result( 13 ), rwork( * )
394 COMPLEX a( lda, * ), h( lda, * ), ht( lda, * ),
395 $ vs( ldvs, * ), w( * ), work( * ), wt( * )
402 parameter( czero = ( 0.0e+0, 0.0e+0 ) )
404 parameter( cone = ( 1.0e+0, 0.0e+0 ) )
406 parameter( zero = 0.0e+0, one = 1.0e+0 )
408 parameter( maxtyp = 21 )
414 INTEGER i, iinfo, imode, isort, itype, iwk, j, jcol,
415 $ jsize, jtype, knteig, lwork, mtypes, n,
416 $ nerrs, nfail, nmax, nnwork, ntest, ntestf,
418 REAL anorm, cond, conds, ovfl, rtulp, rtulpi, ulp,
422 INTEGER idumma( 1 ), ioldsd( 4 ), kconds( maxtyp ),
423 $ kmagn( maxtyp ), kmode( maxtyp ),
429 REAL selwi( 20 ), selwr( 20 )
432 INTEGER seldim, selopt
435 common / sslct / selopt, seldim, selval, selwr, selwi
447 INTRINSIC abs, cmplx, max, min, sqrt
450 DATA ktype / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
451 DATA kmagn / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
453 DATA kmode / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
454 $ 1, 5, 5, 5, 4, 3, 1 /
455 DATA kconds / 3*0, 5*0, 4*1, 6*2, 3*0 /
459 path( 1: 1 ) =
'Complex precision'
474 nmax = max( nmax, nn( j ) )
481 IF( nsizes.LT.0 )
THEN
483 ELSE IF( badnn )
THEN
485 ELSE IF( ntypes.LT.0 )
THEN
487 ELSE IF( thresh.LT.zero )
THEN
489 ELSE IF( nounit.LE.0 )
THEN
491 ELSE IF( lda.LT.1 .OR. lda.LT.nmax )
THEN
493 ELSE IF( ldvs.LT.1 .OR. ldvs.LT.nmax )
THEN
495 ELSE IF( 5*nmax+2*nmax**2.GT.nwork )
THEN
500 CALL
xerbla(
'CDRVES', -info )
506 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
511 unfl =
slamch(
'Safe minimum' )
514 ulp =
slamch(
'Precision' )
523 DO 240 jsize = 1, nsizes
525 IF( nsizes.NE.1 )
THEN
526 mtypes = min( maxtyp, ntypes )
528 mtypes = min( maxtyp+1, ntypes )
531 DO 230 jtype = 1, mtypes
532 IF( .NOT.dotype( jtype ) )
538 ioldsd( j ) = iseed( j )
557 IF( mtypes.GT.maxtyp )
560 itype = ktype( jtype )
561 imode = kmode( jtype )
565 go to( 30, 40, 50 )kmagn( jtype )
581 CALL
claset(
'Full', lda, n, czero, czero, a, lda )
587 IF( itype.EQ.1 )
THEN
593 ELSE IF( itype.EQ.2 )
THEN
598 a( jcol, jcol ) = cmplx( anorm )
601 ELSE IF( itype.EQ.3 )
THEN
606 a( jcol, jcol ) = cmplx( anorm )
608 $ a( jcol, jcol-1 ) = cone
611 ELSE IF( itype.EQ.4 )
THEN
615 CALL
clatms( n, n,
'S', iseed,
'H', rwork, imode, cond,
616 $ anorm, 0, 0,
'N', a, lda, work( n+1 ),
619 ELSE IF( itype.EQ.5 )
THEN
623 CALL
clatms( n, n,
'S', iseed,
'H', rwork, 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
639 CALL
clatme( n,
'D', iseed, work, imode, cond, cone,
640 $
'T',
'T',
'T', rwork, 4, conds, n, n, anorm,
641 $ a, lda, work( 2*n+1 ), iinfo )
643 ELSE IF( itype.EQ.7 )
THEN
647 CALL
clatmr( n, n,
'D', iseed,
'N', work, 6, one, cone,
648 $
'T',
'N', work( n+1 ), 1, one,
649 $ work( 2*n+1 ), 1, one,
'N', idumma, 0, 0,
650 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
652 ELSE IF( itype.EQ.8 )
THEN
656 CALL
clatmr( n, n,
'D', iseed,
'H', work, 6, one, cone,
657 $
'T',
'N', work( n+1 ), 1, one,
658 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
659 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
661 ELSE IF( itype.EQ.9 )
THEN
665 CALL
clatmr( n, n,
'D', iseed,
'N', work, 6, one, cone,
666 $
'T',
'N', work( n+1 ), 1, one,
667 $ work( 2*n+1 ), 1, one,
'N', idumma, n, n,
668 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
670 CALL
claset(
'Full', 2, n, czero, czero, a, lda )
671 CALL
claset(
'Full', n-3, 1, czero, czero, a( 3, 1 ),
673 CALL
claset(
'Full', n-3, 2, czero, czero,
675 CALL
claset(
'Full', 1, n, czero, czero, a( n, 1 ),
679 ELSE IF( itype.EQ.10 )
THEN
683 CALL
clatmr( n, n,
'D', iseed,
'N', work, 6, one, cone,
684 $
'T',
'N', work( n+1 ), 1, one,
685 $ work( 2*n+1 ), 1, one,
'N', idumma, n, 0,
686 $ zero, anorm,
'NO', a, lda, iwork, iinfo )
693 IF( iinfo.NE.0 )
THEN
694 WRITE( nounit, fmt = 9992 )
'Generator', iinfo, n, jtype,
708 nnwork = 5*n + 2*n**2
710 nnwork = max( nnwork, 1 )
721 IF( isort.EQ.0 )
THEN
731 CALL
clacpy(
'F', n, n, a, lda, h, lda )
732 CALL
cgees(
'V', sort,
cslect, n, h, lda, sdim, w, vs,
733 $ ldvs, work, nnwork, rwork, bwork, iinfo )
734 IF( iinfo.NE.0 )
THEN
735 result( 1+rsub ) = ulpinv
736 WRITE( nounit, fmt = 9992 )
'CGEES1', iinfo, n,
744 result( 1+rsub ) = zero
747 IF( h( i, j ).NE.zero )
748 $ result( 1+rsub ) = ulpinv
754 lwork = max( 1, 2*n*n )
755 CALL
chst01( n, 1, n, a, lda, h, lda, vs, ldvs, work,
756 $ lwork, rwork, res )
757 result( 2+rsub ) = res( 1 )
758 result( 3+rsub ) = res( 2 )
762 result( 4+rsub ) = zero
764 IF( h( i, i ).NE.w( i ) )
765 $ result( 4+rsub ) = ulpinv
770 CALL
clacpy(
'F', n, n, a, lda, ht, lda )
771 CALL
cgees(
'N', sort,
cslect, n, ht, lda, sdim, wt,
772 $ vs, ldvs, work, nnwork, rwork, bwork,
774 IF( iinfo.NE.0 )
THEN
775 result( 5+rsub ) = ulpinv
776 WRITE( nounit, fmt = 9992 )
'CGEES2', iinfo, n,
782 result( 5+rsub ) = zero
785 IF( h( i, j ).NE.ht( i, j ) )
786 $ result( 5+rsub ) = ulpinv
792 result( 6+rsub ) = zero
794 IF( w( i ).NE.wt( i ) )
795 $ result( 6+rsub ) = ulpinv
800 IF( isort.EQ.1 )
THEN
805 $ knteig = knteig + 1
807 IF(
cslect( w( i+1 ) ) .AND.
808 $ ( .NOT.
cslect( w( i ) ) ) )result( 13 )
813 $ result( 13 ) = ulpinv
825 IF( result( j ).GE.zero )
827 IF( result( j ).GE.thresh )
832 $ ntestf = ntestf + 1
833 IF( ntestf.EQ.1 )
THEN
834 WRITE( nounit, fmt = 9999 )path
835 WRITE( nounit, fmt = 9998 )
836 WRITE( nounit, fmt = 9997 )
837 WRITE( nounit, fmt = 9996 )
838 WRITE( nounit, fmt = 9995 )thresh
839 WRITE( nounit, fmt = 9994 )
844 IF( result( j ).GE.thresh )
THEN
845 WRITE( nounit, fmt = 9993 )n, iwk, ioldsd, jtype,
850 nerrs = nerrs + nfail
851 ntestt = ntestt + ntest
859 CALL
slasum( path, nounit, nerrs, ntestt )
861 9999 format( / 1x, a3,
' -- Complex Schur Form Decomposition Driver',
862 $ /
' Matrix types (see CDRVES for details): ' )
864 9998 format( /
' Special Matrices:', /
' 1=Zero matrix. ',
865 $
' ',
' 5=Diagonal: geometr. spaced entries.',
866 $ /
' 2=Identity matrix. ',
' 6=Diagona',
867 $
'l: clustered entries.', /
' 3=Transposed Jordan block. ',
868 $
' ',
' 7=Diagonal: large, evenly spaced.', /
' ',
869 $
'4=Diagonal: evenly spaced entries. ',
' 8=Diagonal: s',
870 $
'mall, evenly spaced.' )
871 9997 format(
' Dense, Non-Symmetric Matrices:', /
' 9=Well-cond., ev',
872 $
'enly spaced eigenvals.',
' 14=Ill-cond., geomet. spaced e',
873 $
'igenals.', /
' 10=Well-cond., geom. spaced eigenvals. ',
874 $
' 15=Ill-conditioned, clustered e.vals.', /
' 11=Well-cond',
875 $
'itioned, clustered e.vals. ',
' 16=Ill-cond., random comp',
876 $
'lex ', a6, /
' 12=Well-cond., random complex ', a6,
' ',
877 $
' 17=Ill-cond., large rand. complx ', a4, /
' 13=Ill-condi',
878 $
'tioned, evenly spaced. ',
' 18=Ill-cond., small rand.',
880 9996 format(
' 19=Matrix with random O(1) entries. ',
' 21=Matrix ',
881 $
'with small random entries.', /
' 20=Matrix with large ran',
882 $
'dom entries. ', / )
883 9995 format(
' Tests performed with test threshold =', f8.2,
884 $ /
' ( A denotes A on input and T denotes A on output)',
885 $ / /
' 1 = 0 if T in Schur form (no sort), ',
886 $
' 1/ulp otherwise', /
887 $
' 2 = | A - VS T transpose(VS) | / ( n |A| ulp ) (no sort)',
888 $ /
' 3 = | I - VS transpose(VS) | / ( n ulp ) (no sort) ',
889 $ /
' 4 = 0 if W are eigenvalues of T (no sort),',
890 $
' 1/ulp otherwise', /
891 $
' 5 = 0 if T same no matter if VS computed (no sort),',
892 $
' 1/ulp otherwise', /
893 $
' 6 = 0 if W same no matter if VS computed (no sort)',
894 $
', 1/ulp otherwise' )
895 9994 format(
' 7 = 0 if T in Schur form (sort), ',
' 1/ulp otherwise',
896 $ /
' 8 = | A - VS T transpose(VS) | / ( n |A| ulp ) (sort)',
897 $ /
' 9 = | I - VS transpose(VS) | / ( n ulp ) (sort) ',
898 $ /
' 10 = 0 if W are eigenvalues of T (sort),',
899 $
' 1/ulp otherwise', /
900 $
' 11 = 0 if T same no matter if VS computed (sort),',
901 $
' 1/ulp otherwise', /
902 $
' 12 = 0 if W same no matter if VS computed (sort),',
903 $
' 1/ulp otherwise', /
904 $
' 13 = 0 if sorting succesful, 1/ulp otherwise', / )
905 9993 format(
' N=', i5,
', IWK=', i2,
', seed=', 4( i4,
',' ),
906 $
' type ', i2,
', test(', i2,
')=', g10.3 )
907 9992 format(
' CDRVES: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
908 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )