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 )
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 )