506 SUBROUTINE dchkgg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
507 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
508 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHR1, ALPHI1,
509 $ BETA1, ALPHR3, ALPHI3, BETA3, EVECTL, EVECTR,
510 $ WORK, LWORK, LLWORK, RESULT, INFO )
518 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
519 DOUBLE PRECISION THRESH
522 LOGICAL DOTYPE( * ), LLWORK( * )
523 INTEGER ISEED( 4 ), NN( * )
524 DOUBLE PRECISION A( LDA, * ), ( * ), ALPHI3( * ),
525 $ ALPHR1( * ), ALPHR3( * ), B( LDA, * ),
526 $ beta1( * ), beta3( * ), evectl( ldu, * ),
527 $ evectr( ldu, * ), h( lda, * ), p1( lda, * )
528 $ p2( lda, * ), q( ldu, * ), result(
530 $ u( ldu, * ), v( ldu, * ), work( * ),
537 DOUBLE PRECISION , ONE
540 PARAMETER ( MAXTYP = 26 )
544 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
545 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
547 DOUBLE PRECISION ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2,
551 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
552 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
553 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
554 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
555 $ kbzero( maxtyp ), kclass( maxtyp ),
556 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
557 DOUBLE PRECISION DUMMA( 4 ), RMAGN( 0: 3 )
560 DOUBLE PRECISION DLAMCH, DLANGE, DLARND
561 EXTERNAL dlamch, dlange, dlarnd
569 INTRINSIC abs, dble,
max,
min, sign
572 DATA kclass / 15*1, 10*2, 1*3 /
573 DATA kz1 / 0, 1, 2, 1, 3, 3 /
574 DATA kz2 / 0, 0, 1, 2, 1, 1 /
575 DATA kadd / 0, 0, 0, 0, 3, 2 /
576 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
577 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
578 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
579 $ 1, 1, -4, 2, -4, 8*8, 0 /
580 DATA kazero / 6*1, 2,
582 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
586 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
588 DATA ktrian / 16*0, 10*1 /
589 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
591 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
602 nmax =
max( nmax, nn( j ) )
610 lwkopt =
max( 6*nmax, 2*nmax*nmax, 1 )
614 IF( nsizes.LT.0 )
THEN
616 ELSE IF( badnn )
THEN
618 ELSE IF( ntypes.LT.0 )
THEN
620 ELSE IF( thresh.LT.zero )
THEN
622 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
624 ELSE IF( ldu.LE.1 .OR. ldu.LT.nmax )
THEN
626 ELSE IF( lwkopt.GT.lwork )
THEN
631 CALL xerbla(
'DCHKGG', -info )
637 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
640 safmin = dlamch(
'Safe minimum' )
641 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
642 safmin = safmin / ulp
643 safmax = one / safmin
644 CALL dlabad( safmin, safmax )
658 DO 240 jsize = 1, nsizes
661 rmagn( 2 ) = safmax*ulp / dble( n1 )
662 rmagn( 3 ) = safmin*ulpinv*n1
664 IF( nsizes.NE.1 )
THEN
665 mtypes =
min( maxtyp, ntypes )
667 mtypes =
min( maxtyp+1, ntypes )
670 DO 230 jtype = 1, mtypes
671 IF( .NOT.dotype( jtype ) )
679 ioldsd( j ) = iseed( j )
711 IF( mtypes.GT.maxtyp )
714 IF( kclass( jtype ).LT.3 )
THEN
718 IF( abs( katype( jtype ) ).EQ.3 )
THEN
721 $
CALL dlaset(
'Full', n, n, zero, zero, a, lda )
725 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
726 $ kz2( kazero( jtype ) ), iasign( jtype ),
727 $ rmagn( kamagn( jtype ) ), ulp,
728 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
730 iadd = kadd( kazero( jtype ) )
731 IF( iadd.GT.0 .AND. iadd.LE.n )
732 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
736 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
737 in = 2*( ( n-1 ) / 2 ) + 1
739 $
CALL dlaset(
'Full', n, n, zero, zero, b, lda )
743 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
744 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
745 $ rmagn( kbmagn( jtype ) ), one,
746 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
748 iadd = kadd( kbzero( jtype ) )
749 IF( iadd.NE.0 .AND. iadd.LE.n )
750 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
752 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
761 u( jr, jc ) = dlarnd( 3, iseed )
762 v( jr, jc ) = dlarnd( 3, iseed )
764 CALL dlarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
766 work( 2*n+jc ) = sign( one, u( jc, jc ) )
768 CALL dlarfg( n+1-jc, v( jc, jc ), v( jc+1, jc
770 work( 3*n+jc ) = sign( one, v( jc, jc ) )
775 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
778 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
784 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
786 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
790 CALL dorm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
791 $ lda, work( 2*n+1 ), iinfo )
794 CALL dorm2r(
'R',
'T', n, n, n-1, v, ldu, work( n+1 ),
795 $ a, lda, work( 2*n+1 ), iinfo )
798 CALL dorm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
799 $ lda, work( 2*n+1 ), iinfo )
802 CALL dorm2r(
'R',
'T', 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 = dlange(
'1', n, n, a, lda, work )
822 bnorm = dlange(
'1', n, n, b, lda, work )
826 IF( iinfo.NE.0 )
THEN
827 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
837 CALL dlacpy(
' ', n, n, a, lda, h, lda )
838 CALL dlacpy(
' ', n, n, b, lda, t, lda )
843 IF( iinfo.NE.0 )
THEN
844 WRITE( nounit, fmt = 9999 )
'DGEQR2', iinfo, n, jtype,
850 CALL dorm2r(
'L',
'T', n, n, n, t, lda, work, h, lda,
851 $ work( n+1 ), iinfo )
852 IF( iinfo.NE.0 )
THEN
853 WRITE( nounit, fmt = 9999 )
'DORM2R', iinfo, n, jtype,
859 CALL dlaset(
'Full', n, n, zero, one, u, ldu )
860 CALL dorm2r(
'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 )
'DORM2R', iinfo, n, jtype,
869 CALL dgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
871 IF( iinfo.NE.0 )
THEN
872 WRITE( nounit, fmt = 9999 )
'DGGHRD', iinfo, n, jtype,
881 CALL dget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
883 CALL dget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
885 CALL dget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
887 CALL dget51( 3, n, b, lda, t
896 CALL dlacpy(
' ', n, n, h, lda, s2, lda )
897 CALL dlacpy(
' ', n, n, t, lda, p2, lda )
901 CALL dhgeqz(
'E',
'N',
'N', n, 1, n, s2, lda, p2, lda,
902 $ alphr3, alphi3, beta3, q, ldu, z, ldu, work,
904 IF( iinfo.NE.0 )
THEN
905 WRITE( nounit, fmt = 9999 )
'DHGEQZ(E)', iinfo, n, jtype,
913 CALL dlacpy(
' ', n, n, h, lda, s2, lda )
914 CALL dlacpy(
' ', n, n, t, lda, p2, lda )
916 CALL dhgeqz(
'S',
'N',
'N', n, 1, n, s2, lda, p2, lda,
917 $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
919 IF( iinfo.NE.0 )
THEN
920 WRITE( nounit, fmt = 9999 )
'DHGEQZ(S)', iinfo, n, jtype,
928 CALL dlacpy(
' ', n, n, h, lda, s1, lda )
929 CALL dlacpy(
' ', n, n, t, lda, p1, lda )
931 CALL dhgeqz( 's
', 'i
', 'i
', N, 1, N, S1, LDA, P1, LDA,
932 $ ALPHR1, ALPHI1, BETA1, Q, LDU, Z, LDU, WORK,
934.NE.
IF( IINFO0 ) THEN
935 WRITE( NOUNIT, FMT = 9999 )'dhgeqz(v)
', IINFO, N, JTYPE,
945 CALL DGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
947 CALL DGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
949 CALL DGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
951 CALL DGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
970 LLWORK( J ) = .FALSE.
973 CALL DTGEVC( 'l
', 's
', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
974 $ LDU, DUMMA, LDU, N, IN, WORK, IINFO )
975.NE.
IF( IINFO0 ) THEN
976 WRITE( NOUNIT, FMT = 9999 )'dtgevc(l,s1)
', IINFO, N,
984 LLWORK( J ) = .FALSE.
990 CALL DTGEVC( 'l
', 's
', LLWORK, N, S1, LDA, P1, LDA,
991 $ EVECTL( 1, I1+1 ), LDU, DUMMA, LDU, N, IN,
993.NE.
IF( IINFO0 ) THEN
994 WRITE( NOUNIT, FMT = 9999 )'dtgevc(l,s2)', iinfo, n,
1000 CALL dget52( .true., n, s1, lda, p1, lda, evectl, ldu,
1001 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1002 result( 9 ) = dumma( 1 )
1003 IF( dumma( 2 ).GT.thrshn )
THEN
1004 WRITE( nounit, fmt = 9998 )
'Left',
'DTGEVC(HOWMNY=S)',
1005 $ dumma( 2 ), n, jtype, ioldsd
1012 result( 10 ) = ulpinv
1013 CALL dlacpy( 'f
', N, N, Q, LDU, EVECTL, LDU )
1014 CALL DTGEVC( 'l
', 'b
', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
1015 $ LDU, DUMMA, LDU, N, IN, WORK, IINFO )
1016.NE.
IF( IINFO0 ) THEN
1017 WRITE( NOUNIT, FMT = 9999 )'dtgevc(l,b)
', IINFO, N,
1023 CALL DGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHR1,
1024 $ ALPHI1, BETA1, WORK, DUMMA( 1 ) )
1025 RESULT( 10 ) = DUMMA( 1 )
1026.GT.
IF( DUMMA( 2 )THRSHN ) THEN
1027 WRITE( NOUNIT, FMT = 9998 )'left
', 'dtgevc(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 dtgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, dumma,
1049 $ ldu, evectr, ldu, n, in, work, iinfo )
1050 IF( iinfo.NE.0 )
THEN
1051 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,S1)', iinfo, n,
1059 llwork( j ) = .false.
1061 DO 190 j = i1 + 1, n
1062 llwork( j ) = .true.
1065 CALL dtgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, dumma,
1066 $ ldu, evectr( 1, i1+1 ), ldu, n, in, work,
1068 IF( iinfo.NE.0 )
THEN
1069 WRITE'DTGEVC(R,S2)', iinfo, n,
1075 CALL dget52( .false., n, s1, lda, p1, lda, evectr, ldu,
1076 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1077 result( 11 ) = dumma( 1 )
1078 IF( dumma( 2 ).GT.thresh )
THEN
1079 WRITE( nounit, fmt = 9998 )
'Right',
'DTGEVC(HOWMNY=S)',
1080 $ dumma( 2 ), n, jtype, ioldsd
1087 result( 12 ) = ulpinv
1088 CALL dlacpy(
'F', n, n, z, ldu, evectr, ldu )
1089 CALL dtgevc(
'R',
'B', llwork, n, s1, lda, p1, lda, dumma,
1090 $ ldu, evectr, ldu, n, in, work, iinfo )
1091 IF( iinfo.NE.0 )
THEN
1092 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,B)', iinfo, n,
1098 CALL dget52( .false., n, h, lda, t, lda, evectr, ldu,
1099 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1100 result( 12 ) = dumma( 1 )
1101 IF( dumma( 2 ).GT.thresh )
THEN
1102 WRITE( nounit, fmt = 9998 )
'Right',
'DTGEVC(HOWMNY=B)',
1103 $ dumma( 2 ), n, jtype, ioldsd
1112 CALL dget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1113 $ work, result( 13 ) )
1114 CALL dget51( 2, n, p1, lda, p2, lda, q, ldu, z, ldu,
1115 $ work, result( 14 ) )
1122 temp1 =
max( temp1, abs( alphr1( j )-alphr3( j ) )+
1123 $ abs( alphi1( j )-alphi3( 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 )
'DGG'
1157 WRITE( nounit, fmt = 9996 )
1158 WRITE( nounit, fmt = 9995 )
1159 WRITE( nounit, fmt = 9994 )
'Orthogonal'
1163 WRITE( nounit, fmt = 9993 )
'orthogonal',
'''',
1164 $
'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(
'DGG', nounit, nerrs, ntestt )
1186 9999
FORMAT(
' DCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1187 $ i6,
', JTYPE=', i6, ', iseed=(
', 3( I5, ',
' ), I5, ')
' )
1189 9998 FORMAT( ' dchkgg:
', 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, ' -- real generalized eigenvalue problem
' )
1196 9996 FORMAT( ' matrix types(see
dchkgg 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 )