397 SUBROUTINE cdrgev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
398 $ nounit, a, lda, b, s, t, q, ldq, z, qe, ldqe,
399 $ alpha, beta, alpha1, beta1, work, lwork, rwork,
408 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
414 INTEGER ISEED( 4 ), NN( * )
415 REAL RESULT( * ), RWORK( * )
416 COMPLEX A( lda, * ), ALPHA( * ), ALPHA1( * ),
417 $ b( lda, * ), beta( * ), beta1( * ),
418 $ q( ldq, * ), qe( ldqe, * ), s( lda, * ),
419 $ t( lda, * ), work( * ), z( ldq, * )
426 parameter ( zero = 0.0e+0, one = 1.0e+0 )
428 parameter ( czero = ( 0.0e+0, 0.0e+0 ),
429 $ cone = ( 1.0e+0, 0.0e+0 ) )
431 parameter ( maxtyp = 26 )
435 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
436 $ maxwrk, minwrk, mtypes, n, n1, nb, nerrs,
437 $ nmats, nmax, ntestt
438 REAL SAFMAX, SAFMIN, ULP, ULPINV
442 LOGICAL LASIGN( maxtyp ), LBSIGN( maxtyp )
443 INTEGER 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 )
454 EXTERNAL ilaenv, slamch, clarnd
461 INTRINSIC abs, conjg, max, min,
REAL, 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 lasign / 6*.false., .true., .false., 2*.true.,
482 $ 2*.false., 3*.true., .false., .true.,
483 $ 3*.false., 5*.true., .false. /
484 DATA lbsign / 7*.false., .true., 2*.false.,
485 $ 2*.true., 2*.false., .true., .false., .true.,
497 nmax = max( nmax, nn( j ) )
502 IF( nsizes.LT.0 )
THEN
504 ELSE IF( badnn )
THEN
506 ELSE IF( ntypes.LT.0 )
THEN
508 ELSE IF( thresh.LT.zero )
THEN
510 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
512 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax )
THEN
514 ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax )
THEN
526 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
527 minwrk = nmax*( nmax+1 )
528 nb = max( 1, ilaenv( 1,
'CGEQRF',
' ', nmax, nmax, -1, -1 ),
529 $ ilaenv( 1,
'CUNMQR',
'LC', nmax, nmax, nmax, -1 ),
530 $ ilaenv( 1,
'CUNGQR',
' ', nmax, nmax, nmax, -1 ) )
531 maxwrk = max( 2*nmax, nmax*( nb+1 ), nmax*( nmax+1 ) )
535 IF( lwork.LT.minwrk )
539 CALL xerbla(
'CDRGEV', -info )
545 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
548 ulp = slamch(
'Precision' )
549 safmin = slamch(
'Safe minimum' )
550 safmin = safmin / ulp
551 safmax = one / safmin
552 CALL slabad( safmin, safmax )
566 DO 220 jsize = 1, nsizes
569 rmagn( 2 ) = safmax*ulp /
REAL( 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 )
610 IF( mtypes.GT.maxtyp )
613 IF( kclass( jtype ).LT.3 )
THEN
617 IF( abs( katype( jtype ) ).EQ.3 )
THEN
618 in = 2*( ( n-1 ) / 2 ) + 1
620 $
CALL claset(
'Full', n, n, czero, czero, a, lda )
624 CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
625 $ kz2( kazero( jtype ) ), lasign( jtype ),
626 $ rmagn( kamagn( jtype ) ), ulp,
627 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
629 iadd = kadd( kazero( jtype ) )
630 IF( iadd.GT.0 .AND. iadd.LE.n )
631 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
635 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
636 in = 2*( ( n-1 ) / 2 ) + 1
638 $
CALL claset(
'Full', n, n, czero, czero, b, lda )
642 CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
643 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
644 $ rmagn( kbmagn( jtype ) ), one,
645 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
647 iadd = kadd( kbzero( jtype ) )
648 IF( iadd.NE.0 .AND. iadd.LE.n )
649 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
651 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
660 q( jr, jc ) = clarnd( 3, iseed )
661 z( jr, jc ) = clarnd( 3, iseed )
663 CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
665 work( 2*n+jc ) = sign( one,
REAL( Q( JC, JC ) ) )
667 CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
669 work( 3*n+jc ) = sign( one,
REAL( Z( JC, JC ) ) )
672 ctemp = clarnd( 3, iseed )
675 work( 3*n ) = ctemp / abs( ctemp )
676 ctemp = clarnd( 3, iseed )
679 work( 4*n ) = ctemp / abs( ctemp )
685 a( jr, jc ) = work( 2*n+jr )*
686 $ conjg( work( 3*n+jc ) )*
688 b( jr, jc ) = work( 2*n+jr )*
689 $ conjg( work( 3*n+jc ) )*
693 CALL cunm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
694 $ lda, work( 2*n+1 ), ierr )
697 CALL cunm2r(
'R',
'C', n, n, n-1, z, ldq, work( n+1 ),
698 $ a, lda, work( 2*n+1 ), ierr )
701 CALL cunm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
702 $ lda, work( 2*n+1 ), ierr )
705 CALL cunm2r(
'R',
'C', n, n, n-1, z, ldq, work( n+1 ),
706 $ b, lda, work( 2*n+1 ), ierr )
716 a( jr, jc ) = rmagn( kamagn( jtype ) )*
718 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
727 WRITE( nounit, fmt = 9999 )
'Generator', ierr, n, jtype,
741 CALL clacpy(
' ', n, n, a, lda, s, lda )
742 CALL clacpy(
' ', n, n, b, lda, t, lda )
743 CALL cggev(
'V',
'V', n, s, lda, t, lda, alpha, beta, q,
744 $ ldq, z, ldq, work, lwork, rwork, ierr )
745 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
747 WRITE( nounit, fmt = 9999 )
'CGGEV1', ierr, n, jtype,
755 CALL cget52( .true., n, a, lda, b, lda, q, ldq, alpha, beta,
756 $ work, rwork, result( 1 ) )
757 IF( result( 2 ).GT.thresh )
THEN
758 WRITE( nounit, fmt = 9998 )
'Left',
'CGGEV1',
759 $ result( 2 ), n, jtype, ioldsd
764 CALL cget52( .false., n, a, lda, b, lda, z, ldq, alpha,
765 $ beta, work, rwork, result( 3 ) )
766 IF( result( 4 ).GT.thresh )
THEN
767 WRITE( nounit, fmt = 9998 )
'Right',
'CGGEV1',
768 $ result( 4 ), n, jtype, ioldsd
773 CALL clacpy(
' ', n, n, a, lda, s, lda )
774 CALL clacpy(
' ', n, n, b, lda, t, lda )
775 CALL cggev(
'N',
'N', n, s, lda, t, lda, alpha1, beta1, q,
776 $ ldq, z, ldq, work, lwork, rwork, ierr )
777 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
779 WRITE( nounit, fmt = 9999 )
'CGGEV2', ierr, n, jtype,
786 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
787 $ beta1( j ) )result( 5 ) = ulpinv
793 CALL clacpy(
' ', n, n, a, lda, s, lda )
794 CALL clacpy(
' ', n, n, b, lda, t, lda )
795 CALL cggev(
'V',
'N', n, s, lda, t, lda, alpha1, beta1, qe,
796 $ ldqe, z, ldq, work, lwork, rwork, ierr )
797 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
799 WRITE( nounit, fmt = 9999 )
'CGGEV3', ierr, n, jtype,
806 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
807 $ beta1( j ) )result( 6 ) = ulpinv
812 IF( q( j, jc ).NE.qe( j, jc ) )
813 $ result( 6 ) = ulpinv
820 CALL clacpy(
' ', n, n, a, lda, s, lda )
821 CALL clacpy(
' ', n, n, b, lda, t, lda )
822 CALL cggev(
'N',
'V', n, s, lda, t, lda, alpha1, beta1, q,
823 $ ldq, qe, ldqe, work, lwork, rwork, ierr )
824 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
826 WRITE( nounit, fmt = 9999 )
'CGGEV4', ierr, n, jtype,
833 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
834 $ beta1( j ) )result( 7 ) = ulpinv
839 IF( z( j, jc ).NE.qe( j, jc ) )
840 $ result( 7 ) = ulpinv
853 IF( result( jr ).GE.thresh )
THEN
858 IF( nerrs.EQ.0 )
THEN
859 WRITE( nounit, fmt = 9997 )
'CGV'
863 WRITE( nounit, fmt = 9996 )
864 WRITE( nounit, fmt = 9995 )
865 WRITE( nounit, fmt = 9994 )
'Orthogonal'
869 WRITE( nounit, fmt = 9993 )
873 IF( result( jr ).LT.10000.0 )
THEN
874 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
877 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
888 CALL alasvm(
'CGV', nounit, nerrs, ntestt, 0 )
894 9999
FORMAT(
' CDRGEV: ', a,
' returned INFO=', i6,
'.', / 3x,
'N=',
895 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
897 9998
FORMAT(
' CDRGEV: ', a,
' Eigenvectors from ', a,
' incorrectly ',
898 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 3x,
899 $
'N=', i4,
', JTYPE=', i3,
', ISEED=(', 3( i4,
',' ), i5,
902 9997
FORMAT( / 1x, a3,
' -- Complex Generalized eigenvalue problem ',
905 9996
FORMAT(
' Matrix types (see CDRGEV for details): ' )
907 9995
FORMAT(
' Special Matrices:', 23x,
908 $
'(J''=transposed Jordan block)',
909 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
910 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
911 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
912 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
913 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
914 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
915 9994
FORMAT(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
916 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
917 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
918 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
919 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
920 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
921 $
'23=(small,large) 24=(small,small) 25=(large,large)',
922 $ /
' 26=random O(1) matrices.' )
924 9993
FORMAT( /
' Tests performed: ',
925 $ /
' 1 = max | ( b A - a B )''*l | / const.,',
926 $ /
' 2 = | |VR(i)| - 1 | / ulp,',
927 $ /
' 3 = max | ( b A - a B )*r | / const.',
928 $ /
' 4 = | |VL(i)| - 1 | / ulp,',
929 $ /
' 5 = 0 if W same no matter if r or l computed,',
930 $ /
' 6 = 0 if l same no matter if l computed,',
931 $ /
' 7 = 0 if r same no matter if r computed,', / 1x )
932 9992
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
933 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
934 9991
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
935 $ 4( i4,
',' ),
' result ', i2,
' is', 1p, e10.3 )
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine cunm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
CUNM2R multiplies a general matrix by the unitary matrix from a QR factorization determined by cgeqrf...
subroutine cget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA, WORK, RWORK, RESULT)
CGET52
subroutine clatm4(ITYPE, N, NZ1, NZ2, RSIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
CLATM4
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine cdrgev(NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH, NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE, ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK, RESULT, INFO)
CDRGEV
subroutine cggev(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...