401 SUBROUTINE ddrges( 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 )
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(
'DDRGES', -info )
540 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
543 safmin =
dlamch(
'Safe minimum' )
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
dgges(
'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 )
'DGGES', 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(
' DDRGES: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
940 $ i6,
', JTYPE=', i6,
', ISEED=(', 4( i4,
',' ), i5,
')' )
942 9998 format(
' DDRGES: DGET53 returned INFO=', i1,
' for eigenvalue ',
943 $ i6,
'.', / 9x,
'N=', i6,
', JTYPE=', i6,
', ISEED=(',
944 $ 4( i4,
',' ), i5,
')' )
946 9997 format(
' DDRGES: 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 DDRGES 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 )