401 SUBROUTINE ddrges3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
402 $ nounit, a, lda, b, s, t, q, ldq, z, alphar,
403 $ alphai, beta, work, lwork, result, bwork,
412 INTEGER INFO, LDA, LDQ, LWORK, NOUNIT, NSIZES, NTYPES
413 DOUBLE PRECISION THRESH
416 LOGICAL BWORK( * ), DOTYPE( * )
417 INTEGER ISEED( 4 ), NN( * )
418 DOUBLE PRECISION A( lda, * ), ALPHAI( * ), ALPHAR( * ),
419 $ b( lda, * ), beta( * ), q( ldq, * ),
420 $ result( 13 ), s( lda, * ), t( lda, * ),
421 $ work( * ), z( ldq, * )
427 DOUBLE PRECISION ZERO, ONE
428 parameter ( zero = 0.0d+0, one = 1.0d+0 )
430 parameter ( maxtyp = 26 )
433 LOGICAL BADNN, ILABAD
435 INTEGER I, I1, IADD, IERR, IINFO, IN, ISORT, J, JC, JR,
436 $ jsize, jtype, knteig, maxwrk, minwrk, mtypes,
437 $ n, n1, nb, nerrs, nmats, nmax, ntest, ntestt,
439 DOUBLE PRECISION SAFMAX, SAFMIN, TEMP1, TEMP2, ULP, ULPINV
442 INTEGER IASIGN( maxtyp ), IBSIGN( maxtyp ),
443 $ ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
444 $ katype( maxtyp ), kazero( maxtyp ),
445 $ kbmagn( maxtyp ), kbtype( maxtyp ),
446 $ kbzero( maxtyp ), kclass( maxtyp ),
447 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
448 DOUBLE PRECISION RMAGN( 0: 3 )
453 DOUBLE PRECISION DLAMCH, DLARND
454 EXTERNAL dlctes, ilaenv, dlamch, dlarnd
461 INTRINSIC abs, dble, max, min, sign
464 DATA kclass / 15*1, 10*2, 1*3 /
465 DATA kz1 / 0, 1, 2, 1, 3, 3 /
466 DATA kz2 / 0, 0, 1, 2, 1, 1 /
467 DATA kadd / 0, 0, 0, 0, 3, 2 /
468 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
469 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
470 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
471 $ 1, 1, -4, 2, -4, 8*8, 0 /
472 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
474 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
476 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
478 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
480 DATA ktrian / 16*0, 10*1 /
481 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
483 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
494 nmax = max( nmax, nn( j ) )
499 IF( nsizes.LT.0 )
THEN
501 ELSE IF( badnn )
THEN
503 ELSE IF( ntypes.LT.0 )
THEN
505 ELSE IF( thresh.LT.zero )
THEN
507 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
509 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax )
THEN
521 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
522 minwrk = max( 10*( nmax+1 ), 3*nmax*nmax )
523 nb = max( 1, ilaenv( 1,
'DGEQRF',
' ', nmax, nmax, -1, -1 ),
524 $ ilaenv( 1,
'DORMQR',
'LT', nmax, nmax, nmax, -1 ),
525 $ ilaenv( 1,
'DORGQR',
' ', nmax, nmax, nmax, -1 ) )
526 maxwrk = max( 10*( nmax+1 ), 2*nmax+nmax*nb, 3*nmax*nmax )
530 IF( lwork.LT.minwrk )
534 CALL xerbla(
'DDRGES3', -info )
540 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
543 safmin = dlamch(
'Safe minimum' )
544 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
545 safmin = safmin / ulp
546 safmax = one / safmin
547 CALL dlabad( safmin, safmax )
561 DO 190 jsize = 1, nsizes
564 rmagn( 2 ) = safmax*ulp / dble( n1 )
565 rmagn( 3 ) = safmin*ulpinv*dble( n1 )
567 IF( nsizes.NE.1 )
THEN
568 mtypes = min( maxtyp, ntypes )
570 mtypes = min( maxtyp+1, ntypes )
575 DO 180 jtype = 1, mtypes
576 IF( .NOT.dotype( jtype ) )
584 ioldsd( j ) = iseed( j )
616 IF( mtypes.GT.maxtyp )
619 IF( kclass( jtype ).LT.3 )
THEN
623 IF( abs( katype( jtype ) ).EQ.3 )
THEN
624 in = 2*( ( n-1 ) / 2 ) + 1
626 $
CALL dlaset(
'Full', n, n, zero, zero, a, lda )
630 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
631 $ kz2( kazero( jtype ) ), iasign( jtype ),
632 $ rmagn( kamagn( jtype ) ), ulp,
633 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
635 iadd = kadd( kazero( jtype ) )
636 IF( iadd.GT.0 .AND. iadd.LE.n )
637 $ a( iadd, iadd ) = one
641 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
642 in = 2*( ( n-1 ) / 2 ) + 1
644 $
CALL dlaset(
'Full', n, n, zero, zero, b, lda )
648 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
649 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
650 $ rmagn( kbmagn( jtype ) ), one,
651 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
653 iadd = kadd( kbzero( jtype ) )
654 IF( iadd.NE.0 .AND. iadd.LE.n )
655 $ b( iadd, iadd ) = one
657 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
666 q( jr, jc ) = dlarnd( 3, iseed )
667 z( jr, jc ) = dlarnd( 3, iseed )
669 CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
671 work( 2*n+jc ) = sign( one, q( jc, jc ) )
673 CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
675 work( 3*n+jc ) = sign( one, z( jc, jc ) )
680 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
683 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
689 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
691 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
695 CALL dorm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
696 $ lda, work( 2*n+1 ), iinfo )
699 CALL dorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
700 $ a, lda, work( 2*n+1 ), iinfo )
703 CALL dorm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
704 $ lda, work( 2*n+1 ), iinfo )
707 CALL dorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
708 $ b, lda, work( 2*n+1 ), iinfo )
718 a( jr, jc ) = rmagn( kamagn( jtype ) )*
720 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
728 IF( iinfo.NE.0 )
THEN
729 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
744 IF( isort.EQ.0 )
THEN
754 CALL dlacpy(
'Full', n, n, a, lda, s, lda )
755 CALL dlacpy(
'Full', n, n, b, lda, t, lda )
756 ntest = 1 + rsub + isort
757 result( 1+rsub+isort ) = ulpinv
758 CALL dgges3(
'V',
'V', sort, dlctes, n, s, lda, t, lda,
759 $ sdim, alphar, alphai, beta, q, ldq, z, ldq,
760 $ work, lwork, bwork, iinfo )
761 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
762 result( 1+rsub+isort ) = ulpinv
763 WRITE( nounit, fmt = 9999 )
'DGGES3', iinfo, n, jtype,
773 IF( isort.EQ.0 )
THEN
774 CALL dget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
775 $ work, result( 1 ) )
776 CALL dget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
777 $ work, result( 2 ) )
779 CALL dget54( n, a, lda, b, lda, s, lda, t, lda, q,
780 $ ldq, z, ldq, work, result( 7 ) )
782 CALL dget51( 3, n, a, lda, t, lda, q, ldq, q, ldq, work,
784 CALL dget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
796 IF( alphai( j ).EQ.zero )
THEN
797 temp2 = ( abs( alphar( j )-s( j, j ) ) /
798 $ max( safmin, abs( alphar( j ) ), abs( s( j,
799 $ j ) ) )+abs( beta( j )-t( j, j ) ) /
800 $ max( safmin, abs( beta( j ) ), abs( t( j,
804 IF( s( j+1, j ).NE.zero )
THEN
806 result( 5+rsub ) = ulpinv
810 IF( s( j, j-1 ).NE.zero )
THEN
812 result( 5+rsub ) = ulpinv
817 IF( alphai( j ).GT.zero )
THEN
822 IF( i1.LE.0 .OR. i1.GE.n )
THEN
824 ELSE IF( i1.LT.n-1 )
THEN
825 IF( s( i1+2, i1+1 ).NE.zero )
THEN
827 result( 5+rsub ) = ulpinv
829 ELSE IF( i1.GT.1 )
THEN
830 IF( s( i1, i1-1 ).NE.zero )
THEN
832 result( 5+rsub ) = ulpinv
835 IF( .NOT.ilabad )
THEN
836 CALL dget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
837 $ beta( j ), alphar( j ),
838 $ alphai( j ), temp2, ierr )
840 WRITE( nounit, fmt = 9998 )ierr, j, n,
849 temp1 = max( temp1, temp2 )
851 WRITE( nounit, fmt = 9997 )j, n, jtype, ioldsd
854 result( 6+rsub ) = temp1
856 IF( isort.GE.1 )
THEN
864 IF( dlctes( alphar( i ), alphai( i ),
865 $ beta( i ) ) .OR. dlctes( alphar( i ),
866 $ -alphai( i ), beta( i ) ) )
THEN
870 IF( ( dlctes( alphar( i+1 ), alphai( i+1 ),
871 $ beta( i+1 ) ) .OR. dlctes( alphar( i+1 ),
872 $ -alphai( i+1 ), beta( i+1 ) ) ) .AND.
873 $ ( .NOT.( dlctes( alphar( i ), alphai( i ),
874 $ beta( i ) ) .OR. dlctes( alphar( i ),
875 $ -alphai( i ), beta( i ) ) ) ) .AND.
876 $ iinfo.NE.n+2 )
THEN
877 result( 12 ) = ulpinv
881 IF( sdim.NE.knteig )
THEN
882 result( 12 ) = ulpinv
892 ntestt = ntestt + ntest
897 IF( result( jr ).GE.thresh )
THEN
902 IF( nerrs.EQ.0 )
THEN
903 WRITE( nounit, fmt = 9996 )
'DGS'
907 WRITE( nounit, fmt = 9995 )
908 WRITE( nounit, fmt = 9994 )
909 WRITE( nounit, fmt = 9993 )
'Orthogonal'
913 WRITE( nounit, fmt = 9992 )
'orthogonal',
'''',
914 $
'transpose', (
'''', j = 1, 8 )
918 IF( result( jr ).LT.10000.0d0 )
THEN
919 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
922 WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
933 CALL alasvm(
'DGS', nounit, nerrs, ntestt, 0 )
939 9999
FORMAT(
' DDRGES3: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
940 $ i6,
', JTYPE=', i6,
', ISEED=(', 4( i4,
',' ), i5,
')' )
942 9998
FORMAT(
' DDRGES3: DGET53 returned INFO=', i1,
' for eigenvalue ',
943 $ i6,
'.', / 9x,
'N=', i6,
', JTYPE=', i6,
', ISEED=(',
944 $ 4( i4,
',' ), i5,
')' )
946 9997
FORMAT(
' DDRGES3: S not in Schur form at eigenvalue ', i6,
'.',
947 $ / 9x,
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
950 9996
FORMAT( / 1x, a3,
' -- Real Generalized Schur form driver' )
952 9995
FORMAT(
' Matrix types (see DDRGES3 for details): ' )
954 9994
FORMAT(
' Special Matrices:', 23x,
955 $
'(J''=transposed Jordan block)',
956 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
957 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
958 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
959 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
960 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
961 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
962 9993
FORMAT(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
963 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
964 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
965 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
966 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
967 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
968 $
'23=(small,large) 24=(small,small) 25=(large,large)',
969 $ /
' 26=random O(1) matrices.' )
971 9992
FORMAT( /
' Tests performed: (S is Schur, T is triangular, ',
972 $
'Q and Z are ', a,
',', / 19x,
973 $
'l and r are the appropriate left and right', / 19x,
974 $
'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
975 $
' means ', a,
'.)', /
' Without ordering: ',
976 $ /
' 1 = | A - Q S Z', a,
977 $
' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
978 $
' | / ( |B| n ulp )', /
' 3 = | I - QQ', a,
979 $
' | / ( n ulp ) 4 = | I - ZZ', a,
980 $
' | / ( n ulp )', /
' 5 = A is in Schur form S',
981 $ /
' 6 = difference between (alpha,beta)',
982 $
' and diagonals of (S,T)', /
' With ordering: ',
983 $ /
' 7 = | (A,B) - Q (S,T) Z', a,
984 $
' | / ( |(A,B)| n ulp ) ', /
' 8 = | I - QQ', a,
985 $
' | / ( n ulp ) 9 = | I - ZZ', a,
986 $
' | / ( n ulp )', /
' 10 = A is in Schur form S',
987 $ /
' 11 = difference between (alpha,beta) and diagonals',
988 $
' of (S,T)', /
' 12 = SDIM is the correct number of ',
989 $
'selected eigenvalues', / )
990 9991
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
991 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
992 9990
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
993 $ 4( i4,
',' ),
' result ', i2,
' is', 1p, d10.3 )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine ddrges3(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, ALPHAR, ALPHAI, BETA, WORK, LWORK, RESULT, BWORK, INFO)
DDRGES3
subroutine dgges3(JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, BWORK, INFO)
DGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
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 dget54(N, A, LDA, B, LDB, S, LDS, T, LDT, U, LDU, V, LDV, WORK, RESULT)
DGET54
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...
subroutine dlatm4(ITYPE, N, NZ1, NZ2, ISIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
DLATM4
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
DGET51
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
DGET53
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).