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 DOTYPE( * ), LLWORK( * )
515 INTEGER ISEED( 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, * ),
523 $ work( * ), z( 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,
547 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
548 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( ),
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 )
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 ) )
677 ioldsd( j ) = iseed( j )
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 ), 1,
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
839 CALL clacpy(
' ', n, n, b, lda, t, lda )
843 CALL cgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
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.NE.
IF( IINFO0 ) 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.NE.
IF( IINFO0 ) 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.NE.
IF( IINFO0 ) 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.NE.
IF( IINFO0 ) 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.NE.
IF( IINFO0 ) 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.NE.
IF( IINFO0 ) THEN
1052 WRITE( NOUNIT, FMT = 9999 )'ctgevc(r,s1)
', IINFO, N,
1060 LLWORK( J ) = .FALSE.
1062 DO 190 J = I1 + 1, N
1063 LLWORK( J ) = .TRUE.
1066 CALL CTGEVC( 'r
', 's
', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
1067 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
1069.NE.
IF( IINFO0 ) THEN
1070 WRITE( NOUNIT, FMT = 9999 )'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.GT.
IF( DUMMA( 2 )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.NE.
IF( IINFO0 ) 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.GT.
IF( DUMMA( 2 )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.GE.
IF( RESULT( JR )THRESH ) THEN
1152.EQ.
IF( NERRS0 ) 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.LT.
IF( RESULT( JR )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
''', / '',
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 )