380 SUBROUTINE cdrges( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
381 $ nounit, a, lda, b, s, t, q, ldq, z, alpha,
382 $ beta, work, lwork, rwork, result, bwork, info )
390 INTEGER info, lda, ldq, lwork, nounit, nsizes, ntypes
394 LOGICAL bwork( * ), dotype( * )
395 INTEGER iseed( 4 ), nn( * )
396 REAL result( 13 ), rwork( * )
397 COMPLEX a( lda, * ), alpha( * ), b( lda, * ),
398 $ beta( * ), q( ldq, * ), s( lda, * ),
399 $ t( lda, * ), work( * ), z( ldq, * )
406 parameter( zero = 0.0e+0, one = 1.0e+0 )
408 parameter( czero = ( 0.0e+0, 0.0e+0 ),
409 $ cone = ( 1.0e+0, 0.0e+0 ) )
411 parameter( maxtyp = 26 )
414 LOGICAL badnn, ilabad
416 INTEGER i, iadd, iinfo, in, isort, j, jc, jr, jsize,
417 $ jtype, knteig, maxwrk, minwrk, mtypes, n, n1,
418 $ nb, nerrs, nmats, nmax, ntest, ntestt, rsub,
420 REAL safmax, safmin, temp1, temp2, ulp, ulpinv
424 LOGICAL lasign( maxtyp ), lbsign( maxtyp )
425 INTEGER ioldsd( 4 ), kadd( 6 ), kamagn( maxtyp ),
426 $ katype( maxtyp ), kazero( maxtyp ),
427 $ kbmagn( maxtyp ), kbtype( maxtyp ),
428 $ kbzero( maxtyp ), kclass( maxtyp ),
429 $ ktrian( maxtyp ), kz1( 6 ), kz2( 6 )
444 INTRINSIC abs, aimag, conjg, max, min,
REAL, sign
450 abs1( x ) = abs(
REAL( X ) ) + abs( aimag( x ) )
453 DATA kclass / 15*1, 10*2, 1*3 /
454 DATA kz1 / 0, 1, 2, 1, 3, 3 /
455 DATA kz2 / 0, 0, 1, 2, 1, 1 /
456 DATA kadd / 0, 0, 0, 0, 3, 2 /
457 DATA katype / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
458 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
459 DATA kbtype / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
460 $ 1, 1, -4, 2, -4, 8*8, 0 /
461 DATA kazero / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
463 DATA kbzero / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
465 DATA kamagn / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
467 DATA kbmagn / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
469 DATA ktrian / 16*0, 10*1 /
470 DATA lasign / 6*.false., .true., .false., 2*.true.,
471 $ 2*.false., 3*.true., .false., .true.,
472 $ 3*.false., 5*.true., .false. /
473 DATA lbsign / 7*.false., .true., 2*.false.,
474 $ 2*.true., 2*.false., .true., .false., .true.,
486 nmax = max( nmax, nn( j ) )
491 IF( nsizes.LT.0 )
THEN
493 ELSE IF( badnn )
THEN
495 ELSE IF( ntypes.LT.0 )
THEN
497 ELSE IF( thresh.LT.zero )
THEN
499 ELSE IF( lda.LE.1 .OR. lda.LT.nmax )
THEN
501 ELSE IF( ldq.LE.1 .OR. ldq.LT.nmax )
THEN
513 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
515 nb = max( 1,
ilaenv( 1,
'CGEQRF',
' ', nmax, nmax, -1, -1 ),
516 $
ilaenv( 1,
'CUNMQR',
'LC', nmax, nmax, nmax, -1 ),
517 $
ilaenv( 1,
'CUNGQR',
' ', nmax, nmax, nmax, -1 ) )
518 maxwrk = max( nmax+nmax*nb, 3*nmax*nmax )
522 IF( lwork.LT.minwrk )
526 CALL
xerbla(
'CDRGES', -info )
532 IF( nsizes.EQ.0 .OR. ntypes.EQ.0 )
535 ulp =
slamch(
'Precision' )
536 safmin =
slamch(
'Safe minimum' )
537 safmin = safmin / ulp
538 safmax = one / safmin
539 CALL
slabad( safmin, safmax )
553 DO 190 jsize = 1, nsizes
556 rmagn( 2 ) = safmax*ulp /
REAL( n1 )
557 rmagn( 3 ) = safmin*ulpinv*
REAL( n1 )
559 IF( nsizes.NE.1 )
THEN
560 mtypes = min( maxtyp, ntypes )
562 mtypes = min( maxtyp+1, ntypes )
567 DO 180 jtype = 1, mtypes
568 IF( .NOT.dotype( jtype ) )
576 ioldsd( j ) = iseed( j )
606 IF( mtypes.GT.maxtyp )
609 IF( kclass( jtype ).LT.3 )
THEN
613 IF( abs( katype( jtype ) ).EQ.3 )
THEN
614 in = 2*( ( n-1 ) / 2 ) + 1
616 $ CALL
claset(
'Full', n, n, czero, czero, a, lda )
620 CALL
clatm4( katype( jtype ), in, kz1( kazero( jtype ) ),
621 $ kz2( kazero( jtype ) ), lasign( jtype ),
622 $ rmagn( kamagn( jtype ) ), ulp,
623 $ rmagn( ktrian( jtype )*kamagn( jtype ) ), 2,
625 iadd = kadd( kazero( jtype ) )
626 IF( iadd.GT.0 .AND. iadd.LE.n )
627 $ a( iadd, iadd ) = rmagn( kamagn( jtype ) )
631 IF( abs( kbtype( jtype ) ).EQ.3 )
THEN
632 in = 2*( ( n-1 ) / 2 ) + 1
634 $ CALL
claset(
'Full', n, n, czero, czero, b, lda )
638 CALL
clatm4( kbtype( jtype ), in, kz1( kbzero( jtype ) ),
639 $ kz2( kbzero( jtype ) ), lbsign( jtype ),
640 $ rmagn( kbmagn( jtype ) ), one,
641 $ rmagn( ktrian( jtype )*kbmagn( jtype ) ), 2,
643 iadd = kadd( kbzero( jtype ) )
644 IF( iadd.NE.0 .AND. iadd.LE.n )
645 $ b( iadd, iadd ) = rmagn( kbmagn( jtype ) )
647 IF( kclass( jtype ).EQ.2 .AND. n.GT.0 )
THEN
656 q( jr, jc ) =
clarnd( 3, iseed )
657 z( jr, jc ) =
clarnd( 3, iseed )
659 CALL
clarfg( n+1-jc, q( jc, jc ), q( jc+1, jc ), 1,
661 work( 2*n+jc ) = sign( one,
REAL( Q( JC, JC ) ) )
663 CALL
clarfg( n+1-jc, z( jc, jc ), z( jc+1, jc ), 1,
665 work( 3*n+jc ) = sign( one,
REAL( Z( JC, JC ) ) )
668 ctemp =
clarnd( 3, iseed )
671 work( 3*n ) = ctemp / abs( ctemp )
672 ctemp =
clarnd( 3, iseed )
675 work( 4*n ) = ctemp / abs( ctemp )
681 a( jr, jc ) = work( 2*n+jr )*
682 $ conjg( work( 3*n+jc ) )*
684 b( jr, jc ) = work( 2*n+jr )*
685 $ conjg( work( 3*n+jc ) )*
689 CALL
cunm2r(
'L',
'N', n, n, n-1, q, ldq, work, a,
690 $ lda, work( 2*n+1 ), iinfo )
693 CALL
cunm2r(
'R',
'C', n, n, n-1, z, ldq, work( n+1 ),
694 $ a, lda, work( 2*n+1 ), iinfo )
697 CALL
cunm2r(
'L',
'N', n, n, n-1, q, ldq, work, b,
698 $ lda, work( 2*n+1 ), iinfo )
701 CALL
cunm2r(
'R',
'C', n, n, n-1, z, ldq, work( n+1 ),
702 $ b, lda, work( 2*n+1 ), iinfo )
712 a( jr, jc ) = rmagn( kamagn( jtype ) )*
714 b( jr, jc ) = rmagn( kbmagn( jtype ) )*
722 IF( iinfo.NE.0 )
THEN
723 WRITE( nounit, fmt = 9999 )
'Generator', iinfo, n, jtype,
738 IF( isort.EQ.0 )
THEN
748 CALL
clacpy(
'Full', n, n, a, lda, s, lda )
749 CALL
clacpy(
'Full', n, n, b, lda, t, lda )
750 ntest = 1 + rsub + isort
751 result( 1+rsub+isort ) = ulpinv
752 CALL
cgges(
'V',
'V', sort,
clctes, n, s, lda, t, lda,
753 $ sdim, alpha, beta, q, ldq, z, ldq, work,
754 $ lwork, rwork, bwork, iinfo )
755 IF( iinfo.NE.0 .AND. iinfo.NE.n+2 )
THEN
756 result( 1+rsub+isort ) = ulpinv
757 WRITE( nounit, fmt = 9999 )
'CGGES', iinfo, n, jtype,
767 IF( isort.EQ.0 )
THEN
768 CALL
cget51( 1, n, a, lda, s, lda, q, ldq, z, ldq,
769 $ work, rwork, result( 1 ) )
770 CALL
cget51( 1, n, b, lda, t, lda, q, ldq, z, ldq,
771 $ work, rwork, result( 2 ) )
773 CALL
cget54( n, a, lda, b, lda, s, lda, t, lda, q,
774 $ ldq, z, ldq, work, result( 2+rsub ) )
777 CALL
cget51( 3, n, b, lda, t, lda, q, ldq, q, ldq, work,
778 $ rwork, result( 3+rsub ) )
779 CALL
cget51( 3, n, b, lda, t, lda, z, ldq, z, ldq, work,
780 $ rwork, result( 4+rsub ) )
791 temp2 = ( abs1( alpha( j )-s( j, j ) ) /
792 $ max( safmin, abs1( alpha( j ) ), abs1( s( j,
793 $ j ) ) )+abs1( beta( j )-t( j, j ) ) /
794 $ max( safmin, abs1( beta( j ) ), abs1( t( j,
798 IF( s( j+1, j ).NE.zero )
THEN
800 result( 5+rsub ) = ulpinv
804 IF( s( j, j-1 ).NE.zero )
THEN
806 result( 5+rsub ) = ulpinv
809 temp1 = max( temp1, temp2 )
811 WRITE( nounit, fmt = 9998 )j, n, jtype, ioldsd
814 result( 6+rsub ) = temp1
816 IF( isort.GE.1 )
THEN
824 IF(
clctes( alpha( i ), beta( i ) ) )
825 $ knteig = knteig + 1
828 $ result( 13 ) = ulpinv
837 ntestt = ntestt + ntest
842 IF( result( jr ).GE.thresh )
THEN
847 IF( nerrs.EQ.0 )
THEN
848 WRITE( nounit, fmt = 9997 )
'CGS'
852 WRITE( nounit, fmt = 9996 )
853 WRITE( nounit, fmt = 9995 )
854 WRITE( nounit, fmt = 9994 )
'Unitary'
858 WRITE( nounit, fmt = 9993 )
'unitary',
'''',
859 $
'transpose', (
'''', j = 1, 8 )
863 IF( result( jr ).LT.10000.0 )
THEN
864 WRITE( nounit, fmt = 9992 )n, jtype, ioldsd, jr,
867 WRITE( nounit, fmt = 9991 )n, jtype, ioldsd, jr,
878 CALL
alasvm(
'CGS', nounit, nerrs, ntestt, 0 )
884 9999 format(
' CDRGES: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
885 $ i6,
', JTYPE=', i6,
', ISEED=(', 4( i4,
',' ), i5,
')' )
887 9998 format(
' CDRGES: S not in Schur form at eigenvalue ', i6,
'.',
888 $ / 9x,
'N=', i6,
', JTYPE=', i6,
', ISEED=(', 3( i5,
',' ),
891 9997 format( / 1x, a3,
' -- Complex Generalized Schur from problem ',
894 9996 format(
' Matrix types (see CDRGES for details): ' )
896 9995 format(
' Special Matrices:', 23x,
897 $
'(J''=transposed Jordan block)',
898 $ /
' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
899 $
'6=(diag(J'',I), diag(I,J''))', /
' Diagonal Matrices: ( ',
900 $
'D=diag(0,1,2,...) )', /
' 7=(D,I) 9=(large*D, small*I',
901 $
') 11=(large*I, small*D) 13=(large*D, large*I)', /
902 $
' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
903 $
' 14=(small*D, small*I)', /
' 15=(D, reversed D)' )
904 9994 format(
' Matrices Rotated by Random ', a,
' Matrices U, V:',
905 $ /
' 16=Transposed Jordan Blocks 19=geometric ',
906 $
'alpha, beta=0,1', /
' 17=arithm. alpha&beta ',
907 $
' 20=arithmetic alpha, beta=0,1', /
' 18=clustered ',
908 $
'alpha, beta=0,1 21=random alpha, beta=0,1',
909 $ /
' Large & Small Matrices:', /
' 22=(large, small) ',
910 $
'23=(small,large) 24=(small,small) 25=(large,large)',
911 $ /
' 26=random O(1) matrices.' )
913 9993 format( /
' Tests performed: (S is Schur, T is triangular, ',
914 $
'Q and Z are ', a,
',', / 19x,
915 $
'l and r are the appropriate left and right', / 19x,
916 $
'eigenvectors, resp., a is alpha, b is beta, and', / 19x, a,
917 $
' means ', a,
'.)', /
' Without ordering: ',
918 $ /
' 1 = | A - Q S Z', a,
919 $
' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
920 $
' | / ( |B| n ulp )', /
' 3 = | I - QQ', a,
921 $
' | / ( n ulp ) 4 = | I - ZZ', a,
922 $
' | / ( n ulp )', /
' 5 = A is in Schur form S',
923 $ /
' 6 = difference between (alpha,beta)',
924 $
' and diagonals of (S,T)', /
' With ordering: ',
925 $ /
' 7 = | (A,B) - Q (S,T) Z', a,
' | / ( |(A,B)| n ulp )',
926 $ /
' 8 = | I - QQ', a,
927 $
' | / ( n ulp ) 9 = | I - ZZ', a,
928 $
' | / ( n ulp )', /
' 10 = A is in Schur form S',
929 $ /
' 11 = difference between (alpha,beta) and diagonals',
930 $
' of (S,T)', /
' 12 = SDIM is the correct number of ',
931 $
'selected eigenvalues', / )
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 )