395 SUBROUTINE cdrgev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
396 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
397 $ ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK,
405 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
411 INTEGER ISEED( 4 ), NN( * )
412 REAL RESULT( * ), RWORK( * )
413 COMPLEX A( LDA, * ), ALPHA( * ), ALPHA1( * ),
414 $ b( lda, * ), beta( * ), beta1( * ),
415 $ q( ldq, * ), qe( ldqe, * ), s( lda, * ),
416 $ t( lda, * ), work( * ), z( ldq, * )
423 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
425 parameter( czero = ( 0.0e+0, 0.0e+0 ),
426 $ cone = ( 1.0e+0, 0.0e+0 ) )
428 parameter( maxtyp = 26 )
432 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
433 $ MAXWRK, MINWRK, MTYPES, N, N1, NB, NERRS,
434 $ nmats, nmax, ntestt
435 REAL SAFMAX, SAFMIN, ULP, ULPINV
439 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
440 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
441 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
442 $ kbmagn( maxtyp ), kbtype( maxtyp ),
443 $ kbzero( maxtyp ), kclass( maxtyp ),
444 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
451 EXTERNAL ilaenv, slamch, clarnd
458 INTRINSIC abs, conjg, max, min, real, sign
461 DATA kclass / 15*1, 10*2, 1*3 /
462 DATA kz1 / 0, 1, 2, 1, 3, 3 /
463 DATA kz2 / 0, 0, 1, 2, 1, 1 /
464 DATA kadd / 0, 0, 0, 0, 3, 2 /
465 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
466 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
467 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
468 $ 1, 1, -4, 2, -4, 8*8, 0 /
469 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
471 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
473 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
475 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
477 DATA ktrian / 16*0, 10*1 /
478 DATA lasign / 6*.false., .true., .false., 2*.true.,
479 $ 2*.false., 3*.true., .false., .true.,
480 $ 3*.false., 5*.true., .false. /
481 DATA lbsign / 7*.false., .true., 2*.false.,
482 $ 2*.true., 2*.false., .true., .false., .true.,
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
511 ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax )
THEN
523 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
524 minwrk = nmax*( nmax+1 )
525 nb = max( 1, ilaenv( 1,
'CGEQRF',
' ', nmax, nmax, -1, -1 ),
526 $ ilaenv( 1,
'CUNMQR',
'LC', nmax, nmax, nmax, -1 ),
527 $ ilaenv( 1,
'CUNGQR',
' ', nmax, nmax, nmax, -1 ) )
528 maxwrk = max( 2*nmax, nmax*( nb+1 ), nmax*( nmax+1 ) )
532 IF( lwork.LT.minwrk )
536 CALL xerbla(
'CDRGEV', -info )
542 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
545 ulp = slamch(
'Precision' )
546 safmin = slamch(
'Safe minimum' )
547 safmin = safmin / ulp
548 safmax = one / safmin
562 DO 220 jsize = 1, nsizes
565 rmagn( 2 ) = safmax*ulp / real( n1 )
566 rmagn( 3 ) = safmin*ulpinv*n1
568 IF( nsizes.NE.1 )
THEN
569 mtypes = min( maxtyp, ntypes )
571 mtypes = min( maxtyp+1, ntypes )
574 DO 210 jtype = 1, mtypes
575 IF( .NOT.dotype( jtype ) )
582 ioldsd( j ) = iseed( j )
606 IF( mtypes.GT.maxtyp )
609 IF( kclass( jtype ).LT.3 )
THEN
613 IF( abs( katype( jtype ) ).EQ.3 )
THEN
614 in = 2*( ( n-1 ) / 2 ) + 1
616 $
CALL claset(
'Full', n, n, czero, czero, a, lda )
620 CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
621 $ kz2( kazero( jtype ) ), lasign( jtype ),
622 $ rmagn( kamagn( jtype ) ), ulp,
623 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
625 iadd = kadd( kazero( jtype ) )
626 IF( iadd.GT.0 .AND. iadd.LE.n )
627 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
631 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
632 in = 2*( ( n-1 ) / 2 ) + 1
634 $
CALL claset(
'Full', n, n, czero, czero, b, lda )
638 CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
639 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
640 $ rmagn( kbmagn( jtype ) ), one,
641 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
643 iadd = kadd( kbzero( jtype ) )
644 IF( iadd.NE.0 .AND. iadd.LE.n )
645 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
647 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
656 q( jr, jc ) = clarnd( 3, iseed )
657 z( jr, jc ) = clarnd( 3, iseed )
659 CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
661 work( 2*n+jc ) = sign( one, real( q( jc, jc ) ) )
663 CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
665 work( 3*n+jc ) = sign( one, real( z( jc, jc ) ) )
668 ctemp = clarnd( 3, iseed )
671 work( 3*n ) = ctemp / abs( ctemp )
672 ctemp = clarnd( 3, iseed )
675 work( 4*n ) = ctemp / abs( ctemp )
681 a( jr, jc ) = work( 2*n+jr )*
682 $ conjg( work( 3*n+jc ) )*
684 b( jr, jc ) = work( 2*n+jr )*
685 $ conjg( work( 3*n+jc ) )*
689 CALL cunm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
690 $ lda, work( 2*n+1 ), ierr )
693 CALL cunm2r(
'R',
'C', n, n, n-1, z, ldq, work( n+1 ),
694 $ a, lda, work( 2*n+1 ), ierr )
697 CALL cunm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
698 $ lda, work( 2*n+1 ), ierr )
701 CALL cunm2r(
'R',
'C', n, n, n-1, z, ldq, work( n+1 ),
702 $ b, lda, work( 2*n+1 ), ierr )
712 a( jr, jc ) = rmagn( kamagn( jtype ) )*
714 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
723 WRITE( nounit, fmt = 9999 )
'Generator', ierr, n, jtype,
737 CALL clacpy(
' ', n, n, a, lda, s, lda )
738 CALL clacpy(
' ', n, n, b, lda, t, lda )
739 CALL cggev(
'V',
'V', n, s, lda, t, lda, alpha, beta, q,
740 $ ldq, z, ldq, work, lwork, rwork, ierr )
741 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
743 WRITE( nounit, fmt = 9999 )
'CGGEV1', ierr, n, jtype,
751 CALL cget52( .true., n, a, lda, b, lda, q, ldq, alpha, beta,
752 $ work, rwork, result( 1 ) )
753 IF( result( 2 ).GT.thresh )
THEN
754 WRITE( nounit, fmt = 9998 )
'Left',
'CGGEV1',
755 $ result( 2 ), n, jtype, ioldsd
760 CALL cget52( .false., n, a, lda, b, lda, z, ldq, alpha,
761 $ beta, work, rwork, result( 3 ) )
762 IF( result( 4 ).GT.thresh )
THEN
763 WRITE( nounit, fmt = 9998 )
'Right',
'CGGEV1',
764 $ result( 4 ), n, jtype, ioldsd
769 CALL clacpy(
' ', n, n, a, lda, s, lda )
770 CALL clacpy(
' ', n, n, b, lda, t, lda )
771 CALL cggev(
'N',
'N', n, s, lda, t, lda, alpha1, beta1, q,
772 $ ldq, z, ldq, work, lwork, rwork, ierr )
773 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
775 WRITE( nounit, fmt = 9999 )
'CGGEV2', ierr, n, jtype,
782 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
783 $ beta1( j ) )result( 5 ) = ulpinv
789 CALL clacpy(
' ', n, n, a, lda, s, lda )
790 CALL clacpy(
' ', n, n, b, lda, t, lda )
791 CALL cggev(
'V',
'N', n, s, lda, t, lda, alpha1, beta1, qe,
792 $ ldqe, z, ldq, work, lwork, rwork, ierr )
793 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
795 WRITE( nounit, fmt = 9999 )
'CGGEV3', ierr, n, jtype,
802 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
803 $ beta1( j ) )result( 6 ) = ulpinv
808 IF( q( j, jc ).NE.qe( j, jc ) )
809 $ result( 6 ) = ulpinv
816 CALL clacpy(
' ', n, n, a, lda, s, lda )
817 CALL clacpy(
' ', n, n, b, lda, t, lda )
818 CALL cggev(
'N',
'V', n, s, lda, t, lda, alpha1, beta1, q,
819 $ ldq, qe, ldqe, work, lwork, rwork, ierr )
820 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
822 WRITE( nounit, fmt = 9999 )
'CGGEV4', ierr, n, jtype,
829 IF( alpha( j ).NE.alpha1( j ) .OR. beta( j ).NE.
830 $ beta1( j ) )result( 7 ) = ulpinv
835 IF( z( j, jc ).NE.qe( j, jc ) )
836 $ result( 7 ) = ulpinv
849 IF( result( jr ).GE.thresh )
THEN
854 IF( nerrs.EQ.0 )
THEN
855 WRITE( nounit, fmt = 9997 )
'CGV'
859 WRITE( nounit, fmt = 9996 )
860 WRITE( nounit, fmt = 9995 )
861 WRITE( nounit, fmt = 9994 )
'Orthogonal'
865 WRITE( nounit, fmt = 9993 )
869 IF( result( jr ).LT.10000.0 )
THEN
870 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
873 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
884 CALL alasvm(
'CGV', nounit, nerrs, ntestt, 0 )
890 9999
FORMAT(
' CDRGEV: ', a,
' returned INFO=', i6,
'.', / 3x,
'N=',
891 $ i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ), i5,
')' )
893 9998
FORMAT(
' CDRGEV: ', a,
' Eigenvectors from ', a,
' incorrectly ',
894 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 3x,
895 $
'N=', i4,
', JTYPE=', i3,
', ISEED=(', 3( i4,
',' ), i5,
898 9997
FORMAT( / 1x, a3,
' -- Complex Generalized eigenvalue problem ',
901 9996
FORMAT(
' Matrix types (see CDRGEV for details): ' )
903 9995
FORMAT(
' Special Matrices:', 23x,
904 $
'(J''=transposed Jordan block)',
905 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
906 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
907 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
908 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
909 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
910 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
911 9994
FORMAT(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
912 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
913 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
914 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
915 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
916 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
917 $
'23=(small,large) 24=(small,small) 25=(large,large)',
918 $ /
' 26=random O(1) matrices.' )
920 9993
FORMAT( /
' Tests performed: ',
921 $ /
' 1 = max | ( b A - a B )''*l | / const.,',
922 $ /
' 2 = | |VR(i)| - 1 | / ulp,',
923 $ /
' 3 = max | ( b A - a B )*r | / const.',
924 $ /
' 4 = | |VL(i)| - 1 | / ulp,',
925 $ /
' 5 = 0 if W same no matter if r or l computed,',
926 $ /
' 6 = 0 if l same no matter if l computed,',
927 $ /
' 7 = 0 if r same no matter if r computed,', / 1x )
928 9992
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
929 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
930 9991
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
931 $ 4( i4,
',' ),
' result ', i2,
' is', 1p, e10.3 )