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 )
396 DOUBLE PRECISION DLAMCH, ZLANGE
397 EXTERNAL zlctsx, ilaenv, dlamch, zlange
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 )
546 CALL zggesx(
'V',
'V',
'S', zlctsx, sense, mplusn, ai,
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 )
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine zlatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
ZLATM5
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
ZLAKF2
subroutine zgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
ZGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine zdrgsx(NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
ZDRGSX
subroutine zggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
ZGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
subroutine zget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
ZGET51