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 )
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 alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine clatm4(ITYPE, N, NZ1, NZ2, RSIGN, AMAGN, RCOND, TRIANG, IDIST, ISEED, A, LDA)
CLATM4
subroutine xerbla(SRNAME, INFO)
XERBLA
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
complex function clarnd(IDIST, ISEED)
CLARND
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 clarfg(N, ALPHA, X, INCX, TAU)
CLARFG generates an elementary reflector (Householder matrix).
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
subroutine cget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA, WORK, RWORK, RESULT)
CGET52
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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 ...
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...