406 SUBROUTINE sdrgev( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
407 $ nounit, a, lda, b, s, t, q, ldq, z, qe, ldqe,
408 $ alphar, alphai, beta, alphr1, alphi1, beta1,
409 $ work, lwork, result, info )
417 INTEGER info, lda, ldq, ldqe, lwork, nounit, nsizes,
423 INTEGER iseed( 4 ), nn( * )
424 REAL a( lda, * ), alphai( * ), alphi1( * ),
425 $ alphar( * ), alphr1( * ), b( lda, * ),
426 $ beta( * ), beta1( * ), q( ldq, * ),
427 $ qe( ldqe, * ), result( * ), s( lda, * ),
428 $ t( lda, * ), work( * ), z( ldq, * )
435 parameter( zero = 0.0e+0, one = 1.0e+0 )
437 parameter( maxtyp = 26 )
441 INTEGER i, iadd, ierr, in, j, jc, jr, jsize, jtype,
442 $ maxwrk, minwrk, mtypes, n, n1, nerrs, nmats,
444 REAL safmax, safmin, ulp, ulpinv
447 INTEGER iasign( maxtyp ), ibsign( maxtyp ),
448 $ ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
449 $ katype( maxtyp ), kazero( maxtyp ),
450 $ kbmagn( maxtyp ), kbtype( maxtyp ),
451 $ kbzero( maxtyp ), kclass( maxtyp ),
452 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
465 INTRINSIC abs, max, min,
REAL, sign
468 DATA kclass / 15*1, 10*2, 1*3 /
469 DATA kz1 / 0, 1, 2, 1, 3, 3 /
470 DATA kz2 / 0, 0, 1, 2, 1, 1 /
471 DATA kadd / 0, 0, 0, 0, 3, 2 /
472 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
473 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
474 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
475 $ 1, 1, -4, 2, -4, 8*8, 0 /
476 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
478 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
480 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
482 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
484 DATA ktrian / 16*0, 10*1 /
485 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
487 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
498 nmax = max( nmax, nn( j ) )
503 IF( nsizes.LT.0 )
THEN
505 ELSE IF( badnn )
THEN
507 ELSE IF( ntypes.LT.0 )
THEN
509 ELSE IF( thresh.LT.zero )
THEN
511 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
513 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax )
THEN
515 ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax )
THEN
527 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
528 minwrk = max( 1, 8*nmax, nmax*( nmax+1 ) )
529 maxwrk = 7*nmax + nmax*
ilaenv( 1,
'SGEQRF',
' ', nmax, 1, nmax,
531 maxwrk = max( maxwrk, nmax*( nmax+1 ) )
535 IF( lwork.LT.minwrk )
539 CALL
xerbla(
'SDRGEV', -info )
545 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
548 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 )
612 IF( mtypes.GT.maxtyp )
615 IF( kclass( jtype ).LT.3 )
THEN
619 IF( abs( katype( jtype ) ).EQ.3 )
THEN
620 in = 2*( ( n-1 ) / 2 ) + 1
622 $ CALL
slaset(
'Full', n, n, zero, zero, a, lda )
626 CALL
slatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
627 $ kz2( kazero( jtype ) ), iasign( jtype ),
628 $ rmagn( kamagn( jtype ) ), ulp,
629 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
631 iadd = kadd( kazero( jtype ) )
632 IF( iadd.GT.0 .AND. iadd.LE.n )
633 $ a( iadd, iadd ) = one
637 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
638 in = 2*( ( n-1 ) / 2 ) + 1
640 $ CALL
slaset(
'Full', n, n, zero, zero, b, lda )
644 CALL
slatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
645 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
646 $ rmagn( kbmagn( jtype ) ), one,
647 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
649 iadd = kadd( kbzero( jtype ) )
650 IF( iadd.NE.0 .AND. iadd.LE.n )
651 $ b( iadd, iadd ) = one
653 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
662 q( jr, jc ) =
slarnd( 3, iseed )
663 z( jr, jc ) =
slarnd( 3, iseed )
665 CALL
slarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
667 work( 2*n+jc ) = sign( one, q( jc, jc ) )
669 CALL
slarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
671 work( 3*n+jc ) = sign( one, z( jc, jc ) )
676 work( 3*n ) = sign( one,
slarnd( 2, iseed ) )
679 work( 4*n ) = sign( one,
slarnd( 2, iseed ) )
685 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
687 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
691 CALL
sorm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
692 $ lda, work( 2*n+1 ), ierr )
695 CALL
sorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
696 $ a, lda, work( 2*n+1 ), ierr )
699 CALL
sorm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
700 $ lda, work( 2*n+1 ), ierr )
703 CALL
sorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
704 $ b, lda, work( 2*n+1 ), ierr )
714 a( jr, jc ) = rmagn( kamagn( jtype ) )*
716 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
725 WRITE( nounit, fmt = 9999 )
'Generator', ierr, n, jtype,
739 CALL
slacpy(
' ', n, n, a, lda, s, lda )
740 CALL
slacpy(
' ', n, n, b, lda, t, lda )
741 CALL
sggev(
'V',
'V', n, s, lda, t, lda, alphar, alphai,
742 $ beta, q, ldq, z, ldq, work, lwork, ierr )
743 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
745 WRITE( nounit, fmt = 9999 )
'SGGEV1', ierr, n, jtype,
753 CALL
sget52( .true., n, a, lda, b, lda, q, ldq, alphar,
754 $ alphai, beta, work, result( 1 ) )
755 IF( result( 2 ).GT.thresh )
THEN
756 WRITE( nounit, fmt = 9998 )
'Left',
'SGGEV1',
757 $ result( 2 ), n, jtype, ioldsd
762 CALL
sget52( .false., n, a, lda, b, lda, z, ldq, alphar,
763 $ alphai, beta, work, result( 3 ) )
764 IF( result( 4 ).GT.thresh )
THEN
765 WRITE( nounit, fmt = 9998 )
'Right',
'SGGEV1',
766 $ result( 4 ), n, jtype, ioldsd
771 CALL
slacpy(
' ', n, n, a, lda, s, lda )
772 CALL
slacpy(
' ', n, n, b, lda, t, lda )
773 CALL
sggev(
'N',
'N', n, s, lda, t, lda, alphr1, alphi1,
774 $ beta1, q, ldq, z, ldq, work, lwork, ierr )
775 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
777 WRITE( nounit, fmt = 9999 )
'SGGEV2', ierr, n, jtype,
784 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
785 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
786 $ result( 5 ) = ulpinv
792 CALL
slacpy(
' ', n, n, a, lda, s, lda )
793 CALL
slacpy(
' ', n, n, b, lda, t, lda )
794 CALL
sggev(
'V',
'N', n, s, lda, t, lda, alphr1, alphi1,
795 $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
796 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
798 WRITE( nounit, fmt = 9999 )
'SGGEV3', ierr, n, jtype,
805 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
806 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
807 $ result( 6 ) = ulpinv
812 IF( q( j, jc ).NE.qe( j, jc ) )
813 $ result( 6 ) = ulpinv
820 CALL
slacpy(
' ', n, n, a, lda, s, lda )
821 CALL
slacpy(
' ', n, n, b, lda, t, lda )
822 CALL
sggev(
'N',
'V', n, s, lda, t, lda, alphr1, alphi1,
823 $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
824 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
826 WRITE( nounit, fmt = 9999 )
'SGGEV4', ierr, n, jtype,
833 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
834 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
835 $ result( 7 ) = ulpinv
840 IF( z( j, jc ).NE.qe( j, jc ) )
841 $ result( 7 ) = ulpinv
854 IF( result( jr ).GE.thresh )
THEN
859 IF( nerrs.EQ.0 )
THEN
860 WRITE( nounit, fmt = 9997 )
'SGV'
864 WRITE( nounit, fmt = 9996 )
865 WRITE( nounit, fmt = 9995 )
866 WRITE( nounit, fmt = 9994 )
'Orthogonal'
870 WRITE( nounit, fmt = 9993 )
874 IF( result( jr ).LT.10000.0 )
THEN
875 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
878 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
889 CALL
alasvm(
'SGV', nounit, nerrs, ntestt, 0 )
895 9999 format(
' SDRGEV: ', a,
' returned INFO=', i6,
'.', / 3x,
'N=',
896 $ i6,
', JTYPE=', i6,
', ISEED=(', 4( i4,
',' ), i5,
')' )
898 9998 format(
' SDRGEV: ', a,
' Eigenvectors from ', a,
' incorrectly ',
899 $
'normalized.', /
' Bits of error=', 0p, g10.3,
',', 3x,
900 $
'N=', i4,
', JTYPE=', i3,
', ISEED=(', 4( i4,
',' ), i5,
903 9997 format( / 1x, a3,
' -- Real Generalized eigenvalue problem driver'
906 9996 format(
' Matrix types (see SDRGEV for details): ' )
908 9995 format(
' Special Matrices:', 23x,
909 $
'(J''=transposed Jordan block)',
910 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
911 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
912 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
913 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
914 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
915 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
916 9994 format(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
917 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
918 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
919 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
920 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
921 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
922 $
'23=(small,large) 24=(small,small) 25=(large,large)',
923 $ /
' 26=random O(1) matrices.' )
925 9993 format( /
' Tests performed: ',
926 $ /
' 1 = max | ( b A - a B )''*l | / const.,',
927 $ /
' 2 = | |VR(i)| - 1 | / ulp,',
928 $ /
' 3 = max | ( b A - a B )*r | / const.',
929 $ /
' 4 = | |VL(i)| - 1 | / ulp,',
930 $ /
' 5 = 0 if W same no matter if r or l computed,',
931 $ /
' 6 = 0 if l same no matter if l computed,',
932 $ /
' 7 = 0 if r same no matter if r computed,', / 1x )
933 9992 format(
' Matrix order=', i5,
', type=', i2,
', seed=',
934 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
935 9991 format(
' Matrix order=', i5,
', type=', i2,
', seed=',
936 $ 4( i4,
',' ),
' result ', i2,
' is', 1p, e10.3 )