508 SUBROUTINE dchkgg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
509 $ tstdif, thrshn, nounit, a, lda, b, h, t, s1,
510 $ s2, p1, p2, u, ldu, v, q, z, alphr1, alphi1,
511 $ beta1, alphr3, alphi3, beta3, evectl, evectr,
512 $ work, lwork, llwork, result, info )
521 INTEGER info, lda, ldu, lwork, nounit, nsizes, ntypes
522 DOUBLE PRECISION thresh, thrshn
525 LOGICAL dotype( * ), llwork( * )
526 INTEGER iseed( 4 ), nn( * )
527 DOUBLE PRECISION a( lda, * ), alphi1( * ), alphi3( * ),
528 $ alphr1( * ), alphr3( * ), b( lda, * ),
529 $ beta1( * ), beta3( * ), evectl( ldu, * ),
530 $ evectr( ldu, * ), h( lda, * ), p1( lda, * ),
531 $ p2( lda, * ), q( ldu, * ), result( 15 ),
532 $ s1( lda, * ), s2( lda, * ), t( lda, * ),
533 $ u( ldu, * ), v( ldu, * ), work( * ),
540 DOUBLE PRECISION zero, one
541 parameter( zero = 0.0d0, one = 1.0d0 )
543 parameter( maxtyp = 26 )
547 INTEGER i1, iadd, iinfo, in, j, jc, jr, jsize, jtype,
548 $ lwkopt, mtypes, n, n1, nerrs, nmats, nmax,
550 DOUBLE PRECISION anorm, bnorm, safmax, safmin, temp1, temp2,
554 INTEGER iasign( maxtyp ), ibsign( maxtyp ),
555 $ ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
556 $ katype( maxtyp ), kazero( maxtyp ),
557 $ kbmagn( maxtyp ), kbtype( maxtyp ),
558 $ kbzero( maxtyp ), kclass( maxtyp ),
559 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
560 DOUBLE PRECISION dumma( 4 ), rmagn( 0: 3 )
572 INTRINSIC abs, dble, max, min, sign
575 DATA kclass / 15*1, 10*2, 1*3 /
576 DATA kz1 / 0, 1, 2, 1, 3, 3 /
577 DATA kz2 / 0, 0, 1, 2, 1, 1 /
578 DATA kadd / 0, 0, 0, 0, 3, 2 /
579 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
580 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
581 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
582 $ 1, 1, -4, 2, -4, 8*8, 0 /
583 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
585 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
587 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
589 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
591 DATA ktrian / 16*0, 10*1 /
592 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
594 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
605 nmax = max( nmax, nn( j ) )
613 lwkopt = max( 6*nmax, 2*nmax*nmax, 1 )
617 IF( nsizes.LT.0 )
THEN
619 ELSE IF( badnn )
THEN
621 ELSE IF( ntypes.LT.0 )
THEN
623 ELSE IF( thresh.LT.zero )
THEN
625 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
627 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax )
THEN
629 ELSE IF( lwkopt.GT.lwork )
THEN
634 CALL
xerbla(
'DCHKGG', -info )
640 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
643 safmin =
dlamch(
'Safe minimum' )
645 safmin = safmin / ulp
646 safmax = one / safmin
647 CALL
dlabad( safmin, safmax )
661 DO 240 jsize = 1, nsizes
664 rmagn( 2 ) = safmax*ulp / dble( n1 )
665 rmagn( 3 ) = safmin*ulpinv*n1
667 IF( nsizes.NE.1 )
THEN
668 mtypes = min( maxtyp, ntypes )
670 mtypes = min( maxtyp+1, ntypes )
673 DO 230 jtype = 1, mtypes
674 IF( .NOT.dotype( jtype ) )
682 ioldsd( j ) = iseed( j )
714 IF( mtypes.GT.maxtyp )
717 IF( kclass( jtype ).LT.3 )
THEN
721 IF( abs( katype( jtype ) ).EQ.3 )
THEN
722 in = 2*( ( n-1 ) / 2 ) + 1
724 $ CALL
dlaset(
'Full', n, n, zero, zero, a, lda )
728 CALL
dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
729 $ kz2( kazero( jtype ) ), iasign( jtype ),
730 $ rmagn( kamagn( jtype ) ), ulp,
731 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
733 iadd = kadd( kazero( jtype ) )
734 IF( iadd.GT.0 .AND. iadd.LE.n )
735 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
739 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
740 in = 2*( ( n-1 ) / 2 ) + 1
742 $ CALL
dlaset(
'Full', n, n, zero, zero, b, lda )
746 CALL
dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
747 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
748 $ rmagn( kbmagn( jtype ) ), one,
749 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
751 iadd = kadd( kbzero( jtype ) )
752 IF( iadd.NE.0 .AND. iadd.LE.n )
753 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
755 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
764 u( jr, jc ) =
dlarnd( 3, iseed )
765 v( jr, jc ) =
dlarnd( 3, iseed )
767 CALL
dlarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
769 work( 2*n+jc ) = sign( one, u( jc, jc ) )
771 CALL
dlarfg( n+1-jc, v( jc, jc ), v( jc+1, jc ), 1,
773 work( 3*n+jc ) = sign( one, v( jc, jc ) )
778 work( 3*n ) = sign( one,
dlarnd( 2, iseed ) )
781 work( 4*n ) = sign( one,
dlarnd( 2, iseed ) )
787 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
789 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
793 CALL
dorm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
794 $ lda, work( 2*n+1 ), iinfo )
797 CALL
dorm2r(
'R',
'T', n, n, n-1, v, ldu, work( n+1 ),
798 $ a, lda, work( 2*n+1 ), iinfo )
801 CALL
dorm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
802 $ lda, work( 2*n+1 ), iinfo )
805 CALL
dorm2r(
'R',
'T', n, n, n-1, v, ldu, work( n+1 ),
806 $ b, lda, work( 2*n+1 ), iinfo )
816 a( jr, jc ) = rmagn( kamagn( jtype ) )*
818 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
824 anorm =
dlange(
'1', n, n, a, lda, work )
825 bnorm =
dlange(
'1', n, n, b, lda, work )
829 IF( iinfo.NE.0 )
THEN
830 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
840 CALL
dlacpy(
' ', n, n, a, lda, h, lda )
841 CALL
dlacpy(
' ', n, n, b, lda, t, lda )
845 CALL
dgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
846 IF( iinfo.NE.0 )
THEN
847 WRITE( nounit, fmt = 9999 )
'DGEQR2', iinfo, n, jtype,
853 CALL
dorm2r(
'L',
'T', n, n, n, t, lda, work, h, lda,
854 $ work( n+1 ), iinfo )
855 IF( iinfo.NE.0 )
THEN
856 WRITE( nounit, fmt = 9999 )
'DORM2R', iinfo, n, jtype,
862 CALL
dlaset(
'Full', n, n, zero, one, u, ldu )
863 CALL
dorm2r(
'R',
'N', n, n, n, t, lda, work, u, ldu,
864 $ work( n+1 ), iinfo )
865 IF( iinfo.NE.0 )
THEN
866 WRITE( nounit, fmt = 9999 )
'DORM2R', iinfo, n, jtype,
872 CALL
dgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
874 IF( iinfo.NE.0 )
THEN
875 WRITE( nounit, fmt = 9999 )
'DGGHRD', iinfo, n, jtype,
884 CALL
dget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
886 CALL
dget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
888 CALL
dget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
890 CALL
dget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
899 CALL
dlacpy(
' ', n, n, h, lda, s2, lda )
900 CALL
dlacpy(
' ', n, n, t, lda, p2, lda )
904 CALL
dhgeqz(
'E',
'N',
'N', n, 1, n, s2, lda, p2, lda,
905 $ alphr3, alphi3, beta3, q, ldu, z, ldu, work,
907 IF( iinfo.NE.0 )
THEN
908 WRITE( nounit, fmt = 9999 )
'DHGEQZ(E)', iinfo, n, jtype,
916 CALL
dlacpy(
' ', n, n, h, lda, s2, lda )
917 CALL
dlacpy(
' ', n, n, t, lda, p2, lda )
919 CALL
dhgeqz(
'S',
'N',
'N', n, 1, n, s2, lda, p2, lda,
920 $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
922 IF( iinfo.NE.0 )
THEN
923 WRITE( nounit, fmt = 9999 )
'DHGEQZ(S)', iinfo, n, jtype,
931 CALL
dlacpy(
' ', n, n, h, lda, s1, lda )
932 CALL
dlacpy(
' ', n, n, t, lda, p1, lda )
934 CALL
dhgeqz(
'S',
'I',
'I', n, 1, n, s1, lda, p1, lda,
935 $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
937 IF( iinfo.NE.0 )
THEN
938 WRITE( nounit, fmt = 9999 )
'DHGEQZ(V)', iinfo, n, jtype,
948 CALL
dget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
950 CALL
dget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
952 CALL
dget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
954 CALL
dget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
973 llwork( j ) = .false.
976 CALL
dtgevc(
'L',
'S', llwork, n, s1, lda, p1, lda, evectl,
977 $ ldu, dumma, ldu, n, in, work, iinfo )
978 IF( iinfo.NE.0 )
THEN
979 WRITE( nounit, fmt = 9999 )
'DTGEVC(L,S1)', iinfo, n,
987 llwork( j ) = .false.
993 CALL
dtgevc(
'L',
'S', llwork, n, s1, lda, p1, lda,
994 $ evectl( 1, i1+1 ), ldu, dumma, ldu, n, in,
996 IF( iinfo.NE.0 )
THEN
997 WRITE( nounit, fmt = 9999 )
'DTGEVC(L,S2)', iinfo, n,
1003 CALL
dget52( .true., n, s1, lda, p1, lda, evectl, ldu,
1004 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1005 result( 9 ) = dumma( 1 )
1006 IF( dumma( 2 ).GT.thrshn )
THEN
1007 WRITE( nounit, fmt = 9998 )
'Left',
'DTGEVC(HOWMNY=S)',
1008 $ dumma( 2 ), n, jtype, ioldsd
1015 result( 10 ) = ulpinv
1016 CALL
dlacpy(
'F', n, n, q, ldu, evectl, ldu )
1017 CALL
dtgevc(
'L',
'B', llwork, n, s1, lda, p1, lda, evectl,
1018 $ ldu, dumma, ldu, n, in, work, iinfo )
1019 IF( iinfo.NE.0 )
THEN
1020 WRITE( nounit, fmt = 9999 )
'DTGEVC(L,B)', iinfo, n,
1026 CALL
dget52( .true., n, h, lda, t, lda, evectl, ldu, alphr1,
1027 $ alphi1, beta1, work, dumma( 1 ) )
1028 result( 10 ) = dumma( 1 )
1029 IF( dumma( 2 ).GT.thrshn )
THEN
1030 WRITE( nounit, fmt = 9998 )
'Left',
'DTGEVC(HOWMNY=B)',
1031 $ dumma( 2 ), n, jtype, ioldsd
1038 result( 11 ) = ulpinv
1045 llwork( j ) = .true.
1047 DO 170 j = i1 + 1, n
1048 llwork( j ) = .false.
1051 CALL
dtgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, dumma,
1052 $ ldu, evectr, ldu, n, in, work, iinfo )
1053 IF( iinfo.NE.0 )
THEN
1054 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,S1)', iinfo, n,
1062 llwork( j ) = .false.
1064 DO 190 j = i1 + 1, n
1065 llwork( j ) = .true.
1068 CALL
dtgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, dumma,
1069 $ ldu, evectr( 1, i1+1 ), ldu, n, in, work,
1071 IF( iinfo.NE.0 )
THEN
1072 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,S2)', iinfo, n,
1078 CALL
dget52( .false., n, s1, lda, p1, lda, evectr, ldu,
1079 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1080 result( 11 ) = dumma( 1 )
1081 IF( dumma( 2 ).GT.thresh )
THEN
1082 WRITE( nounit, fmt = 9998 )
'Right',
'DTGEVC(HOWMNY=S)',
1083 $ dumma( 2 ), n, jtype, ioldsd
1090 result( 12 ) = ulpinv
1091 CALL
dlacpy(
'F', n, n, z, ldu, evectr, ldu )
1092 CALL
dtgevc(
'R',
'B', llwork, n, s1, lda, p1, lda, dumma,
1093 $ ldu, evectr, ldu, n, in, work, iinfo )
1094 IF( iinfo.NE.0 )
THEN
1095 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,B)', iinfo, n,
1101 CALL
dget52( .false., n, h, lda, t, lda, evectr, ldu,
1102 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1103 result( 12 ) = dumma( 1 )
1104 IF( dumma( 2 ).GT.thresh )
THEN
1105 WRITE( nounit, fmt = 9998 )
'Right',
'DTGEVC(HOWMNY=B)',
1106 $ dumma( 2 ), n, jtype, ioldsd
1115 CALL
dget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1116 $ work, result( 13 ) )
1117 CALL
dget51( 2, n, p1, lda, p2, lda, q, ldu, z, ldu,
1118 $ work, result( 14 ) )
1125 temp1 = max( temp1, abs( alphr1( j )-alphr3( j ) )+
1126 $ abs( alphi1( j )-alphi3( 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 )
'DGG'
1160 WRITE( nounit, fmt = 9996 )
1161 WRITE( nounit, fmt = 9995 )
1162 WRITE( nounit, fmt = 9994 )
'Orthogonal'
1166 WRITE( nounit, fmt = 9993 )
'orthogonal',
'''',
1167 $
'transpose', (
'''', j = 1, 10 )
1171 IF( result( jr ).LT.10000.0d0 )
THEN
1172 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1175 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1186 CALL
dlasum(
'DGG', nounit, nerrs, ntestt )
1189 9999 format(
' DCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1190 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1192 9998 format(
' DCHKGG: ', 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,
' -- Real Generalized eigenvalue problem' )
1199 9996 format(
' Matrix types (see DCHKGG 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, d10.3 )