506 SUBROUTINE schkgg( 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
522 LOGICAL DOTYPE( * ), LLWORK( * )
523 INTEGER ISEED( 4 ), NN( * )
524 REAL 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( * ),
538 PARAMETER ( ZERO = 0.0, one = 1.0 )
540 PARAMETER ( MAXTYP = 26 )
544 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
545 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
547 REAL 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 REAL DUMMA( 4 ), RMAGN( 0: 3 )
560 REAL SLAMCH, SLANGE, SLARND
561 EXTERNAL slamch, slange, slarnd
569 INTRINSIC abs, max, min, real, 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(
'SCHKGG', -info )
637 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
640 safmin = slamch(
'Safe minimum' )
641 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
642 safmin = safmin / ulp
643 safmax = one / safmin
644 CALL slabad( safmin, safmax )
658 DO 240 jsize = 1, nsizes
661 rmagn( 2 ) = safmax*ulp / real( 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
719 in = 2*( ( n-1 ) / 2 ) + 1
721 $
CALL slaset(
'Full', n, n, zero, zero, a, lda )
725 CALL slatm4( 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 slaset(
'Full', n, n, zero, zero, b, lda )
743 CALL slatm4( 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 ) = slarnd( 3, iseed )
762 v( jr, jc ) = slarnd( 3, iseed )
764 CALL slarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
766 work( 2*n+jc ) = sign( one, u( jc, jc ) )
768 CALL slarfg( n+1-jc, v( jc, jc ), v( jc+1, jc ), 1,
770 work( 3*n+jc ) = sign( one, v( jc, jc ) )
775 work( 3*n ) = sign( one, slarnd( 2, iseed ) )
778 work( 4*n ) = sign( one, slarnd( 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 sorm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
791 $ lda, work( 2*n+1 ), iinfo )
794 CALL sorm2r(
'R',
'T', n, n, n-1, v, ldu, work( n+1 ),
795 $ a, lda, work( 2*n+1 ), iinfo )
798 CALL sorm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
799 $ lda, work( 2*n+1 ), iinfo )
802 CALL sorm2r(
'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 = slange(
'1', n, n, a, lda, work )
822 bnorm = slange(
'1', n, n, b, lda, work )
826 IF( iinfo.NE.0 )
THEN
827 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
837 CALL slacpy(
' ', n, n, a, lda, h, lda )
838 CALL slacpy(
' ', n, n, b, lda, t, lda )
842 CALL sgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
843 IF( iinfo.NE.0 )
THEN
844 WRITE( nounit, fmt = 9999 )
'SGEQR2', iinfo, n, jtype,
850 CALL sorm2r(
'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 )
'SORM2R', iinfo, n, jtype,
859 CALL slaset(
'Full', n, n, zero, one, u, ldu )
860 CALL sorm2r(
'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 )
'SORM2R', iinfo, n, jtype,
869 CALL sgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
871 IF( iinfo.NE.0 )
THEN
872 WRITE( nounit, fmt = 9999 )
'SGGHRD', iinfo, n, jtype,
881 CALL sget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
883 CALL sget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
885 CALL sget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
887 CALL sget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
896 CALL slacpy(
' ', n, n, h, lda, s2, lda )
897 CALL slacpy(
' ', n, n, t, lda, p2, lda )
901 CALL shgeqz(
'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 )
'SHGEQZ(E)', iinfo, n, jtype,
913 CALL slacpy(
' ', n, n, h, lda, s2, lda )
914 CALL slacpy(
' ', n, n, t, lda, p2, lda )
916 CALL shgeqz(
'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 )
'SHGEQZ(S)', iinfo, n, jtype,
928 CALL slacpy(
' ', n, n, h, lda, s1, lda )
929 CALL slacpy(
' ', n, n, t, lda, p1, lda )
931 CALL shgeqz(
'S',
'I',
'I', n, 1, n, s1, lda, p1, lda,
932 $ alphr1, alphi1, beta1, q, ldu, z, ldu, work,
934 IF( iinfo.NE.0 )
THEN
935 WRITE( nounit, fmt = 9999 )
'SHGEQZ(V)', iinfo, n, jtype,
945 CALL sget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
947 CALL sget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
949 CALL sget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
951 CALL sget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
970 llwork( j ) = .false.
973 CALL stgevc(
'L',
'S', llwork, n, s1, lda, p1, lda, evectl,
974 $ ldu, dumma, ldu, n, in, work, iinfo )
975 IF( iinfo.NE.0 )
THEN
976 WRITE( nounit, fmt = 9999 )
'STGEVC(L,S1)', iinfo, n,
984 llwork( j ) = .false.
990 CALL stgevc(
'L',
'S', llwork, n, s1, lda, p1, lda,
991 $ evectl( 1, i1+1 ), ldu, dumma, ldu, n, in,
993 IF( iinfo.NE.0 )
THEN
994 WRITE( nounit, fmt = 9999 )
'STGEVC(L,S2)', iinfo, n,
1000 CALL sget52( .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',
'STGEVC(HOWMNY=S)',
1005 $ dumma( 2 ), n, jtype, ioldsd
1012 result( 10 ) = ulpinv
1013 CALL slacpy(
'F', n, n, q, ldu, evectl, ldu )
1014 CALL stgevc(
'L',
'B', llwork, n, s1, lda, p1, lda, evectl,
1015 $ ldu, dumma, ldu, n, in, work, iinfo )
1016 IF( iinfo.NE.0 )
THEN
1017 WRITE( nounit, fmt = 9999 )
'STGEVC(L,B)', iinfo, n,
1023 CALL sget52( .true., n, h, lda, t, lda, evectl, ldu, alphr1,
1024 $ alphi1, beta1, work, dumma( 1 ) )
1025 result( 10 ) = dumma( 1 )
1026 IF( dumma( 2 ).GT.thrshn )
THEN
1027 WRITE( nounit, fmt = 9998 )
'Left',
'STGEVC(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 stgevc(
'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 )
'STGEVC(R,S1)', iinfo, n,
1059 llwork( j ) = .false.
1061 DO 190 j = i1 + 1, n
1062 llwork( j ) = .true.
1065 CALL stgevc(
'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( nounit, fmt = 9999 )
'STGEVC(R,S2)', iinfo, n,
1075 CALL sget52( .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',
'STGEVC(HOWMNY=S)',
1080 $ dumma( 2 ), n, jtype, ioldsd
1087 result( 12 ) = ulpinv
1088 CALL slacpy(
'F', n, n, z, ldu, evectr, ldu )
1089 CALL stgevc(
'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 )
'STGEVC(R,B)', iinfo, n,
1098 CALL sget52( .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',
'STGEVC(HOWMNY=B)',
1103 $ dumma( 2 ), n, jtype, ioldsd
1112 CALL sget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1113 $ work, result( 13 ) )
1114 CALL sget51( 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 )
'SGG'
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.0 )
THEN
1169 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1172 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1183 CALL slasum(
'SGG', nounit, nerrs, ntestt )
1186 9999
FORMAT(
' SCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1187 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1189 9998
FORMAT(
' SCHKGG: ', 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 SCHKGG 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, e10.3 )
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
subroutine sgeqr2(M, N, A, LDA, TAU, WORK, INFO)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
subroutine slarfg(N, ALPHA, X, INCX, TAU)
SLARFG generates an elementary reflector (Householder matrix).
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
subroutine schkgg(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)
SCHKGG
subroutine sget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
SGET51
subroutine sget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
SGET52
subroutine slatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
SLATM4
subroutine slasum(TYPE, IOUNIT, IE, NRUN)
SLASUM