395 SUBROUTINE cdrgev3( 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,
398 $ RWORK, RESULT, INFO )
405 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
411 INTEGER ISEED( 4 ), NN( * )
412 REAL RESULT( * ), RWORK( * )
413 COMPLEX ( LDA, * ), ( * ), 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
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(
'CDRGEV3', -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
549 CALL slabad( safmin, safmax )
563 DO 220 jsize = 1, nsizes
566 rmagn( 2 ) = safmax*ulp / real( n1 )
567 rmagn( 3 ) = safmin*ulpinv*n1
569 IF( nsizes.NE.1 )
THEN
570 mtypes =
min( maxtyp, ntypes )
572 mtypes =
min( maxtyp+1, ntypes )
575 DO 210 jtype = 1, mtypes
576 IF( .NOT.dotype( jtype ) )
583 ioldsd( j ) = iseed( j )
607 IF( mtypes.GT.maxtyp )
610 IF( kclass( jtype ).LT.3 )
THEN
614 IF( abs( katype( jtype ) ).EQ.3 )
THEN
615 in = 2*( ( n-1 ) / 2 ) + 1
617 $
CALL claset(
'Full', n, n, czero, czero, a, lda )
621 CALL clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
622 $ kz2( kazero( jtype ) ), lasign( jtype ),
623 $ rmagn( kamagn( jtype ) ), ulp,
624 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
626 iadd = kadd( kazero( jtype ) )
627 IF( iadd.GT.0 .AND. iadd.LE.n )
628 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
632 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
633 in = 2*( ( n-1 ) / 2 ) + 1
639 CALL clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
640 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
641 $ rmagn( kbmagn( jtype ) ), one,
642 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
644 iadd = kadd( kbzero( jtype ) )
645 IF( iadd.NE.0 .AND. iadd.LE.n )
646 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
648 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
657 q( jr, jc ) = clarnd( 3, iseed )
658 z( jr, jc ) = clarnd( 3, iseed )
660 CALL clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
662 work( 2*n+jc ) = sign( one, real( q( jc, jc ) ) )
664 CALL clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
666 work( 3*n+jc ) = sign( one, real( z( jc, jc ) ) )
669 ctemp = clarnd( 3, iseed )
672 work( 3*n ) = ctemp / abs( ctemp )
673 ctemp = clarnd( 3, iseed )
676 work( 4*n ) = ctemp / abs( ctemp )
682 a( jr, jc ) = work( 2*n+jr )*
683 $ conjg( work( 3*n+jc ) )*
685 b( jr, jc ) = work( 2*n+jr )*
686 $ conjg( work( 3*n+jc ) )*
690 CALL cunm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
691 $ lda, work( 2*n+1 ), ierr )
694 CALL cunm2r(
'R',
'C', n, n, n-1, z, ldq, work( n+1 ),
695 $ a, lda, work( 2*n+1 ), ierr )
698 CALL cunm2r( 'l
', 'n
', N, N, N-1, Q, LDQ, WORK, B,
699 $ LDA, WORK( 2*N+1 ), IERR )
702 CALL CUNM2R( 'r
', 'c
', N, N, N-1, Z, LDQ, WORK( N+1 ),
703 $ B, LDA, WORK( 2*N+1 ), IERR )
713 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
715 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
724 WRITE( NOUNIT, FMT = 9999 )'generator
', IERR, N, JTYPE,
738 CALL XLAENV( 12, 10 )
739 CALL XLAENV( 13, 12 )
740 CALL XLAENV( 14, 13 )
742 CALL XLAENV( 17, 10 )
746 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
747 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
748 CALL CGGEV3( 'v
', 'v
', N, S, LDA, T, LDA, ALPHA, BETA, Q,
749 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
750.NE..AND..NE.
IF( IERR0 IERRN+1 ) THEN
752 WRITE( NOUNIT, FMT = 9999 )'cggev31
', IERR, N, JTYPE,
760 CALL CGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHA, BETA,
761 $ WORK, RWORK, RESULT( 1 ) )
762.GT.
IF( RESULT( 2 )THRESH ) THEN
763 WRITE( NOUNIT, FMT = 9998 )'left
', 'cggev31
',
764 $ RESULT( 2 ), N, JTYPE, IOLDSD
769 CALL CGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHA,
770 $ BETA, WORK, RWORK, RESULT( 3 ) )
771.GT.
IF( RESULT( 4 )THRESH ) THEN
772 WRITE( NOUNIT, FMT = 9998 )'right
', 'cggev31
',
773 $ RESULT( 4 ), N, JTYPE, IOLDSD
778 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
779 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
780 CALL CGGEV3( 'n
', 'n
', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
781 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
782.NE..AND..NE.
IF( IERR0 IERRN+1 ) THEN
784 WRITE( NOUNIT, FMT = 9999 )'cggev32
', IERR, N, JTYPE,
791.NE..OR..NE.
IF( ALPHA( J )ALPHA1( J ) BETA( J )
792 $ BETA1( J ) ) RESULT( 5 ) = ULPINV
798 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
799 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
800 CALL CGGEV3( 'v
', 'n
', N, S, LDA, T, LDA, ALPHA1, BETA1, QE,
801 $ LDQE, Z, LDQ, WORK, LWORK, RWORK, IERR )
802.NE..AND..NE.
IF( IERR0 IERRN+1 ) THEN
804 WRITE( NOUNIT, FMT = 9999 )'cggev33
', IERR, N, JTYPE,
812.NE..OR.
IF( ALPHA( J )ALPHA1( J )
813.NE.
$ BETA( J )BETA1( J ) ) THEN
820.NE.
IF( Q( J, JC )QE( J, JC ) ) THEN
829 CALL CLACPY( ' ', N, N, A, LDA, S, LDA )
830 CALL CLACPY( ' ', N, N, B, LDA, T, LDA )
831 CALL CGGEV3( 'n
', 'v
', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
832 $ LDQ, QE, LDQE, WORK, LWORK, RWORK, IERR )
833.NE..AND..NE.
IF( IERR0 IERRN+1 ) THEN
835 WRITE( NOUNIT, FMT = 9999 )'cggev34
', IERR, N, JTYPE,
842.NE..OR..NE.
IF( ALPHA( J )ALPHA1( J ) BETA( J )
843 $ BETA1( J ) )RESULT( 7 ) = ULPINV
848.NE.
IF( Z( J, JC )QE( J, JC ) )
849 $ RESULT( 7 ) = ULPINV
862.GE.
IF( RESULT( JR )THRESH ) THEN
867.EQ.
IF( NERRS0 ) THEN
868 WRITE( NOUNIT, FMT = 9997 )'cgv
'
872 WRITE( NOUNIT, FMT = 9996 )
873 WRITE( NOUNIT, FMT = 9995 )
874 WRITE( NOUNIT, FMT = 9994 )'orthogonal
'
878 WRITE( NOUNIT, FMT = 9993 )
882.LT.
IF( RESULT( JR )10000.0 ) THEN
883 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
886 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
897 CALL ALASVM( 'cgv3
', NOUNIT, NERRS, NTESTT, 0 )
903 9999 FORMAT( ' cdrgev3:
', A, ' returned info=
', I6, '.
', / 3X, 'n=
',
904 $ I6, ', jtype=
', I6, ', iseed=(
', 3( I5, ',
' ), I5, ')
' )
906 9998 FORMAT( ' cdrgev3:
', A, ' eigenvectors from
', A,
907 $ ' incorrectly normalized.
', / ' bits of error=
', 0P, G10.3,
908 $ ',
', 3X, 'n=
', I4, ', jtype=
', I3, ', iseed=(
',
909 $ 3( I4, ',
' ), I5, ')
' )
911 9997 FORMAT( / 1X, A3, ' --
Complex Generalized eigenvalue problem
',
914 9996 FORMAT( ' Matrix types (see CDRGEV3 for details):
' )
916 9995 FORMAT( ' Matrices:
', 23X,
917 $ '(J
''=transposed jordan block)
',
918 $ / ' 1=(0,0) 2=(i,0) 3=(0,i) 4=(i,i) 5=(j
'',j
'')
',
919 $ '6=(diag(j
'',i), diag(i,j
''))
', / ' diagonal matrices: (
',
920 $ 'd=diag(0,1,2,...) )
', / ' 7=(d,i) 9=(large*d, small*i
',
921 $ ') 11=(large*i, small*d) 13=(large*d, large*i)
', /
922 $ ' 8=(i,d) 10=(small*d, large*i) 12=(small*i, large*d)
',
923 $ ' 14=(small*d, small*i)
', / ' 15=(d, reversed d)
' )
924 9994 FORMAT( ' matrices rotated by random
', A, ' matrices u, v:
',
925 $ / ' 16=transposed jordan blocks 19=geometric
',
926 $ 'alpha, beta=0,1
', / ' 17=arithm.
alpha&beta
',
927 $ ' 20=arithmetic
alpha, beta=0,1
', / ' 18=clustered
',
928 $ 'alpha, beta=0,1 21=random
alpha, beta=0,1
',
929 $ / ' large & small matrices:
', / ' 22=(large, small)
',
930 $ '23=(small,large) 24=(small,small) 25=(large,large)
',
931 $ / ' 26=random o(1) matrices.
' )
933 9993 FORMAT( / ' tests performed:
',
934 $ / ' 1 =
max | ( b a - a b )
''*l | / const.,
',
935 $ / ' 2 = | |vr(i)| - 1 | / ulp,
',
936 $ / ' 3 =
max | ( b a - a b )*r | / const.
',
937 $ / ' 4 = | |vl(i)| - 1 | / ulp,
',
938 $ / ' 5 = 0
if w same no matter
if r or l computed,
',
939 $ / ' 6 = 0
if l same no matter
if l computed,
',
940 $ / ' 7 = 0
if r same no matter
if r computed,
', / 1X )
941 9992 FORMAT( ' matrix order=
', I5, ', type=
', I2, ',
seed=
',
942 $ 4( I4, ',
' ), ' result
', I2, ' is
', 0P, F8.2 )
943 9991 FORMAT( ' matrix order=
', I5, ', type=
', I2, ',
seed=
',
944 $ 4( I4, ',
' ), ' result
', I2, ' is
', 1P, E10.3 )