498 SUBROUTINE cchkgg( 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 INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
514 LOGICAL DOTYPE( * ), LLWORK( * )
515 INTEGER ISEED( 4 ), NN( * )
516 REAL RESULT( 15 ), RWORK( * )
517 COMPLEX A( LDA, * ), ALPHA1( * ), ALPHA3( * ),
518 $ b( lda, * ), beta1( * ), beta3( * ),
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, * )
530 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
532 PARAMETER ( CZERO = ( 0.0e+0, 0.0e+0 ),
533 $ cone = ( 1.0e+0, 0.0e+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 REAL 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 ),
550 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
551 $ kbzero( maxtyp ), kclass( maxtyp ),
552 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
553 REAL DUMMA( 4 ), RMAGN( 0: 3 )
559 EXTERNAL CLANGE, SLAMCH, CLARND
567 INTRINSIC abs, conjg, max, min, real, 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 /
576 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
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(
'CCHKGG', -info )
635 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
638 safmin = slamch(
'Safe minimum' )
639 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
640 safmin = safmin / ulp
641 safmax = one / safmin
655 DO 240 jsize = 1, nsizes
658 rmagn( 2 ) = safmax*ulp / real( n1 )
659 rmagn( 3 ) = safmin*ulpinv*n1
661 IF( nsizes.NE.1 )
THEN
662 mtypes = min( maxtyp, ntypes )
664 mtypes = min( maxtyp+1, ntypes )
667 DO 230 jtype = 1, mtypes
668 IF( .NOT.dotype( jtype ) )
676 ioldsd( j ) = iseed( j )
706 IF( mtypes.GT.maxtyp )
709 IF( kclass( jtype ).LT.3 )
THEN
713 IF( abs( katype( jtype ) ).EQ.3 )
THEN
714 in = 2*( ( n-1 ) / 2 ) + 1
716 $
CALL claset(
'Full', n, n, czero, czero, a, lda )
720 CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
721 $ kz2( kazero( jtype ) ), lasign( jtype ),
722 $ rmagn( kamagn( jtype ) ), ulp,
723 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 4,
725 iadd = kadd( kazero( jtype ) )
726 IF( iadd.GT.0 .AND. iadd.LE.n )
727 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
731 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
732 in = 2*( ( n-1 ) / 2 ) + 1
734 $
CALL claset(
'Full', n, n, czero, czero, b, lda )
738 CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
739 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
740 $ rmagn( kbmagn( jtype ) ), one,
741 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 4,
743 iadd = kadd( kbzero( jtype ) )
745 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
747 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
757 u( jr, jc ) = clarnd( 3, iseed )
758 v( jr, jc ) = clarnd( 3, iseed )
760 CALL clarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
762 work( 2*n+jc ) = sign( one, real( u( jc, jc ) ) )
764 CALL clarfg( n+1-jc, v( jc, jc ), v( jc+1, jc ), 1,
766 work( 3*n+jc ) = sign( one, real( v( jc, jc ) ) )
769 ctemp = clarnd( 3, iseed )
772 work( 3*n ) = ctemp / abs( ctemp )
773 ctemp = clarnd( 3, iseed )
776 work( 4*n ) = ctemp / abs( ctemp )
782 a( jr, jc ) = work( 2*n+jr )*
783 $ conjg( work( 3*n+jc ) )*
785 b( jr, jc ) = work( 2*n+jr )*
786 $ conjg( work( 3*n+jc ) )*
790 CALL cunm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
791 $ lda, work( 2*n+1 ), iinfo )
794 CALL cunm2r(
'R',
'C', n, n, n-1, v, ldu, work( n+1 ),
795 $ a, lda, work( 2*n+1 ), iinfo )
798 CALL cunm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
799 $ lda, work( 2*n+1 ), iinfo )
802 CALL cunm2r(
'R',
'C', n, n, n-1, v, ldu, work( n+1 ),
803 $ b, lda, work( 2*n+1 ), iinfo )
813 a( jr, jc ) = rmagn( kamagn( jtype ) )*
815 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
821 anorm = clange(
'1', n, n, a, lda, rwork )
822 bnorm = clange(
'1', n, n, b, lda, rwork )
826 IF( iinfo.NE.0 )
THEN
827 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
837 CALL clacpy(
' ', n, n, a, lda, h, lda )
838 CALL clacpy(
' ', n, n, b, lda, t, lda )
842 CALL cgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
843 IF( iinfo.NE.0 )
THEN
844 WRITE( nounit, fmt = 9999 )
'CGEQR2', iinfo, n, jtype,
850 CALL cunm2r(
'L',
'C', n, n, n, t, lda, work, h, lda,
851 $ work( n+1 ), iinfo )
852 IF( iinfo.NE.0 )
THEN
853 WRITE( nounit, fmt = 9999 )
'CUNM2R', iinfo, n, jtype,
859 CALL claset(
'Full', n, n, czero, cone, u, ldu )
860 CALL cunm2r(
'R',
'N', n, n, n, t, lda, work, u, ldu,
861 $ work( n+1 ), iinfo )
862 IF( iinfo.NE.0 )
THEN
863 WRITE( nounit, fmt = 9999 )
'CUNM2R', iinfo, n, jtype,
869 CALL cgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
871 IF( iinfo.NE.0 )
THEN
872 WRITE( nounit, fmt = 9999 )
'CGGHRD', iinfo, n, jtype,
881 CALL cget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
882 $ rwork, result( 1 ) )
883 CALL cget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
884 $ rwork, result( 2 ) )
885 CALL cget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
886 $ rwork, result( 3 ) )
887 CALL cget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
888 $ rwork, result( 4 ) )
896 CALL clacpy(
' ', n, n, h, lda, s2, lda )
897 CALL clacpy(
' ', n, n, t, lda, p2, lda )
901 CALL chgeqz(
'E',
'N',
'N', n, 1, n, s2, lda, p2, lda,
902 $ alpha3, beta3, q, ldu, z, ldu, work, lwork,
904 IF( iinfo.NE.0 )
THEN
905 WRITE( nounit, fmt = 9999 )
'CHGEQZ(E)', iinfo, n, jtype,
913 CALL clacpy(
' ', n, n, h, lda, s2, lda )
914 CALL clacpy(
' ', n, n, t, lda, p2, lda )
916 CALL chgeqz(
'S',
'N',
'N', n, 1, n, s2, lda, p2, lda,
917 $ alpha1, beta1, q, ldu, z, ldu, work, lwork,
919 IF( iinfo.NE.0 )
THEN
920 WRITE( nounit, fmt = 9999 )
'CHGEQZ(S)', iinfo, n, jtype,
928 CALL clacpy(
' ', n, n, h, lda, s1, lda )
929 CALL clacpy(
' ', n, n, t, lda, p1, lda )
931 CALL chgeqz(
'S',
'I',
'I', n, 1, n, s1, lda, p1, lda,
932 $ alpha1, beta1, q, ldu, z, ldu, work, lwork,
934 IF( iinfo.NE.0 )
THEN
935 WRITE( nounit, fmt = 9999 )
'CHGEQZ(V)', iinfo, n, jtype,
945 CALL cget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
946 $ rwork, result( 5 ) )
947 CALL cget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
948 $ rwork, result( 6 ) )
949 CALL cget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
950 $ rwork, result( 7 ) )
951 CALL cget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
952 $ rwork, result( 8 ) )
970 llwork( j ) = .false.
973 CALL ctgevc(
'L',
'S', llwork, n, s1, lda, p1, lda, evectl,
974 $ ldu, cdumma, ldu, n, in, work, rwork, iinfo )
975 IF( iinfo.NE.0 )
THEN
976 WRITE( nounit, fmt = 9999 )
'CTGEVC(L,S1)', iinfo, n,
984 llwork( j ) = .false.
990 CALL ctgevc(
'L',
'S', llwork, n, s1, lda, p1, lda,
991 $ evectl( 1, i1+1 ), ldu, cdumma, ldu, n, in,
992 $ work, rwork, iinfo )
993 IF( iinfo.NE.0 )
THEN
994 WRITE( nounit, fmt = 9999 )
'CTGEVC(L,S2)', iinfo, n,
1000 CALL cget52( .true., n, s1, lda, p1, lda, evectl, ldu,
1001 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1002 result( 9 ) = dumma( 1 )
1003 IF( dumma( 2 ).GT.thrshn )
THEN
1004 WRITE( nounit, fmt = 9998 )
'Left',
'CTGEVC(HOWMNY=S)',
1005 $ dumma( 2 ), n, jtype, ioldsd
1012 result( 10 ) = ulpinv
1013 CALL clacpy(
'F', n, n, q, ldu, evectl, ldu )
1014 CALL ctgevc(
'L',
'B', llwork, n, s1, lda, p1, lda, evectl,
1015 $ ldu, cdumma, ldu, n, in, work, rwork, iinfo )
1016 IF( iinfo.NE.0 )
THEN
1017 WRITE( nounit, fmt = 9999 )
'CTGEVC(L,B)', iinfo, n,
1023 CALL cget52( .true., n, h, lda, t, lda, evectl, ldu, alpha1,
1024 $ beta1, work, rwork, dumma( 1 ) )
1025 result( 10 ) = dumma( 1 )
1026 IF( dumma( 2 ).GT.thrshn )
THEN
1027 WRITE( nounit, fmt = 9998 )
'Left',
'CTGEVC(HOWMNY=B)',
1028 $ dumma( 2 ), n, jtype, ioldsd
1035 result( 11 ) = ulpinv
1042 llwork( j ) = .true.
1044 DO 170 j = i1 + 1, n
1045 llwork( j ) = .false.
1048 CALL ctgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, cdumma,
1049 $ ldu, evectr, ldu, n, in, work, rwork, iinfo )
1050 IF( iinfo.NE.0 )
THEN
1051 WRITE( nounit, fmt = 9999 )
'CTGEVC(R,S1)', iinfo, n,
1059 llwork( j ) = .false.
1061 DO 190 j = i1 + 1, n
1062 llwork( j ) = .true.
1065 CALL ctgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, cdumma,
1066 $ ldu, evectr( 1, i1+1 ), ldu, n, in, work,
1068 IF( iinfo.NE.0 )
THEN
1069 WRITE( nounit, fmt = 9999 )
'CTGEVC(R,S2)', iinfo, n,
1075 CALL cget52( .false., n, s1, lda, p1, lda, evectr, ldu,
1076 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1077 result( 11 ) = dumma( 1 )
1078 IF( dumma( 2 ).GT.thresh )
THEN
1079 WRITE( nounit, fmt = 9998 )
'Right',
'CTGEVC(HOWMNY=S)',
1080 $ dumma( 2 ), n, jtype, ioldsd
1087 result( 12 ) = ulpinv
1088 CALL clacpy(
'F', n, n, z, ldu, evectr, ldu )
1089 CALL ctgevc(
'R',
'B', llwork, n, s1, lda, p1, lda, cdumma,
1090 $ ldu, evectr, ldu, n, in, work, rwork, iinfo )
1091 IF( iinfo.NE.0 )
THEN
1092 WRITE( nounit, fmt = 9999 )
'CTGEVC(R,B)', iinfo, n,
1098 CALL cget52( .false., n, h, lda, t, lda, evectr, ldu,
1099 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1100 result( 12 ) = dumma( 1 )
1101 IF( dumma( 2 ).GT.thresh )
THEN
1102 WRITE( nounit, fmt = 9998 )
'Right',
'CTGEVC(HOWMNY=B)',
1103 $ dumma( 2 ), n, jtype, ioldsd
1112 CALL cget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1113 $ work, rwork, result( 13 ) )
1114 CALL cget51( 2, n, p1, lda, p2, lda, q, ldu, z, ldu,
1115 $ work, rwork, result( 14 ) )
1122 temp1 = max( temp1, abs( alpha1( j )-alpha3( j ) ) )
1123 temp2 = max( temp2, abs( beta1( j )-beta3( j ) ) )
1126 temp1 = temp1 / max( safmin, ulp*max( temp1, anorm ) )
1127 temp2 = temp2 / max( safmin, ulp*max( temp2, bnorm ) )
1128 result( 15 ) = max( temp1, temp2 )
1141 ntestt = ntestt + ntest
1145 DO 220 jr = 1, ntest
1146 IF( result( jr ).GE.thresh )
THEN
1151 IF( nerrs.EQ.0 )
THEN
1152 WRITE( nounit, fmt = 9997 )
'CGG'
1156 WRITE( nounit, fmt = 9996 )
1157 WRITE( nounit, fmt = 9995 )
1158 WRITE( nounit, fmt = 9994 )
'Unitary'
1162 WRITE( nounit, fmt = 9993 )
'unitary',
'*',
1163 $
'conjugate transpose', (
'*', j = 1, 10 )
1167 IF( result( jr ).LT.10000.0 )
THEN
1168 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1171 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1182 CALL slasum(
'CGG', nounit, nerrs, ntestt )
1185 9999
FORMAT(
' CCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1186 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1188 9998
FORMAT(
' CCHKGG: ', a,
' Eigenvectors from ', a,
' incorrectly ',
1189 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
1190 $
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
1193 9997
FORMAT( 1x, a3,
' -- Complex Generalized eigenvalue problem' )
1195 9996
FORMAT(
' Matrix types (see CCHKGG for details): ' )
1197 9995
FORMAT(
' Special Matrices:', 23x,
1198 $
'(J''=transposed Jordan block)',
1199 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
1200 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
1201 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
1202 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
1203 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1204 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
1205 9994
FORMAT(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
1206 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
1207 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
1208 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
1209 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
1210 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
1211 $
'23=(small,large) 24=(small,small) 25=(large,large)',
1212 $ /
' 26=random O(1) matrices.' )
1214 9993
FORMAT( /
' Tests performed: (H is Hessenberg, S is Schur, B, ',
1215 $
'T, P are triangular,', / 20x,
'U, V, Q, and Z are ', a,
1216 $
', l and r are the', / 20x,
1217 $
'appropriate left and right eigenvectors, resp., a is',
1218 $ / 20x,
'alpha, b is beta, and ', a,
' means ', a,
'.)',
1219 $ /
' 1 = | A - U H V', a,
1220 $
' | / ( |A| n ulp ) 2 = | B - U T V', a,
1221 $
' | / ( |B| n ulp )', /
' 3 = | I - UU', a,
1222 $
' | / ( n ulp ) 4 = | I - VV', a,
1223 $
' | / ( n ulp )', /
' 5 = | H - Q S Z', a,
1224 $
' | / ( |H| n ulp )', 6x,
'6 = | T - Q P Z', a,
1225 $
' | / ( |T| n ulp )', /
' 7 = | I - QQ', a,
1226 $
' | / ( n ulp ) 8 = | I - ZZ', a,
1227 $
' | / ( n ulp )', /
' 9 = max | ( b S - a P )', a,
1228 $
' l | / const. 10 = max | ( b H - a T )', a,
1229 $
' l | / const.', /
1230 $
' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
1231 $
' - a T ) r | / const.', / 1x )
1233 9992
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
1234 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
1235 9991
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
1236 $ 4( i4,
',' ),
' result ', i2,
' is', 1p, e10.3 )