450 SUBROUTINE ddrvgg( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
451 $ thrshn, nounit, a, lda, b, s, t, s2, t2, q,
452 $ ldq, z, alphr1, alphi1, beta1, alphr2, alphi2,
453 $ beta2, vl, vr, work, lwork, result, info )
461 INTEGER info, lda, ldq, lwork, nounit, nsizes, ntypes
462 DOUBLE PRECISION thresh, thrshn
466 INTEGER iseed( 4 ), nn( * )
467 DOUBLE PRECISION a( lda, * ), alphi1( * ), alphi2( * ),
468 $ alphr1( * ), alphr2( * ), b( lda, * ),
469 $ beta1( * ), beta2( * ), q( ldq, * ),
470 $ result( * ), s( lda, * ), s2( lda, * ),
471 $ t( lda, * ), t2( lda, * ), vl( ldq, * ),
472 $ vr( ldq, * ), work( * ), z( ldq, * )
478 DOUBLE PRECISION zero, one
479 parameter( zero = 0.0d0, one = 1.0d0 )
481 parameter( maxtyp = 26 )
484 LOGICAL badnn, ilabad
485 INTEGER i1, iadd, iinfo, in, j, jc, jr, jsize, jtype,
486 $ lwkopt, mtypes, n, n1, nb, nbz, nerrs, nmats,
487 $ nmax, ns, ntest, ntestt
488 DOUBLE PRECISION safmax, safmin, temp1, temp2, ulp, ulpinv
491 INTEGER iasign( maxtyp ), ibsign( maxtyp ),
492 $ ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
493 $ katype( maxtyp ), kazero( maxtyp ),
494 $ kbmagn( maxtyp ), kbtype( maxtyp ),
495 $ kbzero( maxtyp ), kclass( maxtyp ),
496 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
497 DOUBLE PRECISION dumma( 4 ), rmagn( 0: 3 )
510 INTRINSIC abs, dble, max, min, sign
513 DATA kclass / 15*1, 10*2, 1*3 /
514 DATA kz1 / 0, 1, 2, 1, 3, 3 /
515 DATA kz2 / 0, 0, 1, 2, 1, 1 /
516 DATA kadd / 0, 0, 0, 0, 3, 2 /
517 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
518 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
519 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
520 $ 1, 1, -4, 2, -4, 8*8, 0 /
521 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
523 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
525 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
527 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
529 DATA ktrian / 16*0, 10*1 /
530 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
532 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
543 nmax = max( nmax, nn( j ) )
551 nb = max( 1,
ilaenv( 1,
'DGEQRF',
' ', nmax, nmax, -1, -1 ),
552 $
ilaenv( 1,
'DORMQR',
'LT', nmax, nmax, nmax, -1 ),
553 $
ilaenv( 1,
'DORGQR',
' ', nmax, nmax, nmax, -1 ) )
554 nbz =
ilaenv( 1,
'DHGEQZ',
'SII', nmax, 1, nmax, 0 )
555 ns =
ilaenv( 4,
'DHGEQZ',
'SII', nmax, 1, nmax, 0 )
557 lwkopt = 2*nmax + max( 6*nmax, nmax*( nb+1 ),
558 $ ( 2*i1+nmax+1 )*( i1+1 ) )
562 IF( nsizes.LT.0 )
THEN
564 ELSE IF( badnn )
THEN
566 ELSE IF( ntypes.LT.0 )
THEN
568 ELSE IF( thresh.LT.zero )
THEN
570 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
572 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax )
THEN
574 ELSE IF( lwkopt.GT.lwork )
THEN
579 CALL
xerbla(
'DDRVGG', -info )
585 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
588 safmin =
dlamch(
'Safe minimum' )
590 safmin = safmin / ulp
591 safmax = one / safmin
592 CALL
dlabad( safmin, safmax )
606 DO 170 jsize = 1, nsizes
609 rmagn( 2 ) = safmax*ulp / dble( n1 )
610 rmagn( 3 ) = safmin*ulpinv*n1
612 IF( nsizes.NE.1 )
THEN
613 mtypes = min( maxtyp, ntypes )
615 mtypes = min( maxtyp+1, ntypes )
618 DO 160 jtype = 1, mtypes
619 IF( .NOT.dotype( jtype ) )
627 ioldsd( j ) = iseed( j )
659 IF( mtypes.GT.maxtyp )
662 IF( kclass( jtype ).LT.3 )
THEN
666 IF( abs( katype( jtype ) ).EQ.3 )
THEN
667 in = 2*( ( n-1 ) / 2 ) + 1
669 $ CALL
dlaset(
'Full', n, n, zero, zero, a, lda )
673 CALL
dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
674 $ kz2( kazero( jtype ) ), iasign( jtype ),
675 $ rmagn( kamagn( jtype ) ), ulp,
676 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
678 iadd = kadd( kazero( jtype ) )
679 IF( iadd.GT.0 .AND. iadd.LE.n )
680 $ a( iadd, iadd ) = one
684 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
685 in = 2*( ( n-1 ) / 2 ) + 1
687 $ CALL
dlaset(
'Full', n, n, zero, zero, b, lda )
691 CALL
dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
692 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
693 $ rmagn( kbmagn( jtype ) ), one,
694 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
696 iadd = kadd( kbzero( jtype ) )
697 IF( iadd.NE.0 .AND. iadd.LE.n )
698 $ b( iadd, iadd ) = one
700 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
709 q( jr, jc ) =
dlarnd( 3, iseed )
710 z( jr, jc ) =
dlarnd( 3, iseed )
712 CALL
dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
714 work( 2*n+jc ) = sign( one, q( jc, jc ) )
716 CALL
dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
718 work( 3*n+jc ) = sign( one, z( jc, jc ) )
723 work( 3*n ) = sign( one,
dlarnd( 2, iseed ) )
726 work( 4*n ) = sign( one,
dlarnd( 2, iseed ) )
732 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
734 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
738 CALL
dorm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
739 $ lda, work( 2*n+1 ), iinfo )
742 CALL
dorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
743 $ a, lda, work( 2*n+1 ), iinfo )
746 CALL
dorm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
747 $ lda, work( 2*n+1 ), iinfo )
750 CALL
dorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
751 $ b, lda, work( 2*n+1 ), iinfo )
761 a( jr, jc ) = rmagn( kamagn( jtype ) )*
763 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
771 IF( iinfo.NE.0 )
THEN
772 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
782 CALL
dlacpy(
' ', n, n, a, lda, s, lda )
783 CALL
dlacpy(
' ', n, n, b, lda, t, lda )
787 CALL
dgegs(
'V',
'V', n, s, lda, t, lda, alphr1, alphi1,
788 $ beta1, q, ldq, z, ldq, work, lwork, iinfo )
789 IF( iinfo.NE.0 )
THEN
790 WRITE( nounit, fmt = 9999 )
'DGEGS', iinfo, n, jtype,
800 CALL
dget51( 1, n, a, lda, s, lda, q, ldq, z, ldq, work,
802 CALL
dget51( 1, n, b, lda, t, lda, q, ldq, z, ldq, work,
804 CALL
dget51( 3, n, b, lda, t, lda, q, ldq, q, ldq, work,
806 CALL
dget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
816 IF( alphi1( j ).EQ.zero )
THEN
817 temp2 = ( abs( alphr1( j )-s( j, j ) ) /
818 $ max( safmin, abs( alphr1( j ) ), abs( s( j,
819 $ j ) ) )+abs( beta1( j )-t( j, j ) ) /
820 $ max( safmin, abs( beta1( j ) ), abs( t( j,
823 IF( s( j+1, j ).NE.zero )
827 IF( s( j, j-1 ).NE.zero )
831 IF( alphi1( j ).GT.zero )
THEN
836 IF( i1.LE.0 .OR. i1.GE.n )
THEN
838 ELSE IF( i1.LT.n-1 )
THEN
839 IF( s( i1+2, i1+1 ).NE.zero )
841 ELSE IF( i1.GT.1 )
THEN
842 IF( s( i1, i1-1 ).NE.zero )
845 IF( .NOT.ilabad )
THEN
846 CALL
dget53( s( i1, i1 ), lda, t( i1, i1 ), lda,
847 $ beta1( j ), alphr1( j ), alphi1( j ),
849 IF( iinfo.GE.3 )
THEN
850 WRITE( nounit, fmt = 9997 )iinfo, j, n, jtype,
858 temp1 = max( temp1, temp2 )
860 WRITE( nounit, fmt = 9996 )j, n, jtype, ioldsd
869 CALL
dlacpy(
' ', n, n, a, lda, s2, lda )
870 CALL
dlacpy(
' ', n, n, b, lda, t2, lda )
874 CALL
dgegv(
'V',
'V', n, s2, lda, t2, lda, alphr2, alphi2,
875 $ beta2, vl, ldq, vr, ldq, work, lwork, iinfo )
876 IF( iinfo.NE.0 )
THEN
877 WRITE( nounit, fmt = 9999 )
'DGEGV', iinfo, n, jtype,
887 CALL
dget52( .true., n, a, lda, b, lda, vl, ldq, alphr2,
888 $ alphi2, beta2, work, dumma( 1 ) )
889 result( 6 ) = dumma( 1 )
890 IF( dumma( 2 ).GT.thrshn )
THEN
891 WRITE( nounit, fmt = 9998 )
'Left',
'DGEGV', dumma( 2 ),
895 CALL
dget52( .false., n, a, lda, b, lda, vr, ldq, alphr2,
896 $ alphi2, beta2, work, dumma( 1 ) )
897 result( 7 ) = dumma( 1 )
898 IF( dumma( 2 ).GT.thresh )
THEN
899 WRITE( nounit, fmt = 9998 )
'Right',
'DGEGV', dumma( 2 ),
907 IF( alphi2( j ).GT.zero )
THEN
910 ELSE IF( alphi2( j+1 ).GE.zero )
THEN
913 ELSE IF( alphi2( j ).LT.zero )
THEN
916 ELSE IF( alphi2( j-1 ).LE.zero )
THEN
921 WRITE( nounit, fmt = 9996 )j, n, jtype, ioldsd
929 ntestt = ntestt + ntest
934 IF( result( jr ).GE.thresh )
THEN
939 IF( nerrs.EQ.0 )
THEN
940 WRITE( nounit, fmt = 9995 )
'DGG'
944 WRITE( nounit, fmt = 9994 )
945 WRITE( nounit, fmt = 9993 )
946 WRITE( nounit, fmt = 9992 )
'Orthogonal'
950 WRITE( nounit, fmt = 9991 )
'orthogonal',
'''',
951 $
'transpose', (
'''', j = 1, 5 )
955 IF( result( jr ).LT.10000.0d0 )
THEN
956 WRITE( nounit, fmt = 9990 )n, jtype, ioldsd, jr,
959 WRITE( nounit, fmt = 9989 )n, jtype, ioldsd, jr,
970 CALL
alasvm(
'DGG', nounit, nerrs, ntestt, 0 )
973 9999 format(
' DDRVGG: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
974 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
976 9998 format(
' DDRVGG: ', a,
' Eigenvectors from ', a,
' incorrectly ',
977 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 9x,
978 $
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
981 9997 format(
' DDRVGG: DGET53 returned INFO=', i1,
' for eigenvalue ',
982 $ i6,
'.', / 9x,
'N=', i6,
', JTYPE=', i6,
', ISEED=(',
983 $ 3( i5,
',' ), i5,
')' )
985 9996 format(
' DDRVGG: S not in Schur form at eigenvalue ', i6,
'.',
986 $ / 9x,
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
989 9995 format( / 1x, a3,
' -- Real Generalized eigenvalue problem driver'
992 9994 format(
' Matrix types (see DDRVGG for details): ' )
994 9993 format(
' Special Matrices:', 23x,
995 $
'(J''=transposed Jordan block)',
996 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
997 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
998 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
999 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
1000 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1001 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
1002 9992 format(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
1003 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
1004 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
1005 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
1006 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
1007 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
1008 $
'23=(small,large) 24=(small,small) 25=(large,large)',
1009 $ /
' 26=random O(1) matrices.' )
1011 9991 format( /
' Tests performed: (S is Schur, T is triangular, ',
1012 $
'Q and Z are ', a,
',', / 20x,
1013 $
'l and r are the appropriate left and right', / 19x,
1014 $
'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
1015 $
' means ', a,
'.)', /
' 1 = | A - Q S Z', a,
1016 $
' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
1017 $
' | / ( |B| n ulp )', /
' 3 = | I - QQ', a,
1018 $
' | / ( n ulp ) 4 = | I - ZZ', a,
1019 $
' | / ( n ulp )', /
1020 $
' 5 = difference between (alpha,beta) and diagonals of',
1021 $
' (S,T)', /
' 6 = max | ( b A - a B )', a,
1022 $
' l | / const. 7 = max | ( b A - a B ) r | / const.',
1024 9990 format(
' Matrix order=', i5,
', type=', i2,
', seed=',
1025 $ 4( i4,
',' ),
' result ', i3,
' is', 0p, f8.2 )
1026 9989 format(
' Matrix order=', i5,
', type=', i2,
', seed=',
1027 $ 4( i4,
',' ),
' result ', i3,
' is', 1p, d10.3 )