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, THRSHN
522 LOGICAL DOTYPE( * ), LLWORK( * )
523 INTEGER ISEED( 4 ), NN( * )
524 DOUBLE PRECISION A( LDA, * ), ALPHI1( * ), ALPHI3( * ),
525 $ ALPHR1( * ), ALPHR3( * ), B( LDA, * ),
526 $ beta1( * ), beta3( * ), evectl( ldu, * ),
527 $ evectr( ldu, * ), h( lda, * ), p1( lda, * ),
528 $ p2( lda, * ), q( ldu, * ), result( 15 ),
529 $ s1( lda, * ), s2( lda, * ), t( lda, * ),
530 $ u( ldu, * ), v( ldu, * ), work( * ),
537 DOUBLE PRECISION ZERO, ONE
538 PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
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, 1, 2*2, 2*1, 2*2, 3, 1, 3,
582 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
584 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
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
657 DO 240 jsize = 1, nsizes
660 rmagn( 2 ) = safmax*ulp / dble( n1 )
661 rmagn( 3 ) = safmin*ulpinv*n1
663 IF( nsizes.NE.1 )
THEN
664 mtypes = min( maxtyp, ntypes )
666 mtypes = min( maxtyp+1, ntypes )
669 DO 230 jtype = 1, mtypes
670 IF( .NOT.dotype( jtype ) )
678 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 dlaset(
'Full', n, n, zero, zero, a, lda )
724 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
725 $ kz2( kazero( jtype ) ), iasign( jtype ),
726 $ rmagn( kamagn( jtype ) ), ulp,
727 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
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 dlaset(
'Full', n, n, zero, zero, b, lda )
742 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
743 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
744 $ rmagn( kbmagn( jtype ) ), one,
745 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
747 iadd = kadd( kbzero( jtype ) )
748 IF( iadd.NE.0 .AND. iadd.LE.n )
749 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
751 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
760 u( jr, jc ) = dlarnd( 3, iseed )
761 v( jr, jc ) = dlarnd( 3, iseed )
763 CALL dlarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
765 work( 2*n+jc ) = sign( one, u( jc, jc ) )
767 CALL dlarfg( n+1-jc, v( jc, jc ), v( jc+1, jc ), 1,
769 work( 3*n+jc ) = sign( one, v( jc, jc ) )
774 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
777 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
783 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
785 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
789 CALL dorm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
790 $ lda, work( 2*n+1 ), iinfo )
793 CALL dorm2r(
'R',
'T', n, n, n-1, v, ldu, work( n+1 ),
794 $ a, lda, work( 2*n+1 ), iinfo )
797 CALL dorm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
798 $ lda, work( 2*n+1 ), iinfo )
801 CALL dorm2r(
'R',
'T', n, n, n-1, v, ldu, work( n+1 ),
802 $ b, lda, work( 2*n+1 ), iinfo )
812 a( jr, jc ) = rmagn( kamagn( jtype ) )*
814 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
820 anorm = dlange(
'1', n, n, a, lda, work )
821 bnorm = dlange(
'1', n, n, b, lda, work )
825 IF( iinfo.NE.0 )
THEN
826 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
836 CALL dlacpy(
' ', n, n, a, lda, h, lda )
837 CALL dlacpy(
' ', n, n, b, lda, t, lda )
841 CALL dgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
842 IF( iinfo.NE.0 )
THEN
843 WRITE( nounit, fmt = 9999 )
'DGEQR2', iinfo, n, jtype,
849 CALL dorm2r(
'L',
'T', n, n, n, t, lda, work, h, lda,
850 $ work( n+1 ), iinfo )
851 IF( iinfo.NE.0 )
THEN
852 WRITE( nounit, fmt = 9999 )
'DORM2R', iinfo, n, jtype,
858 CALL dlaset(
'Full', n, n, zero, one, u, ldu )
859 CALL dorm2r(
'R',
'N', n, n, n, t, lda, work, u, ldu,
860 $ work( n+1 ), iinfo )
861 IF( iinfo.NE.0 )
THEN
862 WRITE( nounit, fmt = 9999 )
'DORM2R', iinfo, n, jtype,
868 CALL dgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
870 IF( iinfo.NE.0 )
THEN
871 WRITE( nounit, fmt = 9999 )
'DGGHRD', iinfo, n, jtype,
880 CALL dget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
882 CALL dget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
884 CALL dget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
886 CALL dget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
895 CALL dlacpy(
' ', n, n, h, lda, s2, lda )
896 CALL dlacpy(
' ', n, n, t, lda, p2, lda )
900 CALL dhgeqz(
'E',
'N',
'N', n, 1, n, s2, lda, p2, lda,
901 $ alphr3, alphi3, beta3, q, ldu, z, ldu, work,
903 IF( iinfo.NE.0 )
THEN
904 WRITE( nounit, fmt = 9999 )
'DHGEQZ(E)', iinfo, n, jtype,
912 CALL dlacpy(
' ', n, n, h, lda, s2, lda )
913 CALL dlacpy(
' ', n, n, t, lda, p2, lda )
915 CALL dhgeqz(
'S',
'N',
'N', n, 1, n, s2, lda, p2, lda,
916 $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
918 IF( iinfo.NE.0 )
THEN
919 WRITE( nounit, fmt = 9999 )
'DHGEQZ(S)', iinfo, n, jtype,
927 CALL dlacpy(
' ', n, n, h, lda, s1, lda )
928 CALL dlacpy(
' ', n, n, t, lda, p1, lda )
930 CALL dhgeqz(
'S',
'I',
'I', n, 1, n, s1, lda, p1, lda,
931 $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
933 IF( iinfo.NE.0 )
THEN
934 WRITE( nounit, fmt = 9999 )
'DHGEQZ(V)', iinfo, n, jtype,
944 CALL dget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
946 CALL dget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
948 CALL dget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
950 CALL dget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
969 llwork( j ) = .false.
972 CALL dtgevc(
'L',
'S', llwork, n, s1, lda, p1, lda, evectl,
973 $ ldu, dumma, ldu, n, in, work, iinfo )
974 IF( iinfo.NE.0 )
THEN
975 WRITE( nounit, fmt = 9999 )
'DTGEVC(L,S1)', iinfo, n,
983 llwork( j ) = .false.
989 CALL dtgevc(
'L',
'S', llwork, n, s1, lda, p1, lda,
990 $ evectl( 1, i1+1 ), ldu, dumma, ldu, n, in,
992 IF( iinfo.NE.0 )
THEN
993 WRITE( nounit, fmt = 9999 )
'DTGEVC(L,S2)', iinfo, n,
999 CALL dget52( .true., n, s1, lda, p1, lda, evectl, ldu,
1000 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1001 result( 9 ) = dumma( 1 )
1002 IF( dumma( 2 ).GT.thrshn )
THEN
1003 WRITE( nounit, fmt = 9998 )
'Left',
'DTGEVC(HOWMNY=S)',
1004 $ dumma( 2 ), n, jtype, ioldsd
1011 result( 10 ) = ulpinv
1012 CALL dlacpy(
'F', n, n, q, ldu, evectl, ldu )
1013 CALL dtgevc(
'L',
'B', llwork, n, s1, lda, p1, lda, evectl,
1014 $ ldu, dumma, ldu, n, in, work, iinfo )
1015 IF( iinfo.NE.0 )
THEN
1016 WRITE( nounit, fmt = 9999 )
'DTGEVC(L,B)', iinfo, n,
1022 CALL dget52( .true., n, h, lda, t, lda, evectl, ldu, alphr1,
1023 $ alphi1, beta1, work, dumma( 1 ) )
1024 result( 10 ) = dumma( 1 )
1025 IF( dumma( 2 ).GT.thrshn )
THEN
1026 WRITE( nounit, fmt = 9998 )
'Left',
'DTGEVC(HOWMNY=B)',
1027 $ dumma( 2 ), n, jtype, ioldsd
1034 result( 11 ) = ulpinv
1041 llwork( j ) = .true.
1043 DO 170 j = i1 + 1, n
1044 llwork( j ) = .false.
1047 CALL dtgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, dumma,
1048 $ ldu, evectr, ldu, n, in, work, iinfo )
1049 IF( iinfo.NE.0 )
THEN
1050 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,S1)', iinfo, n,
1058 llwork( j ) = .false.
1060 DO 190 j = i1 + 1, n
1061 llwork( j ) = .true.
1064 CALL dtgevc(
'R',
'S', llwork, n, s1, lda, p1, lda, dumma,
1065 $ ldu, evectr( 1, i1+1 ), ldu, n, in, work,
1067 IF( iinfo.NE.0 )
THEN
1068 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,S2)', iinfo, n,
1074 CALL dget52( .false., n, s1, lda, p1, lda, evectr, ldu,
1075 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1076 result( 11 ) = dumma( 1 )
1077 IF( dumma( 2 ).GT.thresh )
THEN
1078 WRITE( nounit, fmt = 9998 )
'Right',
'DTGEVC(HOWMNY=S)',
1079 $ dumma( 2 ), n, jtype, ioldsd
1086 result( 12 ) = ulpinv
1087 CALL dlacpy(
'F', n, n, z, ldu, evectr, ldu )
1088 CALL dtgevc(
'R',
'B', llwork, n, s1, lda, p1, lda, dumma,
1089 $ ldu, evectr, ldu, n, in, work, iinfo )
1090 IF( iinfo.NE.0 )
THEN
1091 WRITE( nounit, fmt = 9999 )
'DTGEVC(R,B)', iinfo, n,
1097 CALL dget52( .false., n, h, lda, t, lda, evectr, ldu,
1098 $ alphr1, alphi1, beta1, work, dumma( 1 ) )
1099 result( 12 ) = dumma( 1 )
1100 IF( dumma( 2 ).GT.thresh )
THEN
1101 WRITE( nounit, fmt = 9998 )
'Right',
'DTGEVC(HOWMNY=B)',
1102 $ dumma( 2 ), n, jtype, ioldsd
1111 CALL dget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1112 $ work, result( 13 ) )
1113 CALL dget51( 2, n, p1, lda, p2, lda, q, ldu, z, ldu,
1114 $ work, result( 14 ) )
1121 temp1 = max( temp1, abs( alphr1( j )-alphr3( j ) )+
1122 $ abs( alphi1( j )-alphi3( 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 )
'DGG'
1156 WRITE( nounit, fmt = 9996 )
1157 WRITE( nounit, fmt = 9995 )
1158 WRITE( nounit, fmt = 9994 )
'Orthogonal'
1162 WRITE( nounit, fmt = 9993 )
'orthogonal',
'''',
1163 $
'transpose', (
'''', j = 1, 10 )
1167 IF( result( jr ).LT.10000.0d0 )
THEN
1168 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1171 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1182 CALL dlasum(
'DGG', nounit, nerrs, ntestt )
1185 9999
FORMAT(
' DCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1186 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1188 9998
FORMAT(
' DCHKGG: ', 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,
' -- Real Generalized eigenvalue problem' )
1195 9996
FORMAT(
' Matrix types (see DCHKGG 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, d10.3 )
subroutine xerbla(srname, info)
subroutine dchkgg(nsizes, nn, ntypes, dotype, iseed, thresh, tstdif, thrshn, nounit, a, lda, b, h, t, s1, s2, p1, p2, u, ldu, v, q, z, alphr1, alphi1, beta1, alphr3, alphi3, beta3, evectl, evectr, work, lwork, llwork, result, info)
DCHKGG
subroutine dget51(itype, n, a, lda, b, ldb, u, ldu, v, ldv, work, result)
DGET51
subroutine dget52(left, n, a, lda, b, ldb, e, lde, alphar, alphai, beta, work, result)
DGET52
subroutine dlasum(type, iounit, ie, nrun)
DLASUM
subroutine dlatm4(itype, n, nz1, nz2, isign, amagn, rcond, triang, idist, iseed, a, lda)
DLATM4
subroutine dgeqr2(m, n, a, lda, tau, work, info)
DGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
subroutine dgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
DGGHRD
subroutine dhgeqz(job, compq, compz, n, ilo, ihi, h, ldh, t, ldt, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, info)
DHGEQZ
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
subroutine dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine dtgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
DTGEVC
subroutine dorm2r(side, trans, m, n, k, a, lda, tau, c, ldc, work, info)
DORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...