498 SUBROUTINE zchkgg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
499 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
500 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1,
501 $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK,
502 $ RWORK, LLWORK, RESULT, INFO )
510 INTEGER , LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
511 DOUBLE PRECISION THRESH,
514 LOGICAL DOTYPE( * ), LLWORK(
516DOUBLE PRECISION RESULT( 15 ), RWORK( * )
517 COMPLEX*16 A( LDA, * ), ( * ), (
519 $ evectl( ldu, * ), evectr( ldu, * ),
520 $ h( lda, * ), p1( lda, * ), p2( lda, * ),
521 $ q( ldu, * ), s1( lda, * ), s2( lda, * ),
522 $ t( lda, * ), u( ldu, * ), v( ldu, * ),
523 $ work( * ), z( ldu, * )
529 DOUBLE PRECISION ZERO, ONE
530 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
531 COMPLEX*16 CZERO, CONE
532 PARAMETER ( = ( 0.0d+0, 0.0d+0 ),
533 $ cone = ( 1.0d+0, 0.0d+0 ) )
535 parameter( maxtyp = 26 )
539 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
540 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
542 DOUBLE PRECISION ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2,
547 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
548 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
549 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
552 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
553 DOUBLE PRECISION DUMMA( 4 ), RMAGN( 0: 3 )
554 COMPLEX*16 CDUMMA( 4 )
557 DOUBLE PRECISION DLAMCH, ZLANGE
559 EXTERNAL DLAMCH, ZLANGE, ZLARND
567 INTRINSIC abs, dble, dconjg,
max,
min, sign
570 DATA kclass / 15*1, 10*2, 1*3 /
571 DATA kz1 / 0, 1, 2, 1, 3, 3 /
572 DATA kz2 / 0, 0, 1, 2, 1, 1 /
573 DATA kadd / 0, 0, 0, 0, 3, 2 /
574 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
575 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
577 $ 1, 1, -4, 2, -4, 8*8, 0 /
578 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
580 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
582 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
584 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
586 DATA ktrian / 16*0, 10*1 /
587 DATA lasign / 6*.false., .true., .false., 2*.true.,
588 $ 2*.false., 3*.true., .false., .true.,
589 $ 3*.false., 5*.true., .false. /
590 DATA lbsign / 7*.false., .true., 2*.false.,
591 $ 2*.true., 2*.false., .true., .false., .true.,
603 nmax =
max( nmax, nn( j ) )
608 lwkopt =
max( 2*nmax*nmax, 4*nmax, 1 )
612 IF( nsizes.LT.0 )
THEN
614 ELSE IF( badnn )
THEN
616 ELSE IF( ntypes.LT.0 )
THEN
618 ELSE IF( thresh.LT.zero )
THEN
620 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
622 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax )
THEN
624 ELSE IF( lwkopt.GT.lwork )
THEN
629 CALL xerbla(
'ZCHKGG', -info )
635 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
638 safmin = dlamch(
'Safe minimum' )
639 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
640 safmin = safmin / ulp
641 safmax = one / safmin
642 CALL dlabad( safmin, safmax )
656 DO 240 jsize = 1, nsizes
659 rmagn( 2 ) = safmax*ulp / dble( n1 )
660 rmagn( 3 ) = safmin*ulpinv*n1
662 IF( nsizes.NE.1 )
THEN
663 mtypes =
min( maxtyp, ntypes )
665 mtypes =
min( maxtyp+1, ntypes )
668 DO 230 jtype = 1, mtypes
669 IF( .NOT.dotype( jtype ) )
677 ioldsd( j ) = iseed( j )
707 IF( mtypes.GT.maxtyp )
710 IF( kclass( jtype ).LT.3 )
THEN
714 IF( abs( katype( jtype ) ).EQ.3 )
THEN
715 in = 2*( ( n-1 ) / 2 ) + 1
717 $
CALL zlaset( 'full
', N, N, CZERO, CZERO, A, LDA )
721 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
722 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
723 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
724 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4,
726 IADD = KADD( KAZERO( JTYPE ) )
727.GT..AND..LE.
IF( IADD0 IADDN )
728 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
732.EQ.
IF( ABS( KBTYPE( JTYPE ) )3 ) THEN
733 IN = 2*( ( N-1 ) / 2 ) + 1
735 $ CALL ZLASET( 'full
', N, N, CZERO, CZERO, B, LDA )
739 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
740 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
741 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
742 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4,
744 IADD = KADD( KBZERO( JTYPE ) )
746 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
748.EQ..AND..GT.
IF( KCLASS( JTYPE )2 N0 ) THEN
758 U( JR, JC ) = ZLARND( 3, ISEED )
759 V( JR, JC ) = ZLARND( 3, ISEED )
761 CALL ZLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1,
763 WORK( 2*N+JC ) = SIGN( ONE, DBLE( U( JC, JC ) ) )
765 CALL ZLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1,
767 WORK( 3*N+JC ) = SIGN( ONE, DBLE( V( JC, JC ) ) )
770 CTEMP = ZLARND( 3, ISEED )
773 WORK( 3*N ) = CTEMP / ABS( CTEMP )
774 CTEMP = ZLARND( 3, ISEED )
777 WORK( 4*N ) = CTEMP / ABS( CTEMP )
783 A( JR, JC ) = WORK( 2*N+JR )*
784 $ DCONJG( WORK( 3*N+JC ) )*
786 B( JR, JC ) = WORK( 2*N+JR )*
787 $ DCONJG( WORK( 3*N+JC ) )*
791 CALL ZUNM2R( 'l
', 'n
', N, N, N-1, U, LDU, WORK, A,
792 $ LDA, WORK( 2*N+1 ), IINFO )
795 CALL ZUNM2R( 'r
', 'c
', N, N, N-1, V, LDU, WORK( N+1 ),
796 $ A, LDA, WORK( 2*N+1 ), IINFO )
799 CALL ZUNM2R( 'l
', 'n
', N, N, N-1, U, LDU, WORK, B,
800 $ LDA, WORK( 2*N+1 ), IINFO )
803 CALL ZUNM2R( 'r
', 'c
', N, N, N-1, V, LDU, WORK( N+1 ),
804 $ B, LDA, WORK( 2*N+1 ), IINFO )
814 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
816 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
822 ANORM = ZLANGE( '1
', N, N, A, LDA, RWORK )
823 BNORM = ZLANGE( '1
', N, N, B, LDA, RWORK )
827.NE.
IF( IINFO0 ) THEN
828 WRITE( NOUNIT, FMT = 9999 )'generator
', IINFO, N, JTYPE,
838 CALL ZLACPY( ' ', N, N, A, LDA, H, LDA )
839 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
843 CALL ZGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO )
844.NE.
IF( IINFO0 ) THEN
845 WRITE( NOUNIT, FMT = 9999 )'zgeqr2', IINFO, N, JTYPE,
851 CALL ZUNM2R( 'l
', 'c
', N, N, N, T, LDA, WORK, H, LDA,
852 $ WORK( N+1 ), IINFO )
853.NE.
IF( IINFO0 ) THEN
854 WRITE( NOUNIT, FMT = 9999 )'zunm2r', IINFO, N, JTYPE,
860 CALL ZLASET( 'full
', N, N, CZERO, CONE, U, LDU )
861 CALL ZUNM2R( 'r
', 'n
', N, N, N, T, LDA, WORK, U, LDU,
862 $ WORK( N+1 ), IINFO )
863.NE.
IF( IINFO0 ) THEN
864 WRITE( NOUNIT, FMT = 9999 )'zunm2r', IINFO, N, JTYPE,
870 CALL ZGGHRD( 'v
', 'i
', N, 1, N, H, LDA, T, LDA, U, LDU, V,
872.NE.
IF( IINFO0 ) THEN
873 WRITE( NOUNIT, FMT = 9999 )'zgghrd', IINFO, N, JTYPE,
882 CALL ZGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK,
883 $ RWORK, RESULT( 1 ) )
884 CALL ZGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK,
885 $ RWORK, RESULT( 2 ) )
886 CALL ZGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK,
887 $ RWORK, RESULT( 3 ) )
888 CALL ZGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK,
889 $ RWORK, RESULT( 4 ) )
897 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
898 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
902 CALL ZHGEQZ( 'e
', 'n
', 'n
', N, 1, N, S2, LDA, P2, LDA,
903 $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK,
905.NE.
IF( IINFO0 ) THEN
906 WRITE( NOUNIT, FMT = 9999 )'zhgeqz(e)
', IINFO, N, JTYPE,
914 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
915 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
917 CALL ZHGEQZ( 's
', 'n
', 'n
', N, 1, N, S2, LDA, P2, LDA,
918 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
920.NE.
IF( IINFO0 ) THEN
921 WRITE( NOUNIT, FMT = 9999 )'zhgeqz(s)
', IINFO, N, JTYPE,
929 CALL ZLACPY( ' ', N, N, H, LDA, S1, LDA )
930 CALL ZLACPY( ' ', N, N, T, LDA, P1, LDA )
932 CALL ZHGEQZ( 's
', 'i
', 'i
', N, 1, N, S1, LDA, P1, LDA,
933 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
935.NE.
IF( IINFO0 ) THEN
936 WRITE( NOUNIT, FMT = 9999 )'zhgeqz(v)
', IINFO, N, JTYPE,
946 CALL ZGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
947 $ RWORK, RESULT( 5 ) )
948 CALL ZGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
949 $ RWORK, RESULT( 6 ) )
950 CALL ZGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
951 $ RWORK, RESULT( 7 ) )
952 CALL ZGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
953 $ RWORK, RESULT( 8 ) )
971 LLWORK( J ) = .FALSE.
974 CALL ZTGEVC( 'l
', 's
', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
975 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
976.NE.
IF( IINFO0 ) THEN
977 WRITE( NOUNIT, FMT = 9999 )'ztgevc(l,s1)
', IINFO, N,
985 LLWORK( J ) = .FALSE.
991 CALL ZTGEVC( 'l',
'S', llwork, n, s1, lda, p1, lda,
992 $ evectl( 1, i1+1 ), ldu, cdumma, ldu, n, in,
993 $ work, rwork, iinfo )
994 IF( iinfo.NE.0 )
THEN
995 WRITE( nounit, fmt = 9999 )
'ZTGEVC(L,S2)', iinfo, n,
1001 CALL zget52( .true., n, s1, lda, p1, lda, evectl, ldu,
1002 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1003 result( 9 ) = dumma( 1 )
1004 IF( dumma( 2 ).GT.thrshn )
THEN
1005 WRITE( nounit, fmt = 9998 )
'Left',
'ZTGEVC(HOWMNY=S)',
1006 $ dumma( 2 ), n, jtype, ioldsd
1013 result( 10 ) = ulpinv
1014 CALL zlacpy(
'F', n, n, q, ldu, evectl, ldu )
1015 CALL ztgevc(
'L',
'B', llwork, n, s1, lda, p1, lda, evectl,
1016 $ ldu, cdumma, ldu, n, in, work, rwork, iinfo )
1017 IF( iinfo.NE.0 )
THEN
1018 WRITE( nounit, fmt = 9999 )
'ZTGEVC(L,B)', iinfo, n,
1024 CALL zget52( .true., n, h, lda, t, lda, evectl, ldu, alpha1,
1025 $ beta1, work, rwork, dumma( 1 ) )
1026 result( 10 ) = dumma( 1 )
1027 IF( dumma( 2 ).GT.thrshn )
THEN
1028 WRITE( nounit, fmt = 9998 )
'Left',
'ZTGEVC(HOWMNY=B)',
1029 $ dumma( 2 ), n, jtype, ioldsd
1036 result( 11 ) = ulpinv
1043 llwork( j ) = .true.
1045 DO 170 j = i1 + 1, n
1046 llwork( j ) = .false.
1049 CALL ztgevc( 'r
', 's
', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
1050 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
1051.NE.
IF( IINFO0 ) THEN
1052 WRITE( NOUNIT, FMT = 9999 )'ztgevc(r,s1)
', IINFO, N,
1060 LLWORK( J ) = .FALSE.
1062 DO 190 J = I1 + 1, N
1063 LLWORK( J ) = .TRUE.
1066 CALL ZTGEVC( 'r
', 's
', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
1067 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
1069.NE.
IF( IINFO0 ) THEN
1070 WRITE( NOUNIT, FMT = 9999 )'ztgevc(r,s2)
', IINFO, N,
1076 CALL ZGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU,
1077 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
1078 RESULT( 11 ) = DUMMA( 1 )
1079.GT.
IF( DUMMA( 2 )THRESH ) THEN
1080 WRITE( NOUNIT, FMT = 9998 )'right
', 'ztgevc(howmny=s)
',
1081 $ DUMMA( 2 ), N, JTYPE, IOLDSD
1088 RESULT( 12 ) = ULPINV
1089 CALL ZLACPY( 'f
', N, N, Z, LDU, EVECTR, LDU )
1090 CALL ZTGEVC( 'r
', 'b', llwork, n, s1, lda, p1, lda, cdumma,
1091 $ ldu, evectr, ldu, n, in, work, rwork, iinfo )
1092 IF( iinfo.NE.0 )
THEN
1093 WRITE( nounit, fmt = 9999 )
'ZTGEVC(R,B)', iinfo, n,
1100 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1101 result( 12 ) = dumma( 1 )
1102 IF( dumma( 2 ).GT.thresh )
THEN
1103 WRITE( nounit, fmt = 9998 )
'Right',
'ZTGEVC(HOWMNY=B)',
1113 CALL zget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1114 $ work, rwork, result( 13 ) )
1115 CALL zget51( 2, n, p1, lda, p2, lda, q, ldu, z, ldu,
1116 $ work, rwork, result( 14 ) )
1123 temp1 =
max( temp1, abs( alpha1( j )-alpha3( j ) ) )
1124 temp2 =
max( temp2, abs( beta1( j )-beta3( j ) ) )
1127 temp1 = temp1 /
max( safmin, ulp*
max( temp1, anorm ) )
1128 temp2 = temp2 /
max( safmin, ulp*
max( temp2, bnorm ) )
1129 result( 15 ) =
max( temp1, temp2 )
1142 ntestt = ntestt + ntest
1146 DO 220 jr = 1, ntest
1147 IF( result( jr ).GE.thresh )
THEN
1152 IF( nerrs.EQ.0 )
THEN
1153 WRITE( nounit, fmt = 9997 )
'ZGG'
1157 WRITE( nounit, fmt = 9996 )
1158 WRITE( nounit, fmt = 9995 )
1159 WRITE( nounit, fmt = 9994 )
'Unitary'
1163 WRITE( nounit, fmt = 9993 )
'unitary',
'*',
1164 $
'conjugate transpose', (
'*', j = 1, 10 )
1168 IF( result( jr ).LT.10000.0d0 )
THEN
1169 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1172 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1183 CALL dlasum(
'ZGG', nounit, nerrs, ntestt )
1186 9999
FORMAT(
' ZCHKGG: ', a
' returned INFO=''.''N=',
1187 $ i6,
', JTYPE=', i6,
', ISEED=('',' ), i5,
')' )
1189 9998
FORMAT(
' ZCHKGG: ', a,
' Eigenvectors from ', a,
' incorrectly ',
1190 $
'normalized.', /
' Bits of error=', 0p, g10.3,
','
1191 $
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
','
1194 9997
FORMAT( 1x, a3,
' -- Complex Generalized eigenvalue problem'
1196 9996
FORMAT(
' Matrix types (see ZCHKGG for details): ' )
1198 9995
FORMAT(
' Special Matrices:', 23x,
1199 $
'(J''=transposed Jordan block)',
1200 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
1201 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
1202 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
1203 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
1204 $ ' 8=(i,d) 10=(small*d, large*i) 12=(small*i, large*d) ',
1205 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
1206 9994
FORMAT(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
1207 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
1208 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
1209 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
1210 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
1211 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
1212 $
'23=(small,large) 24=(small,small) 25=(large,large)',
1213 $ /
' 26=random O(1) matrices.' )
1215 9993
FORMAT( /
' Tests performed: (H is Hessenberg, S is Schur, B, ',
1216 $
'T, P are triangular,', / 20x,
'U, V, Q, and Z are ', a,
1217 $
', l and r are the', / 20x,
1218 $
'appropriate left and right eigenvectors, resp., a is',
1219 $ / 20x,
'alpha, b is beta, and ', a,
' means ', a,
'.)',
1220 $ /
' 1 = | A - U H V', a,
1221 $
' | / ( |A| n ulp ) 2 = | B - U T V', a,
1222 $
' | / ( |B| n ulp )', /
' 3 = | I - UU', a,
1223 $
' | / ( n ulp ) 4 = | I - VV', a,
1224 $
' | / ( n ulp )', /
' 5 = | H - Q S Z', a,
1225 $
' | / ( |H| n ulp )', 6x,
'6 = | T - Q P Z', a,
1226 $
' | / ( |T| n ulp )', /
' 7 = | I - QQ', a,
1227 $
' | / ( n ulp ) 8 = | I - ZZ', a,
1228 $
' | / ( n ulp )', /
' 9 = max | ( b S - a P )', a,
1229 $
' l | / const. 10 = max | ( b H - a T )', a,
1230 $
' l | / const.', /
1231 $
' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
1232 $
' - a T ) r | / const.', / 1x )
1234 9992
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
1235 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
1236 9991
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
1237 $ 4( i4, ',
' ), ' result
', I2, ' is
', 1P, D10.3 )