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 ( * ), LLWORK( * )
515 INTEGER ( 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, * ),
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 ), ( 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
642 CALL slabad( safmin, safmax )
656 DO 240 jsize = 1, nsizes
659 rmagn( 2 ) = safmax*ulp / real( 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 ) )
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 claset(
'Full', n, n, czero, czero, a, lda )
721 CALL clatm4( 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 IF( iadd.GT.0 .AND. iadd.LE.n )
728 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
732 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
733 in = 2*( ( n-1 ) / 2 ) + 1
735 $
CALL claset(
'Full', n, n, czero, czero, b, lda )
739 CALL clatm4( 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 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
758 u( jr, jc ) = clarnd( 3, iseed )
759 v( jr, jc ) = clarnd( 3, iseed )
761 CALL clarfg( n+1-jc, u( jc, jc ), u( jc+1, jc
763 work( 2*n+jc ) = sign( one, real( u( jc, jc ) ) )
765 CALL clarfg( n+1-jc, v( jc, jc ), v( jc+1, jc ), 1,
767 work( 3*n+jc ) = sign( one, real( v( jc, jc ) ) )
770 ctemp = clarnd( 3, iseed )
773 work( 3*n ) = ctemp / abs( ctemp )
774 ctemp = clarnd( 3, iseed )
777 work( 4*n ) = ctemp / abs( ctemp )
783 a( jr, jc ) = work( 2*n+jr )*
784 $ conjg( work( 3*n+jc ) )*
786 b( jr, jc ) = work( 2*n+jr )*
787 $ conjg( work( 3*n+jc ) )*
791 CALL cunm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
792 $ lda, work( 2*n+1 ), iinfo )
795 CALL cunm2r(
'R',
'C', n, n, n-1, v, ldu, work( n+1 ),
796 $ a, lda, work( 2*n+1 ), iinfo )
799 CALL cunm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
800 $ lda, work( 2*n+1 ), iinfo )
803 CALL cunm2r(
'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 = clange(
'1', n, n, a, lda, rwork )
823 bnorm = clange(
'1', n, n, b, lda, rwork )
827 IF( iinfo.NE.0 )
THEN
828 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
838 CALL clacpy(
' ', n, n, a, lda, h, lda )
839 CALL clacpy(
' ', n, n, b, lda, t, lda )
844 IF( iinfo.NE.0 )
THEN
845 WRITE( nounit, fmt = 9999 )
'CGEQR2', iinfo, n, jtype,
851 CALL cunm2r(
'L',
'C', n, n, n, t, lda, work, h, lda,
852 $ work( n+1 ), iinfo )
853 IF( iinfo.NE.0 )
THEN
854 WRITE( nounit, fmt = 9999 )
'CUNM2R', iinfo, n, jtype,
860 CALL claset(
'Full', n, n, czero, cone, u, ldu )
861 CALL cunm2r(
'R',
'N', n, n, n, t, lda, work, u, ldu
862 $ work( n+1 ), iinfo )
863 IF( iinfo.NE.0 )
THEN
864 WRITE( nounit, fmt = 9999 )
'CUNM2R', iinfo, n, jtype,
870 CALL cgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
872 IF( iinfo.NE.0 )
THEN
873 WRITE( nounit, fmt = 9999 )
'CGGHRD', iinfo, n, jtype,
882 CALL cget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
883 $ rwork, result( 1 ) )
884 CALL cget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
885 $ rwork, result( 2 ) )
886 CALL cget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
887 $ rwork, result( 3 ) )
888 CALL cget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
889 $ rwork, result( 4 ) )
897 CALL clacpy(
' ', n, n, h, lda, s2, lda )
898 CALL clacpy(
' ', n, n, t, lda, p2, lda )
902 CALL chgeqz(
'E',
'N',
'N', n, 1, n, s2, lda, p2, lda,
903 $ alpha3, beta3, q, ldu, z, ldu, work, lwork,
905 IF( iinfo.NE.0 )
THEN
906 WRITE( nounit, fmt = 9999 )
'CHGEQZ(E)', iinfo, n, jtype,
914 CALL clacpy(
' ', n, n, h, lda, s2, lda )
915 CALL clacpy(
' ', n, n, t, lda, p2, lda )
917 CALL chgeqz(
'S',
'N',
'N', n, 1, n, s2, lda, p2, lda,
918 $ alpha1, beta1, q, ldu, z, ldu, work, lwork,
920 IF( iinfo.NE.0 )
THEN
921 WRITE( nounit, fmt = 9999 )
'CHGEQZ(S)', iinfo, n, jtype,
929 CALL clacpy(
' ', n, n, h, lda, s1, lda )
930 CALL clacpy(
' ', n, n, t, lda, p1, lda )
932 CALL chgeqz(
'S',
'I',
'I', n, 1, n, s1, lda, p1, lda,
933 $ alpha1, beta1, q, ldu, z, ldu, work, lwork,
935 IF( iinfo.NE.0 )
THEN
936 WRITE( nounit, fmt = 9999 )
'CHGEQZ(V)', iinfo, n, jtype,
946 CALL cget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
947 $ rwork, result( 5 ) )
948 CALL cget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
949 $ rwork, result( 6 ) )
950 CALL cget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
951 $ rwork, result( 7 ) )
952 CALL cget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
953 $ rwork, result( 8 ) )
971 llwork( j ) = .false.
974 CALL ctgevc(
'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 )'ctgevc(l,s1)
', IINFO, N,
985 LLWORK( J ) = .FALSE.
991 CALL CTGEVC( 'l
', 's
', LLWORK, N, S1, LDA, P1, LDA,
992 $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN,
993 $ WORK, RWORK, IINFO )
994.NE.
IF( IINFO0 ) THEN
995 WRITE( NOUNIT, FMT = 9999 )'ctgevc(l,s2)
', IINFO, N,
1001 CALL CGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU,
1002 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
1003 RESULT( 9 ) = DUMMA( 1 )
1004.GT.
IF( DUMMA( 2 )THRSHN ) THEN
1005 WRITE( NOUNIT, FMT = 9998 )'left
', 'ctgevc(howmny=s)
',
1006 $ DUMMA( 2 ), N, JTYPE, IOLDSD
1013 RESULT( 10 ) = ULPINV
1014 CALL CLACPY( 'f
', N, N, Q, LDU, EVECTL, LDU )
1015 CALL CTGEVC( 'l
', 'b
', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
1016 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
1017.NE.
IF( IINFO0 ) THEN
1018 WRITE( NOUNIT, FMT = 9999 )'ctgevc(l,b)
', IINFO, N,
1024 CALL CGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1,
1025 $ BETA1, WORK, RWORK, DUMMA( 1 ) )
1026 RESULT( 10 ) = DUMMA( 1 )
1027.GT.
IF( DUMMA( 2 )THRSHN ) THEN
1028 WRITE( NOUNIT, FMT = 9998 )'left
', 'ctgevc(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 CTGEVC( 'r
', 's', llwork, n, s1, lda, p1, lda, cdumma,
1050 $ ldu, evectr, ldu, n, in, work, rwork, iinfo )
1051 IF( iinfo.NE.0 )
THEN
1052 WRITE( nounit, fmt = 9999 )
'CTGEVC(R,S1)', iinfo, n,
1060 llwork( j ) = .false.
1062 DO 190 j = i1 + 1, n
1066 CALL ctgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, cdumma,
1067 $ ldu, evectr( 1, i1+1 ), ldu, n, in, work,
1069 IF( iinfo.NE.0 )
THEN
1070 WRITE'CTGEVC(R,S2)', iinfo, n,
1076 CALL cget52( .false., n, s1, lda, p1, lda, evectr, ldu,
1077 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1078 result( 11 ) = dumma( 1 )
1079 IF( dumma( 2 ).GT.thresh )
THEN
1080 WRITE( nounit, fmt = 9998 )
'Right',
'CTGEVC(HOWMNY=S)',
1081 $ dumma( 2 ), n, jtype, ioldsd
1088 result( 12 ) = ulpinv
1089 CALL clacpy(
'F', n, n, z, ldu, evectr, ldu )
1090 CALL ctgevc(
'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 )
'CTGEVC(R,B)', iinfo, n,
1099 CALL cget52( .false., n, h, lda, t, lda, evectr, ldu,
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',
'CTGEVC(HOWMNY=B)',
1104 $ dumma( 2 ), n, jtype, ioldsd
1113 CALL cget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1114 $ work, rwork, result( 13 ) )
1115 CALL cget51( 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 )
'CGG'
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.0 )
THEN
1169 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1172 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1183 CALL slasum(
'CGG', nounit, nerrs, ntestt )
1186 9999
FORMAT(
' CCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1187 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1189 9998
FORMAT(
' CCHKGG: ', a,
' Eigenvectors from ', a,
' incorrectly ',
1190 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
1191 $
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
1194 9997
FORMAT( 1x, a3,
' -- Complex Generalized eigenvalue problem' )
1196 9996
FORMAT(
' Matrix types (see CCHKGG 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, e10.3 )