406 SUBROUTINE ddrgev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
407 $ nounit, a, lda, b, s, t, q, ldq, z, qe, ldqe,
408 $ alphar, alphai, beta, alphr1, alphi1, beta1,
409 $ work, lwork, result, info )
417 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
419 DOUBLE PRECISION THRESH
423 INTEGER ISEED( 4 ), NN( * )
424 DOUBLE PRECISION A( lda, * ), ALPHAI( * ), ALPHAR( * ),
425 $ alphi1( * ), alphr1( * ), b( lda, * ),
426 $ beta( * ), beta1( * ), q( ldq, * ),
427 $ qe( ldqe, * ), result( * ), s( lda, * ),
428 $ t( lda, * ), work( * ), z( ldq, * )
434 DOUBLE PRECISION ZERO, ONE
435 parameter ( zero = 0.0d+0, one = 1.0d+0 )
437 parameter ( maxtyp = 26 )
441 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
442 $ maxwrk, minwrk, mtypes, n, n1, nerrs, nmats,
444 DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV
447 INTEGER IASIGN( maxtyp ), IBSIGN( maxtyp ),
448 $ ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
449 $ katype( maxtyp ), kazero( maxtyp ),
450 $ kbmagn( maxtyp ), kbtype( maxtyp ),
451 $ kbzero( maxtyp ), kclass( maxtyp ),
452 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
453 DOUBLE PRECISION RMAGN( 0: 3 )
457 DOUBLE PRECISION DLAMCH, DLARND
458 EXTERNAL ilaenv, dlamch, dlarnd
465 INTRINSIC abs, dble, max, min, sign
468 DATA kclass / 15*1, 10*2, 1*3 /
469 DATA kz1 / 0, 1, 2, 1, 3, 3 /
470 DATA kz2 / 0, 0, 1, 2, 1, 1 /
471 DATA kadd / 0, 0, 0, 0, 3, 2 /
472 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
473 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
474 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
475 $ 1, 1, -4, 2, -4, 8*8, 0 /
476 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
478 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
480 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
482 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
484 DATA ktrian / 16*0, 10*1 /
485 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
487 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
498 nmax = max( nmax, nn( j ) )
503 IF( nsizes.LT.0 )
THEN
505 ELSE IF( badnn )
THEN
507 ELSE IF( ntypes.LT.0 )
THEN
509 ELSE IF( thresh.LT.zero )
THEN
511 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
513 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax )
THEN
515 ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax )
THEN
527 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
528 minwrk = max( 1, 8*nmax, nmax*( nmax+1 ) )
529 maxwrk = 7*nmax + nmax*ilaenv( 1,
'DGEQRF',
' ', nmax, 1, nmax,
531 maxwrk = max( maxwrk, nmax*( nmax+1 ) )
535 IF( lwork.LT.minwrk )
539 CALL xerbla(
'DDRGEV', -info )
545 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
548 safmin = dlamch(
'Safe minimum' )
549 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
550 safmin = safmin / ulp
551 safmax = one / safmin
552 CALL dlabad( safmin, safmax )
566 DO 220 jsize = 1, nsizes
569 rmagn( 2 ) = safmax*ulp / dble( n1 )
570 rmagn( 3 ) = safmin*ulpinv*n1
572 IF( nsizes.NE.1 )
THEN
573 mtypes = min( maxtyp, ntypes )
575 mtypes = min( maxtyp+1, ntypes )
578 DO 210 jtype = 1, mtypes
579 IF( .NOT.dotype( jtype ) )
586 ioldsd( j ) = iseed( j )
612 IF( mtypes.GT.maxtyp )
615 IF( kclass( jtype ).LT.3 )
THEN
619 IF( abs( katype( jtype ) ).EQ.3 )
THEN
620 in = 2*( ( n-1 ) / 2 ) + 1
622 $
CALL dlaset(
'Full', n, n, zero, zero, a, lda )
626 CALL dlatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
627 $ kz2( kazero( jtype ) ), iasign( jtype ),
628 $ rmagn( kamagn( jtype ) ), ulp,
629 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
631 iadd = kadd( kazero( jtype ) )
632 IF( iadd.GT.0 .AND. iadd.LE.n )
633 $ a( iadd, iadd ) = one
637 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
638 in = 2*( ( n-1 ) / 2 ) + 1
640 $
CALL dlaset(
'Full', n, n, zero, zero, b, lda )
644 CALL dlatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
645 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
646 $ rmagn( kbmagn( jtype ) ), one,
647 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
649 iadd = kadd( kbzero( jtype ) )
650 IF( iadd.NE.0 .AND. iadd.LE.n )
651 $ b( iadd, iadd ) = one
653 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
662 q( jr, jc ) = dlarnd( 3, iseed )
663 z( jr, jc ) = dlarnd( 3, iseed )
665 CALL dlarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
667 work( 2*n+jc ) = sign( one, q( jc, jc ) )
669 CALL dlarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
671 work( 3*n+jc ) = sign( one, z( jc, jc ) )
676 work( 3*n ) = sign( one, dlarnd( 2, iseed ) )
679 work( 4*n ) = sign( one, dlarnd( 2, iseed ) )
685 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
687 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
691 CALL dorm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
692 $ lda, work( 2*n+1 ), ierr )
695 CALL dorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
696 $ a, lda, work( 2*n+1 ), ierr )
699 CALL dorm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
700 $ lda, work( 2*n+1 ), ierr )
703 CALL dorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
704 $ b, lda, work( 2*n+1 ), ierr )
714 a( jr, jc ) = rmagn( kamagn( jtype ) )*
716 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
725 WRITE( nounit, fmt = 9999 )
'Generator', ierr, n, jtype,
739 CALL dlacpy(
' ', n, n, a, lda, s, lda )
740 CALL dlacpy(
' ', n, n, b, lda, t, lda )
741 CALL dggev(
'V',
'V', n, s, lda, t, lda, alphar, alphai,
742 $ beta, q, ldq, z, ldq, work, lwork, ierr )
743 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
745 WRITE( nounit, fmt = 9999 )
'DGGEV1', ierr, n, jtype,
753 CALL dget52( .true., n, a, lda, b, lda, q, ldq, alphar,
754 $ alphai, beta, work, result( 1 ) )
755 IF( result( 2 ).GT.thresh )
THEN
756 WRITE( nounit, fmt = 9998 )
'Left',
'DGGEV1',
757 $ result( 2 ), n, jtype, ioldsd
762 CALL dget52( .false., n, a, lda, b, lda, z, ldq, alphar,
763 $ alphai, beta, work, result( 3 ) )
764 IF( result( 4 ).GT.thresh )
THEN
765 WRITE( nounit, fmt = 9998 )
'Right',
'DGGEV1',
766 $ result( 4 ), n, jtype, ioldsd
771 CALL dlacpy(
' ', n, n, a, lda, s, lda )
772 CALL dlacpy(
' ', n, n, b, lda, t, lda )
773 CALL dggev(
'N',
'N', n, s, lda, t, lda, alphr1, alphi1,
774 $ beta1, q, ldq, z, ldq, work, lwork, ierr )
775 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
777 WRITE( nounit, fmt = 9999 )
'DGGEV2', ierr, n, jtype,
784 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
785 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 5 )
792 CALL dlacpy(
' ', n, n, a, lda, s, lda )
793 CALL dlacpy(
' ', n, n, b, lda, t, lda )
794 CALL dggev(
'V',
'N', n, s, lda, t, lda, alphr1, alphi1,
795 $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
796 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
798 WRITE( nounit, fmt = 9999 )
'DGGEV3', ierr, n, jtype,
805 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
806 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 6 )
812 IF( q( j, jc ).NE.qe( j, jc ) )
813 $ result( 6 ) = ulpinv
820 CALL dlacpy(
' ', n, n, a, lda, s, lda )
821 CALL dlacpy(
' ', n, n, b, lda, t, lda )
822 CALL dggev(
'N',
'V', n, s, lda, t, lda, alphr1, alphi1,
823 $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
824 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
826 WRITE( nounit, fmt = 9999 )
'DGGEV4', ierr, n, jtype,
833 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
834 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )result( 7 )
840 IF( z( j, jc ).NE.qe( j, jc ) )
841 $ result( 7 ) = ulpinv
854 IF( result( jr ).GE.thresh )
THEN
859 IF( nerrs.EQ.0 )
THEN
860 WRITE( nounit, fmt = 9997 )
'DGV'
864 WRITE( nounit, fmt = 9996 )
865 WRITE( nounit, fmt = 9995 )
866 WRITE( nounit, fmt = 9994 )
'Orthogonal'
870 WRITE( nounit, fmt = 9993 )
874 IF( result( jr ).LT.10000.0d0 )
THEN
875 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
878 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
889 CALL alasvm(
'DGV', nounit, nerrs, ntestt, 0 )
895 9999
FORMAT(
' DDRGEV: ', a,
' returned INFO=', i6,
'.', / 3x,
'N=',
896 $ i6,
', JTYPE=', i6,
', ISEED=(', 4( i4,
',' ), i5,
')' )
898 9998
FORMAT(
' DDRGEV: ', a,
' Eigenvectors from ', a,
' incorrectly ',
899 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 3x,
900 $
'N=', i4,
', JTYPE=', i3,
', ISEED=(', 4( i4,
',' ), i5,
903 9997
FORMAT( / 1x, a3,
' -- Real Generalized eigenvalue problem driver'
906 9996
FORMAT(
' Matrix types (see DDRGEV for details): ' )
908 9995
FORMAT(
' Special Matrices:', 23x,
909 $
'(J''=transposed Jordan block)',
910 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
911 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
912 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
913 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
914 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
915 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
916 9994
FORMAT(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
917 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
918 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
919 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
920 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
921 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
922 $
'23=(small,large) 24=(small,small) 25=(large,large)',
923 $ /
' 26=random O(1) matrices.' )
925 9993
FORMAT( /
' Tests performed: ',
926 $ /
' 1 = max | ( b A - a B )''*l | / const.,',
927 $ /
' 2 = | |VR(i)| - 1 | / ulp,',
928 $ /
' 3 = max | ( b A - a B )*r | / const.',
929 $ /
' 4 = | |VL(i)| - 1 | / ulp,',
930 $ /
' 5 = 0 if W same no matter if r or l computed,',
931 $ /
' 6 = 0 if l same no matter if l computed,',
932 $ /
' 7 = 0 if r same no matter if r computed,', / 1x )
933 9992
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
934 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
935 9991
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
936 $ 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 dggev(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
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 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 dlabad(SMALL, LARGE)
DLABAD
subroutine dlarfg(N, ALPHA, X, INCX, TAU)
DLARFG generates an elementary reflector (Householder matrix).
subroutine ddrgev(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1, WORK, LWORK, RESULT, INFO)
DDRGEV
subroutine dget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
DGET52