404 SUBROUTINE sdrgev3( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
405 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
406 $ ALPHAR, ALPHAI, BETA, ALPHR1, ALPHI1, BETA1,
407 $ WORK, LWORK, RESULT, INFO )
414 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
420 INTEGER ISEED( 4 ), NN( * )
421 REAL A( LDA, * ), ALPHAI( * ), ALPHI1( * ),
422 $ alphar( * ), alphr1( * ), b( lda, * ),
423 $ beta( * ), beta1( * ), q( ldq, * ),
424 $ qe( ldqe, * ), result( * ), s( lda, * ),
425 $ t( lda, * ), work( * ), z( ldq, * )
432 PARAMETER ( ZERO = 0.0e+0, one = 1.0e+0 )
434 parameter( maxtyp = 26 )
438 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
439 $ MAXWRK, MINWRK, MTYPES, N, N1, NERRS, NMATS,
441 REAL SAFMAX, SAFMIN, ULP, ULPINV
444 INTEGER IASIGN( MAXTYP ), IBSIGN( MAXTYP ),
445 $ IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
446 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
447 $ kbmagn( maxtyp ), kbtype( maxtyp ),
448 $ kbzero( maxtyp ), kclass( maxtyp ),
449 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
455 EXTERNAL ILAENV, SLAMCH, SLARND
462 INTRINSIC abs, max, min, real, sign
465 DATA kclass / 15*1, 10*2, 1*3 /
466 DATA kz1 / 0, 1, 2, 1, 3, 3 /
467 DATA kz2 / 0, 0, 1, 2, 1, 1 /
468 DATA kadd / 0, 0, 0, 0, 3, 2 /
469 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
470 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
471 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
472 $ 1, 1, -4, 2, -4, 8*8, 0 /
473 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
475 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
477 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
479 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
481 DATA ktrian / 16*0, 10*1 /
482 DATA iasign / 6*0, 2, 0, 2*2, 2*0, 3*2, 0, 2, 3*0,
484 DATA ibsign / 7*0, 2, 2*0, 2*2, 2*0, 2, 0, 2, 9*0 /
495 nmax = max( nmax, nn( j ) )
500 IF( nsizes.LT.0 )
THEN
502 ELSE IF( badnn )
THEN
504 ELSE IF( ntypes.LT.0 )
THEN
506 ELSE IF( thresh.LT.zero )
THEN
508 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
510 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax )
THEN
512 ELSE IF( ldqe.LE.1 .OR. ldqe.LT.nmax )
THEN
524 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
525 minwrk = max( 1, 8*nmax, nmax*( nmax+1 ) )
526 maxwrk = 7*nmax + nmax*ilaenv( 1,
'SGEQRF',
' ', nmax, 1, nmax,
528 maxwrk = max( maxwrk, nmax*( nmax+1 ) )
532 IF( lwork.LT.minwrk )
536 CALL xerbla(
'SDRGEV3', -info )
542 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
545 safmin = slamch(
'Safe minimum' )
546 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
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 )
608 IF( mtypes.GT.maxtyp )
611 IF( kclass( jtype ).LT.3 )
THEN
615 IF( abs( katype( jtype ) ).EQ.3 )
THEN
616 in = 2*( ( n-1 ) / 2 ) + 1
618 $
CALL slaset(
'Full', n, n, zero, zero, a, lda )
622 CALL slatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
623 $ kz2( kazero( jtype ) ), iasign( jtype ),
624 $ rmagn( kamagn( jtype ) ), ulp,
625 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
627 iadd = kadd( kazero( jtype ) )
628 IF( iadd.GT.0 .AND. iadd.LE.n )
629 $ a( iadd, iadd ) = one
633 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
634 in = 2*( ( n-1 ) / 2 ) + 1
636 $
CALL slaset(
'Full', n, n, zero, zero, b, lda )
640 CALL slatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
641 $ kz2( kbzero( jtype ) ), ibsign( jtype ),
642 $ rmagn( kbmagn( jtype ) ), one,
643 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
645 iadd = kadd( kbzero( jtype ) )
646 IF( iadd.NE.0 .AND. iadd.LE.n )
647 $ b( iadd, iadd ) = one
649 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
658 q( jr, jc ) = slarnd( 3, iseed )
659 z( jr, jc ) = slarnd( 3, iseed )
661 CALL slarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
663 work( 2*n+jc ) = sign( one, q( jc, jc ) )
665 CALL slarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
667 work( 3*n+jc ) = sign( one, z( jc, jc ) )
672 work( 3*n ) = sign( one, slarnd( 2, iseed ) )
675 work( 4*n ) = sign( one, slarnd( 2, iseed ) )
681 a( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
683 b( jr, jc ) = work( 2*n+jr )*work( 3*n+jc )*
687 CALL sorm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
688 $ lda, work( 2*n+1 ), ierr )
691 CALL sorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
692 $ a, lda, work( 2*n+1 ), ierr )
695 CALL sorm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
696 $ lda, work( 2*n+1 ), ierr )
699 CALL sorm2r(
'R',
'T', n, n, n-1, z, ldq, work( n+1 ),
700 $ b, lda, work( 2*n+1 ), ierr )
710 a( jr, jc ) = rmagn( kamagn( jtype ) )*
712 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
721 WRITE( nounit, fmt = 9999 )
'Generator', ierr, n, jtype,
743 CALL slacpy(
' ', n, n, a, lda, s, lda )
744 CALL slacpy(
' ', n, n, b, lda, t, lda )
745 CALL sggev3(
'V',
'V', n, s, lda, t, lda, alphar, alphai,
746 $ beta, q, ldq, z, ldq, work, lwork, ierr )
747 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
749 WRITE( nounit, fmt = 9999 )
'SGGEV31', ierr, n, jtype,
757 CALL sget52( .true., n, a, lda, b, lda, q, ldq, alphar,
758 $ alphai, beta, work, result( 1 ) )
759 IF( result( 2 ).GT.thresh )
THEN
760 WRITE( nounit, fmt = 9998 )
'Left',
'SGGEV31',
761 $ result( 2 ), n, jtype, ioldsd
766 CALL sget52( .false., n, a, lda, b, lda, z, ldq, alphar,
767 $ alphai, beta, work, result( 3 ) )
768 IF( result( 4 ).GT.thresh )
THEN
769 WRITE( nounit, fmt = 9998 )
'Right',
'SGGEV31',
770 $ result( 4 ), n, jtype, ioldsd
775 CALL slacpy(
' ', n, n, a, lda, s, lda )
776 CALL slacpy(
' ', n, n, b, lda, t, lda )
777 CALL sggev3(
'N',
'N', n, s, lda, t, lda, alphr1, alphi1,
778 $ beta1, q, ldq, z, ldq, work, lwork, ierr )
779 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
781 WRITE( nounit, fmt = 9999 )
'SGGEV32', ierr, n, jtype,
788 IF( alphar( j ).NE.alphr1( j ) .OR.
789 $ beta( j ).NE. beta1( j ) )
THEN
797 CALL slacpy(
' ', n, n, a, lda, s, lda )
798 CALL slacpy(
' ', n, n, b, lda, t, lda )
799 CALL sggev3(
'V',
'N', n, s, lda, t, lda, alphr1, alphi1,
800 $ beta1, qe, ldqe, z, ldq, work, lwork, ierr )
801 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
803 WRITE( nounit, fmt = 9999 )
'SGGEV33', ierr, n, jtype,
810 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
811 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
812 $ result( 6 ) = ulpinv
817 IF( q( j, jc ).NE.qe( j, jc ) )
818 $ result( 6 ) = ulpinv
825 CALL slacpy(
' ', n, n, a, lda, s, lda )
826 CALL slacpy(
' ', n, n, b, lda, t, lda )
827 CALL sggev3(
'N',
'V', n, s, lda, t, lda, alphr1, alphi1,
828 $ beta1, q, ldq, qe, ldqe, work, lwork, ierr )
829 IF( ierr.NE.0 .AND. ierr.NE.n+1 )
THEN
831 WRITE( nounit, fmt = 9999 )
'SGGEV34', ierr, n, jtype,
838 IF( alphar( j ).NE.alphr1( j ) .OR. alphai( j ).NE.
839 $ alphi1( j ) .OR. beta( j ).NE.beta1( j ) )
840 $ result( 7 ) = ulpinv
845 IF( z( j, jc ).NE.qe( j, jc ) )
846 $ result( 7 ) = ulpinv
859 IF( result( jr ).GE.thresh )
THEN
864 IF( nerrs.EQ.0 )
THEN
865 WRITE( nounit, fmt = 9997 )
'SGV'
869 WRITE( nounit, fmt = 9996 )
870 WRITE( nounit, fmt = 9995 )
871 WRITE( nounit, fmt = 9994 )
'Orthogonal'
875 WRITE( nounit, fmt = 9993 )
879 IF( result( jr ).LT.10000.0 )
THEN
880 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
883 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
894 CALL alasvm(
'SGV', nounit, nerrs, ntestt, 0 )
900 9999
FORMAT(
' SDRGEV3: ', a,
' returned INFO=', i6,
'.', / 3x,
'N=',
901 $ i6,
', JTYPE=', i6,
', ISEED=(', 4( i4,
',' ), i5,
')' )
903 9998
FORMAT(
' SDRGEV3: ', a,
' Eigenvectors from ', a,
904 $
' incorrectly normalized.', /
' Bits of error=', 0p, g10.3,
905 $
',', 3x,
'N=', i4,
', JTYPE=', i3,
', ISEED=(',
906 $ 4( i4,
',' ), i5,
')' )
908 9997
FORMAT( / 1x, a3,
' -- Real Generalized eigenvalue problem driver'
911 9996
FORMAT(
' Matrix types (see SDRGEV3 for details): ' )
913 9995
FORMAT(
' Special Matrices:', 23x,
914 $
'(J''=transposed Jordan block)',
915 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
916 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
917 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
918 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
919 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
920 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
921 9994
FORMAT(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
922 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
923 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
924 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
925 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
926 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
927 $
'23=(small,large) 24=(small,small) 25=(large,large)',
928 $ /
' 26=random O(1) matrices.' )
930 9993
FORMAT( /
' Tests performed: ',
931 $ /
' 1 = max | ( b A - a B )''*l | / const.,',
932 $ /
' 2 = | |VR(i)| - 1 | / ulp,',
933 $ /
' 3 = max | ( b A - a B )*r | / const.',
934 $ /
' 4 = | |VL(i)| - 1 | / ulp,',
935 $ /
' 5 = 0 if W same no matter if r or l computed,',
936 $ /
' 6 = 0 if l same no matter if l computed,',
937 $ /
' 7 = 0 if r same no matter if r computed,', / 1x )
938 9992
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
939 $ 4( i4,
',' ),
' result ', i2,
' is', 0p, f8.2 )
940 9991
FORMAT(
' Matrix order=', i5,
', type=', i2,
', seed=',
941 $ 4( i4,
',' ),
' result ', i2,
' is', 1p, e10.3 )