500 SUBROUTINE cchkgg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
501 $ tstdif, thrshn, nounit, a, lda, b, h, t, s1,
502 $ s2, p1, p2, u, ldu, v, q, z, alpha1, beta1,
503 $ alpha3, beta3, evectl, evectr, work, lwork,
504 $ rwork, llwork, result, info )
513 INTEGER info, lda, ldu, lwork, nounit, nsizes, ntypes
517 LOGICAL dotype( * ), llwork( * )
518 INTEGER iseed( 4 ), nn( * )
519 REAL result( 15 ), rwork( * )
520 COMPLEX a( lda, * ), alpha1( * ), alpha3( * ),
521 $ b( lda, * ), beta1( * ), beta3( * ),
522 $ evectl( ldu, * ), evectr( ldu, * ),
523 $ h( lda, * ), p1( lda, * ), p2( lda, * ),
524 $ q( ldu, * ), s1( lda, * ), s2( lda, * ),
525 $ t( lda, * ), u( ldu, * ), v( ldu, * ),
526 $ work( * ), z( ldu, * )
533 parameter( zero = 0.0e+0, one = 1.0e+0 )
535 parameter( czero = ( 0.0e+0, 0.0e+0 ),
536 $ cone = ( 1.0e+0, 0.0e+0 ) )
538 parameter( maxtyp = 26 )
542 INTEGER i1, iadd, iinfo, in, j, jc, jr, jsize, jtype,
543 $ lwkopt, mtypes, n, n1, nerrs, nmats, nmax,
545 REAL anorm, bnorm, safmax, safmin, temp1, temp2,
550 LOGICAL lasign( maxtyp ), lbsign( maxtyp )
551 INTEGER ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
552 $ katype( maxtyp ), kazero( maxtyp ),
553 $ kbmagn( maxtyp ), kbtype( maxtyp ),
554 $ kbzero( maxtyp ), kclass( maxtyp ),
555 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
556 REAL dumma( 4 ), rmagn( 0: 3 )
570 INTRINSIC abs, conjg, max, min,
REAL, sign
573 DATA kclass / 15*1, 10*2, 1*3 /
574 DATA kz1 / 0, 1, 2, 1, 3, 3 /
575 DATA kz2 / 0, 0, 1, 2, 1, 1 /
576 DATA kadd / 0, 0, 0, 0, 3, 2 /
577 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
578 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
579 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
580 $ 1, 1, -4, 2, -4, 8*8, 0 /
581 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
583 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
585 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
587 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
589 DATA ktrian / 16*0, 10*1 /
590 DATA lasign / 6*.false., .true., .false., 2*.true.,
591 $ 2*.false., 3*.true., .false., .true.,
592 $ 3*.false., 5*.true., .false. /
593 DATA lbsign / 7*.false., .true., 2*.false.,
594 $ 2*.true., 2*.false., .true., .false., .true.,
606 nmax = max( nmax, nn( j ) )
611 lwkopt = max( 2*nmax*nmax, 4*nmax, 1 )
615 IF( nsizes.LT.0 )
THEN
617 ELSE IF( badnn )
THEN
619 ELSE IF( ntypes.LT.0 )
THEN
621 ELSE IF( thresh.LT.zero )
THEN
623 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
625 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax )
THEN
627 ELSE IF( lwkopt.GT.lwork )
THEN
632 CALL
xerbla(
'CCHKGG', -info )
638 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
641 safmin =
slamch(
'Safe minimum' )
643 safmin = safmin / ulp
644 safmax = one / safmin
645 CALL
slabad( safmin, safmax )
659 DO 240 jsize = 1, nsizes
662 rmagn( 2 ) = safmax*ulp /
REAL( n1 )
663 rmagn( 3 ) = safmin*ulpinv*n1
665 IF( nsizes.NE.1 )
THEN
666 mtypes = min( maxtyp, ntypes )
668 mtypes = min( maxtyp+1, ntypes )
671 DO 230 jtype = 1, mtypes
672 IF( .NOT.dotype( jtype ) )
680 ioldsd( j ) = iseed( j )
710 IF( mtypes.GT.maxtyp )
713 IF( kclass( jtype ).LT.3 )
THEN
717 IF( abs( katype( jtype ) ).EQ.3 )
THEN
718 in = 2*( ( n-1 ) / 2 ) + 1
720 $ CALL
claset(
'Full', n, n, czero, czero, a, lda )
724 CALL
clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
725 $ kz2( kazero( jtype ) ), lasign( jtype ),
726 $ rmagn( kamagn( jtype ) ), ulp,
727 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 4,
729 iadd = kadd( kazero( jtype ) )
730 IF( iadd.GT.0 .AND. iadd.LE.n )
731 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
735 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
736 in = 2*( ( n-1 ) / 2 ) + 1
738 $ CALL
claset(
'Full', n, n, czero, czero, b, lda )
742 CALL
clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
743 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
744 $ rmagn( kbmagn( jtype ) ), one,
745 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 4,
747 iadd = kadd( kbzero( jtype ) )
749 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
751 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
761 u( jr, jc ) =
clarnd( 3, iseed )
762 v( jr, jc ) =
clarnd( 3, iseed )
764 CALL
clarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
766 work( 2*n+jc ) = sign( one,
REAL( U( JC, JC ) ) )
768 CALL
clarfg( n+1-jc, v( jc, jc ), v( jc+1, jc ), 1,
770 work( 3*n+jc ) = sign( one,
REAL( V( JC, JC ) ) )
773 ctemp =
clarnd( 3, iseed )
776 work( 3*n ) = ctemp / abs( ctemp )
777 ctemp =
clarnd( 3, iseed )
780 work( 4*n ) = ctemp / abs( ctemp )
786 a( jr, jc ) = work( 2*n+jr )*
787 $ conjg( work( 3*n+jc ) )*
789 b( jr, jc ) = work( 2*n+jr )*
790 $ conjg( work( 3*n+jc ) )*
794 CALL
cunm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
795 $ lda, work( 2*n+1 ), iinfo )
798 CALL
cunm2r(
'R',
'C', n, n, n-1, v, ldu, work( n+1 ),
799 $ a, lda, work( 2*n+1 ), iinfo )
802 CALL
cunm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
803 $ lda, work( 2*n+1 ), iinfo )
806 CALL
cunm2r(
'R',
'C', n, n, n-1, v, ldu, work( n+1 ),
807 $ b, lda, work( 2*n+1 ), iinfo )
817 a( jr, jc ) = rmagn( kamagn( jtype ) )*
819 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
825 anorm =
clange(
'1', n, n, a, lda, rwork )
826 bnorm =
clange(
'1', n, n, b, lda, rwork )
830 IF( iinfo.NE.0 )
THEN
831 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
841 CALL
clacpy(
' ', n, n, a, lda, h, lda )
842 CALL
clacpy(
' ', n, n, b, lda, t, lda )
846 CALL
cgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
847 IF( iinfo.NE.0 )
THEN
848 WRITE( nounit, fmt = 9999 )
'CGEQR2', iinfo, n, jtype,
854 CALL
cunm2r(
'L',
'C', n, n, n, t, lda, work, h, lda,
855 $ work( n+1 ), iinfo )
856 IF( iinfo.NE.0 )
THEN
857 WRITE( nounit, fmt = 9999 )
'CUNM2R', iinfo, n, jtype,
863 CALL
claset(
'Full', n, n, czero, cone, u, ldu )
864 CALL
cunm2r(
'R',
'N', n, n, n, t, lda, work, u, ldu,
865 $ work( n+1 ), iinfo )
866 IF( iinfo.NE.0 )
THEN
867 WRITE( nounit, fmt = 9999 )
'CUNM2R', iinfo, n, jtype,
873 CALL
cgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
875 IF( iinfo.NE.0 )
THEN
876 WRITE( nounit, fmt = 9999 )
'CGGHRD', iinfo, n, jtype,
885 CALL
cget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
886 $ rwork, result( 1 ) )
887 CALL
cget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
888 $ rwork, result( 2 ) )
889 CALL
cget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
890 $ rwork, result( 3 ) )
891 CALL
cget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
892 $ rwork, result( 4 ) )
900 CALL
clacpy(
' ', n, n, h, lda, s2, lda )
901 CALL
clacpy(
' ', n, n, t, lda, p2, lda )
905 CALL
chgeqz(
'E',
'N',
'N', n, 1, n, s2, lda, p2, lda,
906 $ alpha3, beta3, q, ldu, z, ldu, work, lwork,
908 IF( iinfo.NE.0 )
THEN
909 WRITE( nounit, fmt = 9999 )
'CHGEQZ(E)', iinfo, n, jtype,
917 CALL
clacpy(
' ', n, n, h, lda, s2, lda )
918 CALL
clacpy(
' ', n, n, t, lda, p2, lda )
920 CALL
chgeqz(
'S',
'N',
'N', n, 1, n, s2, lda, p2, lda,
921 $ alpha1, beta1, q, ldu, z, ldu, work, lwork,
923 IF( iinfo.NE.0 )
THEN
924 WRITE( nounit, fmt = 9999 )
'CHGEQZ(S)', iinfo, n, jtype,
932 CALL
clacpy(
' ', n, n, h, lda, s1, lda )
933 CALL
clacpy(
' ', n, n, t, lda, p1, lda )
935 CALL
chgeqz(
'S',
'I',
'I', n, 1, n, s1, lda, p1, lda,
936 $ alpha1, beta1, q, ldu, z, ldu, work, lwork,
938 IF( iinfo.NE.0 )
THEN
939 WRITE( nounit, fmt = 9999 )
'CHGEQZ(V)', iinfo, n, jtype,
949 CALL
cget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
950 $ rwork, result( 5 ) )
951 CALL
cget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
952 $ rwork, result( 6 ) )
953 CALL
cget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
954 $ rwork, result( 7 ) )
955 CALL
cget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
956 $ rwork, result( 8 ) )
974 llwork( j ) = .false.
977 CALL
ctgevc(
'L',
'S', llwork, n, s1, lda, p1, lda, evectl,
978 $ ldu, cdumma, ldu, n, in, work, rwork, iinfo )
979 IF( iinfo.NE.0 )
THEN
980 WRITE( nounit, fmt = 9999 )
'CTGEVC(L,S1)', iinfo, n,
988 llwork( j ) = .false.
994 CALL
ctgevc(
'L',
'S', llwork, n, s1, lda, p1, lda,
995 $ evectl( 1, i1+1 ), ldu, cdumma, ldu, n, in,
996 $ work, rwork, iinfo )
997 IF( iinfo.NE.0 )
THEN
998 WRITE( nounit, fmt = 9999 )
'CTGEVC(L,S2)', iinfo, n,
1004 CALL
cget52( .true., n, s1, lda, p1, lda, evectl, ldu,
1005 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1006 result( 9 ) = dumma( 1 )
1007 IF( dumma( 2 ).GT.thrshn )
THEN
1008 WRITE( nounit, fmt = 9998 )
'Left',
'CTGEVC(HOWMNY=S)',
1009 $ dumma( 2 ), n, jtype, ioldsd
1016 result( 10 ) = ulpinv
1017 CALL
clacpy(
'F', n, n, q, ldu, evectl, ldu )
1018 CALL
ctgevc(
'L',
'B', llwork, n, s1, lda, p1, lda, evectl,
1019 $ ldu, cdumma, ldu, n, in, work, rwork, iinfo )
1020 IF( iinfo.NE.0 )
THEN
1021 WRITE( nounit, fmt = 9999 )
'CTGEVC(L,B)', iinfo, n,
1027 CALL
cget52( .true., n, h, lda, t, lda, evectl, ldu, alpha1,
1028 $ beta1, work, rwork, dumma( 1 ) )
1029 result( 10 ) = dumma( 1 )
1030 IF( dumma( 2 ).GT.thrshn )
THEN
1031 WRITE( nounit, fmt = 9998 )
'Left',
'CTGEVC(HOWMNY=B)',
1032 $ dumma( 2 ), n, jtype, ioldsd
1039 result( 11 ) = ulpinv
1046 llwork( j ) = .true.
1048 DO 170 j = i1 + 1, n
1049 llwork( j ) = .false.
1052 CALL
ctgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, cdumma,
1053 $ ldu, evectr, ldu, n, in, work, rwork, iinfo )
1054 IF( iinfo.NE.0 )
THEN
1055 WRITE( nounit, fmt = 9999 )
'CTGEVC(R,S1)', iinfo, n,
1063 llwork( j ) = .false.
1065 DO 190 j = i1 + 1, n
1066 llwork( j ) = .true.
1069 CALL
ctgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, cdumma,
1070 $ ldu, evectr( 1, i1+1 ), ldu, n, in, work,
1072 IF( iinfo.NE.0 )
THEN
1073 WRITE( nounit, fmt = 9999 )
'CTGEVC(R,S2)', iinfo, n,
1079 CALL
cget52( .false., n, s1, lda, p1, lda, evectr, ldu,
1080 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1081 result( 11 ) = dumma( 1 )
1082 IF( dumma( 2 ).GT.thresh )
THEN
1083 WRITE( nounit, fmt = 9998 )
'Right',
'CTGEVC(HOWMNY=S)',
1084 $ dumma( 2 ), n, jtype, ioldsd
1091 result( 12 ) = ulpinv
1092 CALL
clacpy(
'F', n, n, z, ldu, evectr, ldu )
1093 CALL
ctgevc(
'R',
'B', llwork, n, s1, lda, p1, lda, cdumma,
1094 $ ldu, evectr, ldu, n, in, work, rwork, iinfo )
1095 IF( iinfo.NE.0 )
THEN
1096 WRITE( nounit, fmt = 9999 )
'CTGEVC(R,B)', iinfo, n,
1102 CALL
cget52( .false., n, h, lda, t, lda, evectr, ldu,
1103 $ alpha1, beta1, work, rwork, dumma( 1 ) )
1104 result( 12 ) = dumma( 1 )
1105 IF( dumma( 2 ).GT.thresh )
THEN
1106 WRITE( nounit, fmt = 9998 )
'Right',
'CTGEVC(HOWMNY=B)',
1107 $ dumma( 2 ), n, jtype, ioldsd
1116 CALL
cget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1117 $ work, rwork, result( 13 ) )
1118 CALL
cget51( 2, n, p1, lda, p2, lda, q, ldu, z, ldu,
1119 $ work, rwork, result( 14 ) )
1126 temp1 = max( temp1, abs( alpha1( j )-alpha3( j ) ) )
1127 temp2 = max( temp2, abs( beta1( j )-beta3( j ) ) )
1130 temp1 = temp1 / max( safmin, ulp*max( temp1, anorm ) )
1131 temp2 = temp2 / max( safmin, ulp*max( temp2, bnorm ) )
1132 result( 15 ) = max( temp1, temp2 )
1145 ntestt = ntestt + ntest
1149 DO 220 jr = 1, ntest
1150 IF( result( jr ).GE.thresh )
THEN
1155 IF( nerrs.EQ.0 )
THEN
1156 WRITE( nounit, fmt = 9997 )
'CGG'
1160 WRITE( nounit, fmt = 9996 )
1161 WRITE( nounit, fmt = 9995 )
1162 WRITE( nounit, fmt = 9994 )
'Unitary'
1166 WRITE( nounit, fmt = 9993 )
'unitary',
'*',
1167 $
'conjugate transpose', (
'*', j = 1, 10 )
1171 IF( result( jr ).LT.10000.0 )
THEN
1172 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1175 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1186 CALL
slasum(
'CGG', nounit, nerrs, ntestt )
1189 9999 format(
' CCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1190 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1192 9998 format(
' CCHKGG: ', a,
' Eigenvectors from ', a,
' incorrectly ',
1193 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
1194 $
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
1197 9997 format( 1x, a3,
' -- Complex Generalized eigenvalue problem' )
1199 9996 format(
' Matrix types (see CCHKGG for details): ' )
1201 9995 format(
' Special Matrices:', 23x,
1202 $
'(J''=transposed Jordan block)',
1203 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
1204 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
1205 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
1206 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
1207 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1208 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
1209 9994 format(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
1210 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
1211 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
1212 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
1213 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
1214 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
1215 $
'23=(small,large) 24=(small,small) 25=(large,large)',
1216 $ /
' 26=random O(1) matrices.' )
1218 9993 format( /
' Tests performed: (H is Hessenberg, S is Schur, B, ',
1219 $
'T, P are triangular,', / 20x,
'U, V, Q, and Z are ', a,
1220 $
', l and r are the', / 20x,
1221 $
'appropriate left and right eigenvectors, resp., a is',
1222 $ / 20x,
'alpha, b is beta, and ', a,
' means ', a,
'.)',
1223 $ /
' 1 = | A - U H V', a,
1224 $
' | / ( |A| n ulp ) 2 = | B - U T V', a,
1225 $
' | / ( |B| n ulp )', /
' 3 = | I - UU', a,
1226 $
' | / ( n ulp ) 4 = | I - VV', a,
1227 $
' | / ( n ulp )', /
' 5 = | H - Q S Z', a,
1228 $
' | / ( |H| n ulp )', 6x,
'6 = | T - Q P Z', a,
1229 $
' | / ( |T| n ulp )', /
' 7 = | I - QQ', a,
1230 $
' | / ( n ulp ) 8 = | I - ZZ', a,
1231 $
' | / ( n ulp )', /
' 9 = max | ( b S - a P )', a,
1232 $
' l | / const. 10 = max | ( b H - a T )', a,
1233 $
' l | / const.', /
1234 $
' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
1235 $
' - a T ) r | / const.', / 1x )
1237 9992 format(
' Matrix order=', i5,
', type=', i2,
', seed=',
1238 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
1239 9991 format(
' Matrix order=', i5,
', type=', i2,
', seed=',
1240 $ 4( i4,
',' ),
' result ', i2,
' is', 1p, e10.3 )