348 SUBROUTINE zdrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
349 $ bi, z, q, alpha, beta, c, ldc, s, work, lwork,
350 $ rwork, iwork, liwork, bwork, info )
358 INTEGER info, lda, ldc, liwork, lwork, ncmax, nin,
360 DOUBLE PRECISION thresh
365 DOUBLE PRECISION rwork( * ), s( * )
366 COMPLEX*16 a( lda, * ), ai( lda, * ), alpha( * ),
367 $ b( lda, * ), beta( * ), bi( lda, * ),
368 $ c( ldc, * ), q( lda, * ), work( * ),
375 DOUBLE PRECISION zero, one, ten
376 parameter( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1 )
378 parameter( czero = ( 0.0d+0, 0.0d+0 ) )
383 INTEGER bdspac, i, ifunc, j, linfo, maxwrk, minwrk, mm,
384 $ mn2, nerrs, nptknt, ntest, ntestt, prtype, qba,
386 DOUBLE PRECISION abnrm, bignum, diftru, pltru, smlnum, temp1,
387 $ temp2, thrsh2, ulp, ulpinv, weight
391 DOUBLE PRECISION difest( 2 ), pl( 2 ), result( 10 )
405 INTEGER k, m, mplusn, n
408 common / mn / m, n, mplusn, k, fs
411 INTRINSIC abs, dble, dimag, max, sqrt
414 DOUBLE PRECISION abs1
417 abs1( x ) = abs( dble( x ) ) + abs( dimag( x ) )
424 IF( nsize.LT.0 )
THEN
426 ELSE IF( thresh.LT.zero )
THEN
428 ELSE IF( nin.LE.0 )
THEN
430 ELSE IF( nout.LE.0 )
THEN
432 ELSE IF( lda.LT.1 .OR. lda.LT.nsize )
THEN
434 ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 )
THEN
436 ELSE IF( liwork.LT.nsize+2 )
THEN
448 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
449 minwrk = 3*nsize*nsize / 2
453 maxwrk = nsize*( 1+
ilaenv( 1,
'ZGEQRF',
' ', nsize, 1, nsize,
455 maxwrk = max( maxwrk, nsize*( 1+
ilaenv( 1,
'ZUNGQR',
' ',
456 $ nsize, 1, nsize, -1 ) ) )
460 bdspac = 3*nsize*nsize / 2
461 maxwrk = max( maxwrk, nsize*nsize*
462 $ ( 1+
ilaenv( 1,
'ZGEBRD',
' ', nsize*nsize / 2,
463 $ nsize*nsize / 2, -1, -1 ) ) )
464 maxwrk = max( maxwrk, bdspac )
466 maxwrk = max( maxwrk, minwrk )
471 IF( lwork.LT.minwrk )
475 CALL
xerbla(
'ZDRGSX', -info )
483 smlnum =
dlamch(
'S' ) / ulp
484 bignum = one / smlnum
485 CALL
dlabad( smlnum, bignum )
507 DO 40 m = 1, nsize - 1
508 DO 30 n = 1, nsize - m
510 weight = one / weight
518 CALL
zlaset(
'Full', mplusn, mplusn, czero, czero, ai,
520 CALL
zlaset(
'Full', mplusn, mplusn, czero, czero, bi,
523 CALL
zlatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
524 $ lda, ai( 1, m+1 ), lda, bi, lda,
525 $ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
526 $ q, lda, z, lda, weight, qba, qbb )
533 IF( ifunc.EQ.0 )
THEN
535 ELSE IF( ifunc.EQ.1 )
THEN
537 ELSE IF( ifunc.EQ.2 )
THEN
539 ELSE IF( ifunc.EQ.3 )
THEN
543 CALL
zlacpy(
'Full', mplusn, mplusn, ai, lda, a, lda )
544 CALL
zlacpy(
'Full', mplusn, mplusn, bi, lda, b, lda )
547 $ lda, bi, lda, mm, alpha, beta, q, lda, z,
548 $ lda, pl, difest, work, lwork, rwork,
549 $ iwork, liwork, bwork, linfo )
551 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 )
THEN
553 WRITE( nout, fmt = 9999 )
'ZGGESX', linfo, mplusn,
561 CALL
zlacpy(
'Full', mplusn, mplusn, ai, lda, work,
563 CALL
zlacpy(
'Full', mplusn, mplusn, bi, lda,
564 $ work( mplusn*mplusn+1 ), mplusn )
565 abnrm =
zlange(
'Fro', mplusn, 2*mplusn, work, mplusn,
571 CALL
zget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
572 $ lda, work, rwork, result( 1 ) )
573 CALL
zget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
574 $ lda, work, rwork, result( 2 ) )
575 CALL
zget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
576 $ lda, work, rwork, result( 3 ) )
577 CALL
zget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
578 $ lda, work, rwork, result( 4 ) )
590 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
591 $ max( smlnum, abs1( alpha( j ) ),
592 $ abs1( ai( j, j ) ) )+
593 $ abs1( beta( j )-bi( j, j ) ) /
594 $ max( smlnum, abs1( beta( j ) ),
595 $ abs1( bi( j, j ) ) ) ) / ulp
596 IF( j.LT.mplusn )
THEN
597 IF( ai( j+1, j ).NE.zero )
THEN
603 IF( ai( j, j-1 ).NE.zero )
THEN
608 temp1 = max( temp1, temp2 )
610 WRITE( nout, fmt = 9997 )j, mplusn, prtype
619 IF( linfo.EQ.mplusn+3 )
THEN
621 ELSE IF( mm.NE.n )
THEN
630 mn2 = mm*( mplusn-mm )*2
631 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax )
THEN
636 CALL
zlakf2( mm, mplusn-mm, ai, lda,
637 $ ai( mm+1, mm+1 ), bi,
638 $ bi( mm+1, mm+1 ), c, ldc )
640 CALL
zgesvd(
'N',
'N', mn2, mn2, c, ldc, s, work,
641 $ 1, work( 2 ), 1, work( 3 ), lwork-2,
645 IF( difest( 2 ).EQ.zero )
THEN
646 IF( diftru.GT.abnrm*ulp )
647 $ result( 8 ) = ulpinv
648 ELSE IF( diftru.EQ.zero )
THEN
649 IF( difest( 2 ).GT.abnrm*ulp )
650 $ result( 8 ) = ulpinv
651 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
652 $ ( diftru*thrsh2.LT.difest( 2 ) ) )
THEN
653 result( 8 ) = max( diftru / difest( 2 ),
654 $ difest( 2 ) / diftru )
662 IF( linfo.EQ.( mplusn+2 ) )
THEN
663 IF( diftru.GT.abnrm*ulp )
664 $ result( 9 ) = ulpinv
665 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
666 $ result( 9 ) = ulpinv
667 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
668 $ result( 9 ) = ulpinv
672 ntestt = ntestt + ntest
677 IF( result( j ).GE.thresh )
THEN
682 IF( nerrs.EQ.0 )
THEN
683 WRITE( nout, fmt = 9996 )
'ZGX'
687 WRITE( nout, fmt = 9994 )
691 WRITE( nout, fmt = 9993 )
'unitary',
'''',
692 $
'transpose', (
'''', i = 1, 4 )
696 IF( result( j ).LT.10000.0d0 )
THEN
697 WRITE( nout, fmt = 9992 )mplusn, prtype,
698 $ weight, m, j, result( j )
700 WRITE( nout, fmt = 9991 )mplusn, prtype,
701 $ weight, m, j, result( j )
721 READ( nin, fmt = *,
END = 140 )mplusn
724 READ( nin, fmt = *,
END = 140 )n
726 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
729 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
731 READ( nin, fmt = * )pltru, diftru
738 CALL
zlacpy(
'Full', mplusn, mplusn, ai, lda, a, lda )
739 CALL
zlacpy(
'Full', mplusn, mplusn, bi, lda, b, lda )
744 CALL
zggesx(
'V',
'V',
'S',
zlctsx,
'B', mplusn, ai, lda, bi, lda,
745 $ mm, alpha, beta, q, lda, z, lda, pl, difest, work,
746 $ lwork, rwork, iwork, liwork, bwork, linfo )
748 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 )
THEN
750 WRITE( nout, fmt = 9998 )
'ZGGESX', linfo, mplusn, nptknt
757 CALL
zlacpy(
'Full', mplusn, mplusn, ai, lda, work, mplusn )
758 CALL
zlacpy(
'Full', mplusn, mplusn, bi, lda,
759 $ work( mplusn*mplusn+1 ), mplusn )
760 abnrm =
zlange(
'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
764 CALL
zget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
765 $ rwork, result( 1 ) )
766 CALL
zget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
767 $ rwork, result( 2 ) )
768 CALL
zget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
769 $ rwork, result( 3 ) )
770 CALL
zget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
771 $ rwork, result( 4 ) )
783 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
784 $ max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
785 $ abs1( beta( j )-bi( j, j ) ) /
786 $ max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
788 IF( j.LT.mplusn )
THEN
789 IF( ai( j+1, j ).NE.zero )
THEN
795 IF( ai( j, j-1 ).NE.zero )
THEN
800 temp1 = max( temp1, temp2 )
802 WRITE( nout, fmt = 9997 )j, mplusn, nptknt
811 IF( linfo.EQ.mplusn+3 )
812 $ result( 7 ) = ulpinv
818 IF( difest( 2 ).EQ.zero )
THEN
819 IF( diftru.GT.abnrm*ulp )
820 $ result( 8 ) = ulpinv
821 ELSE IF( diftru.EQ.zero )
THEN
822 IF( difest( 2 ).GT.abnrm*ulp )
823 $ result( 8 ) = ulpinv
824 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
825 $ ( diftru*thrsh2.LT.difest( 2 ) ) )
THEN
826 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
833 IF( linfo.EQ.( mplusn+2 ) )
THEN
834 IF( diftru.GT.abnrm*ulp )
835 $ result( 9 ) = ulpinv
836 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
837 $ result( 9 ) = ulpinv
838 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
839 $ result( 9 ) = ulpinv
846 IF( pl( 1 ).EQ.zero )
THEN
847 IF( pltru.GT.abnrm*ulp )
848 $ result( 10 ) = ulpinv
849 ELSE IF( pltru.EQ.zero )
THEN
850 IF( pl( 1 ).GT.abnrm*ulp )
851 $ result( 10 ) = ulpinv
852 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
853 $ ( pltru*thresh.LT.pl( 1 ) ) )
THEN
854 result( 10 ) = ulpinv
857 ntestt = ntestt + ntest
862 IF( result( j ).GE.thresh )
THEN
867 IF( nerrs.EQ.0 )
THEN
868 WRITE( nout, fmt = 9996 )
'ZGX'
872 WRITE( nout, fmt = 9995 )
876 WRITE( nout, fmt = 9993 )
'unitary',
'''',
'transpose',
881 IF( result( j ).LT.10000.0d0 )
THEN
882 WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
884 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
898 CALL
alasvm(
'ZGX', nout, nerrs, ntestt, 0 )
904 9999 format(
' ZDRGSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
905 $ i6,
', JTYPE=', i6,
')' )
907 9998 format(
' ZDRGSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
908 $ i6,
', Input Example #', i2,
')' )
910 9997 format(
' ZDRGSX: S not in Schur form at eigenvalue ', i6,
'.',
911 $ / 9x,
'N=', i6,
', JTYPE=', i6,
')' )
913 9996 format( / 1x, a3,
' -- Complex Expert Generalized Schur form',
914 $
' problem driver' )
916 9995 format(
'Input Example' )
918 9994 format(
' Matrix types: ', /
919 $
' 1: A is a block diagonal matrix of Jordan blocks ',
920 $
'and B is the identity ', /
' matrix, ',
921 $ /
' 2: A and B are upper triangular matrices, ',
922 $ /
' 3: A and B are as type 2, but each second diagonal ',
923 $
'block in A_11 and ', /
924 $
' each third diaongal block in A_22 are 2x2 blocks,',
925 $ /
' 4: A and B are block diagonal matrices, ',
926 $ /
' 5: (A,B) has potentially close or common ',
927 $
'eigenvalues.', / )
929 9993 format( /
' Tests performed: (S is Schur, T is triangular, ',
930 $
'Q and Z are ', a,
',', / 19x,
931 $
' a is alpha, b is beta, and ', a,
' means ', a,
'.)',
932 $ /
' 1 = | A - Q S Z', a,
933 $
' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
934 $
' | / ( |B| n ulp )', /
' 3 = | I - QQ', a,
935 $
' | / ( n ulp ) 4 = | I - ZZ', a,
936 $
' | / ( n ulp )', /
' 5 = 1/ULP if A is not in ',
937 $
'Schur form S', /
' 6 = difference between (alpha,beta)',
938 $
' and diagonals of (S,T)', /
939 $
' 7 = 1/ULP if SDIM is not the correct number of ',
940 $
'selected eigenvalues', /
941 $
' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
942 $
'DIFTRU/DIFEST > 10*THRESH',
943 $ /
' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
944 $
'when reordering fails', /
945 $
' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
946 $
'PLTRU/PLEST > THRESH', /
947 $
' ( Test 10 is only for input examples )', / )
948 9992 format(
' Matrix order=', i2,
', type=', i2,
', a=', d10.3,
949 $
', order(A_11)=', i2,
', result ', i2,
' is ', 0p, f8.2 )
950 9991 format(
' Matrix order=', i2,
', type=', i2,
', a=', d10.3,
951 $
', order(A_11)=', i2,
', result ', i2,
' is ', 0p, d10.3 )
952 9990 format(
' Input example #', i2,
', matrix order=', i4,
',',
953 $
' result ', i2,
' is', 0p, f8.2 )
954 9989 format(
' Input example #', i2,
', matrix order=', i4,
',',
955 $
' result ', i2,
' is', 1p, d10.3 )