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
657 DO 240 jsize = 1, nsizes
660 rmagn( 2 ) = safmax*ulp / real( 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 slaset(
'Full', n, n, zero, zero, a, lda )
724 CALL slatm4( 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 slaset(
'Full', n, n, zero, zero, b, lda )
742 CALL slatm4( 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 ) = slarnd( 3, iseed )
761 v( jr, jc ) = slarnd( 3, iseed )
763 CALL slarfg( n+1-jc, u( jc, jc ), u( jc+1, jc ), 1,
765 work( 2*n+jc ) = sign( one, u( jc, jc ) )
767 CALL slarfg( 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, slarnd( 2, iseed ) )
777 work( 4*n ) = sign( one, slarnd( 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 sorm2r(
'L',
'N', n, n, n-1, u, ldu, work, a,
790 $ lda, work( 2*n+1 ), iinfo )
793 CALL sorm2r(
'R',
'T', n, n, n-1, v, ldu, work( n+1 ),
794 $ a, lda, work( 2*n+1 ), iinfo )
797 CALL sorm2r(
'L',
'N', n, n, n-1, u, ldu, work, b,
798 $ lda, work( 2*n+1 ), iinfo )
801 CALL sorm2r(
'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 = slange(
'1', n, n, a, lda, work )
821 bnorm = slange(
'1', n, n, b, lda, work )
825 IF( iinfo.NE.0 )
THEN
826 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
836 CALL slacpy(
' ', n, n, a, lda, h, lda )
837 CALL slacpy(
' ', n, n, b, lda, t, lda )
841 CALL sgeqr2( n, n, t, lda, work, work( n+1 ), iinfo )
842 IF( iinfo.NE.0 )
THEN
843 WRITE( nounit, fmt = 9999 )
'SGEQR2', iinfo, n, jtype,
849 CALL sorm2r(
'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 )
'SORM2R', iinfo, n, jtype,
858 CALL slaset(
'Full', n, n, zero, one, u, ldu )
859 CALL sorm2r(
'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 )
'SORM2R', iinfo, n, jtype,
868 CALL sgghrd(
'V',
'I', n, 1, n, h, lda, t, lda, u, ldu, v,
870 IF( iinfo.NE.0 )
THEN
871 WRITE( nounit, fmt = 9999 )
'SGGHRD', iinfo, n, jtype,
880 CALL sget51( 1, n, a, lda, h, lda, u, ldu, v, ldu, work,
882 CALL sget51( 1, n, b, lda, t, lda, u, ldu, v, ldu, work,
884 CALL sget51( 3, n, b, lda, t, lda, u, ldu, u, ldu, work,
886 CALL sget51( 3, n, b, lda, t, lda, v, ldu, v, ldu, work,
895 CALL slacpy(
' ', n, n, h, lda, s2, lda )
896 CALL slacpy(
' ', n, n, t, lda, p2, lda )
900 CALL shgeqz(
'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 )
'SHGEQZ(E)', iinfo, n, jtype,
912 CALL slacpy(
' ', n, n, h, lda, s2, lda )
913 CALL slacpy(
' ', n, n, t, lda, p2, lda )
915 CALL shgeqz(
'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 )
'SHGEQZ(S)', iinfo, n, jtype,
927 CALL slacpy(
' ', n, n, h, lda, s1, lda )
928 CALL slacpy(
' ', n, n, t, lda, p1, lda )
930 CALL shgeqz(
'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 )
'SHGEQZ(V)', iinfo, n, jtype,
944 CALL sget51( 1, n, h, lda, s1, lda, q, ldu, z, ldu, work,
946 CALL sget51( 1, n, t, lda, p1, lda, q, ldu, z, ldu, work,
948 CALL sget51( 3, n, t, lda, p1, lda, q, ldu, q, ldu, work,
950 CALL sget51( 3, n, t, lda, p1, lda, z, ldu, z, ldu, work,
969 llwork( j ) = .false.
972 CALL stgevc(
'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 )
'STGEVC(L,S1)', iinfo, n,
983 llwork( j ) = .false.
989 CALL stgevc(
'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 )
'STGEVC(L,S2)', iinfo, n,
999 CALL sget52( .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',
'STGEVC(HOWMNY=S)',
1004 $ dumma( 2 ), n, jtype, ioldsd
1011 result( 10 ) = ulpinv
1012 CALL slacpy(
'F', n, n, q, ldu, evectl, ldu )
1013 CALL stgevc(
'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 )
'STGEVC(L,B)', iinfo, n,
1022 CALL sget52( .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',
'STGEVC(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 stgevc(
'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 )
'STGEVC(R,S1)', iinfo, n,
1058 llwork( j ) = .false.
1060 DO 190 j = i1 + 1, n
1061 llwork( j ) = .true.
1064 CALL stgevc(
'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 )
'STGEVC(R,S2)', iinfo, n,
1074 CALL sget52( .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',
'STGEVC(HOWMNY=S)',
1079 $ dumma( 2 ), n, jtype, ioldsd
1086 result( 12 ) = ulpinv
1087 CALL slacpy(
'F', n, n, z, ldu, evectr, ldu )
1088 CALL stgevc(
'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 )
'STGEVC(R,B)', iinfo, n,
1097 CALL sget52( .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',
'STGEVC(HOWMNY=B)',
1102 $ dumma( 2 ), n, jtype, ioldsd
1111 CALL sget51( 2, n, s1, lda, s2, lda, q, ldu, z, ldu,
1112 $ work, result( 13 ) )
1113 CALL sget51( 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 )
'SGG'
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.0 )
THEN
1168 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
1171 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
1182 CALL slasum(
'SGG', nounit, nerrs, ntestt )
1185 9999
FORMAT(
' SCHKGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
1186 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
1188 9998
FORMAT(
' SCHKGG: ', 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 SCHKGG 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, e10.3 )
subroutine xerbla(srname, info)
subroutine sgeqr2(m, n, a, lda, tau, work, info)
SGEQR2 computes the QR factorization of a general rectangular matrix using an unblocked algorithm.
subroutine sgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
SGGHRD
subroutine shgeqz(job, compq, compz, n, ilo, ihi, h, ldh, t, ldt, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, info)
SHGEQZ
subroutine slacpy(uplo, m, n, a, lda, b, ldb)
SLACPY copies all or part of one two-dimensional array to another.
subroutine slarfg(n, alpha, x, incx, tau)
SLARFG generates an elementary reflector (Householder matrix).
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 stgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
STGEVC
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 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 slasum(type, iounit, ie, nrun)
SLASUM
subroutine slatm4(itype, n, nz1, nz2, isign, amagn, rcond, triang, idist, iseed, a, lda)
SLATM4