348 SUBROUTINE cdrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
349 $ ai, bi, z, q, alpha, beta, c, ldc, s, work,
350 $ lwork, rwork, iwork, liwork, bwork, info )
358 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
365 REAL RWORK( * ), S( * )
366 COMPLEX A( lda, * ), AI( lda, * ), ALPHA( * ),
367 $ b( lda, * ), beta( * ), bi( lda, * ),
368 $ c( ldc, * ), q( lda, * ), work( * ),
376 parameter ( zero = 0.0e+0, one = 1.0e+0, ten = 1.0e+1 )
378 parameter ( czero = ( 0.0e+0, 0.0e+0 ) )
383 INTEGER BDSPAC, I, IFUNC, J, LINFO, MAXWRK, MINWRK, MM,
384 $ mn2, nerrs, nptknt, ntest, ntestt, prtype, qba,
386 REAL ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
387 $ temp2, thrsh2, ulp, ulpinv, weight
391 REAL DIFEST( 2 ), PL( 2 ), RESULT( 10 )
397 EXTERNAL clctsx, ilaenv, clange, slamch
405 INTEGER K, M, MPLUSN, N
408 COMMON / mn / m, n, mplusn, k, fs
411 INTRINSIC abs, aimag, max,
REAL, SQRT
417 abs1( x ) = abs(
REAL( X ) ) + abs( AIMAG( x ) )
423 IF( nsize.LT.0 )
THEN
425 ELSE IF( thresh.LT.zero )
THEN
427 ELSE IF( nin.LE.0 )
THEN
429 ELSE IF( nout.LE.0 )
THEN
431 ELSE IF( lda.LT.1 .OR. lda.LT.nsize )
THEN
433 ELSE IF( ldc.LT.1 .OR. ldc.LT.nsize*nsize / 2 )
THEN
435 ELSE IF( liwork.LT.nsize+2 )
THEN
447 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
448 minwrk = 3*nsize*nsize / 2
452 maxwrk = nsize*( 1+ilaenv( 1,
'CGEQRF',
' ', nsize, 1, nsize,
454 maxwrk = max( maxwrk, nsize*( 1+ilaenv( 1,
'CUNGQR',
' ',
455 $ nsize, 1, nsize, -1 ) ) )
459 bdspac = 3*nsize*nsize / 2
460 maxwrk = max( maxwrk, nsize*nsize*
461 $ ( 1+ilaenv( 1,
'CGEBRD',
' ', nsize*nsize / 2,
462 $ nsize*nsize / 2, -1, -1 ) ) )
463 maxwrk = max( maxwrk, bdspac )
465 maxwrk = max( maxwrk, minwrk )
470 IF( lwork.LT.minwrk )
474 CALL xerbla(
'CDRGSX', -info )
482 smlnum = slamch(
'S' ) / ulp
483 bignum = one / smlnum
484 CALL slabad( smlnum, bignum )
506 DO 40 m = 1, nsize - 1
507 DO 30 n = 1, nsize - m
509 weight = one / weight
517 CALL claset(
'Full', mplusn, mplusn, czero, czero, ai,
519 CALL claset(
'Full', mplusn, mplusn, czero, czero, bi,
522 CALL clatm5( prtype, m, n, ai, lda, ai( m+1, m+1 ),
523 $ lda, ai( 1, m+1 ), lda, bi, lda,
524 $ bi( m+1, m+1 ), lda, bi( 1, m+1 ), lda,
525 $ q, lda, z, lda, weight, qba, qbb )
532 IF( ifunc.EQ.0 )
THEN
534 ELSE IF( ifunc.EQ.1 )
THEN
536 ELSE IF( ifunc.EQ.2 )
THEN
538 ELSE IF( ifunc.EQ.3 )
THEN
542 CALL clacpy(
'Full', mplusn, mplusn, ai, lda, a, lda )
543 CALL clacpy(
'Full', mplusn, mplusn, bi, lda, b, lda )
545 CALL cggesx(
'V',
'V',
'S', clctsx, sense, mplusn, ai,
546 $ lda, bi, lda, mm, alpha, beta, q, lda, z,
547 $ lda, pl, difest, work, lwork, rwork,
548 $ iwork, liwork, bwork, linfo )
550 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 )
THEN
552 WRITE( nout, fmt = 9999 )
'CGGESX', linfo, mplusn,
560 CALL clacpy(
'Full', mplusn, mplusn, ai, lda, work,
562 CALL clacpy(
'Full', mplusn, mplusn, bi, lda,
563 $ work( mplusn*mplusn+1 ), mplusn )
564 abnrm = clange(
'Fro', mplusn, 2*mplusn, work, mplusn,
570 CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
571 $ lda, work, rwork, result( 1 ) )
572 CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
573 $ lda, work, rwork, result( 2 ) )
574 CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
575 $ lda, work, rwork, result( 3 ) )
576 CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
577 $ lda, work, rwork, result( 4 ) )
589 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
590 $ max( smlnum, abs1( alpha( j ) ),
591 $ abs1( ai( j, j ) ) )+
592 $ abs1( beta( j )-bi( j, j ) ) /
593 $ max( smlnum, abs1( beta( j ) ),
594 $ abs1( bi( j, j ) ) ) ) / ulp
595 IF( j.LT.mplusn )
THEN
596 IF( ai( j+1, j ).NE.zero )
THEN
602 IF( ai( j, j-1 ).NE.zero )
THEN
607 temp1 = max( temp1, temp2 )
609 WRITE( nout, fmt = 9997 )j, mplusn, prtype
618 IF( linfo.EQ.mplusn+3 )
THEN
620 ELSE IF( mm.NE.n )
THEN
629 mn2 = mm*( mplusn-mm )*2
630 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax )
THEN
635 CALL clakf2( mm, mplusn-mm, ai, lda,
636 $ ai( mm+1, mm+1 ), bi,
637 $ bi( mm+1, mm+1 ), c, ldc )
639 CALL cgesvd(
'N',
'N', mn2, mn2, c, ldc, s, work,
640 $ 1, work( 2 ), 1, work( 3 ), lwork-2,
644 IF( difest( 2 ).EQ.zero )
THEN
645 IF( diftru.GT.abnrm*ulp )
646 $ result( 8 ) = ulpinv
647 ELSE IF( diftru.EQ.zero )
THEN
648 IF( difest( 2 ).GT.abnrm*ulp )
649 $ result( 8 ) = ulpinv
650 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
651 $ ( diftru*thrsh2.LT.difest( 2 ) ) )
THEN
652 result( 8 ) = max( diftru / difest( 2 ),
653 $ difest( 2 ) / diftru )
661 IF( linfo.EQ.( mplusn+2 ) )
THEN
662 IF( diftru.GT.abnrm*ulp )
663 $ result( 9 ) = ulpinv
664 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
665 $ result( 9 ) = ulpinv
666 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
667 $ result( 9 ) = ulpinv
671 ntestt = ntestt + ntest
676 IF( result( j ).GE.thresh )
THEN
681 IF( nerrs.EQ.0 )
THEN
682 WRITE( nout, fmt = 9996 )
'CGX'
686 WRITE( nout, fmt = 9994 )
690 WRITE( nout, fmt = 9993 )
'unitary',
'''',
691 $
'transpose', (
'''', i = 1, 4 )
695 IF( result( j ).LT.10000.0 )
THEN
696 WRITE( nout, fmt = 9992 )mplusn, prtype,
697 $ weight, m, j, result( j )
699 WRITE( nout, fmt = 9991 )mplusn, prtype,
700 $ weight, m, j, result( j )
720 READ( nin, fmt = *, end = 140 )mplusn
723 READ( nin, fmt = *, end = 140 )n
725 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
728 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
730 READ( nin, fmt = * )pltru, diftru
737 CALL clacpy(
'Full', mplusn, mplusn, ai, lda, a, lda )
738 CALL clacpy(
'Full', mplusn, mplusn, bi, lda, b, lda )
743 CALL cggesx(
'V',
'V',
'S', clctsx,
'B', mplusn, ai, lda, bi, lda,
744 $ mm, alpha, beta, q, lda, z, lda, pl, difest, work,
745 $ lwork, rwork, iwork, liwork, bwork, linfo )
747 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 )
THEN
749 WRITE( nout, fmt = 9998 )
'CGGESX', linfo, mplusn, nptknt
756 CALL clacpy(
'Full', mplusn, mplusn, ai, lda, work, mplusn )
757 CALL clacpy(
'Full', mplusn, mplusn, bi, lda,
758 $ work( mplusn*mplusn+1 ), mplusn )
759 abnrm = clange(
'Fro', mplusn, 2*mplusn, work, mplusn, rwork )
763 CALL cget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
764 $ rwork, result( 1 ) )
765 CALL cget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
766 $ rwork, result( 2 ) )
767 CALL cget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
768 $ rwork, result( 3 ) )
769 CALL cget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
770 $ rwork, result( 4 ) )
782 temp2 = ( abs1( alpha( j )-ai( j, j ) ) /
783 $ max( smlnum, abs1( alpha( j ) ), abs1( ai( j, j ) ) )+
784 $ abs1( beta( j )-bi( j, j ) ) /
785 $ max( smlnum, abs1( beta( j ) ), abs1( bi( j, j ) ) ) )
787 IF( j.LT.mplusn )
THEN
788 IF( ai( j+1, j ).NE.zero )
THEN
794 IF( ai( j, j-1 ).NE.zero )
THEN
799 temp1 = max( temp1, temp2 )
801 WRITE( nout, fmt = 9997 )j, mplusn, nptknt
810 IF( linfo.EQ.mplusn+3 )
811 $ result( 7 ) = ulpinv
817 IF( difest( 2 ).EQ.zero )
THEN
818 IF( diftru.GT.abnrm*ulp )
819 $ result( 8 ) = ulpinv
820 ELSE IF( diftru.EQ.zero )
THEN
821 IF( difest( 2 ).GT.abnrm*ulp )
822 $ result( 8 ) = ulpinv
823 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
824 $ ( diftru*thrsh2.LT.difest( 2 ) ) )
THEN
825 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
832 IF( linfo.EQ.( mplusn+2 ) )
THEN
833 IF( diftru.GT.abnrm*ulp )
834 $ result( 9 ) = ulpinv
835 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
836 $ result( 9 ) = ulpinv
837 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
838 $ result( 9 ) = ulpinv
845 IF( pl( 1 ).EQ.zero )
THEN
846 IF( pltru.GT.abnrm*ulp )
847 $ result( 10 ) = ulpinv
848 ELSE IF( pltru.EQ.zero )
THEN
849 IF( pl( 1 ).GT.abnrm*ulp )
850 $ result( 10 ) = ulpinv
851 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
852 $ ( pltru*thresh.LT.pl( 1 ) ) )
THEN
853 result( 10 ) = ulpinv
856 ntestt = ntestt + ntest
861 IF( result( j ).GE.thresh )
THEN
866 IF( nerrs.EQ.0 )
THEN
867 WRITE( nout, fmt = 9996 )
'CGX'
871 WRITE( nout, fmt = 9995 )
875 WRITE( nout, fmt = 9993 )
'unitary',
'''',
'transpose',
880 IF( result( j ).LT.10000.0 )
THEN
881 WRITE( nout, fmt = 9990 )nptknt, mplusn, j, result( j )
883 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
897 CALL alasvm(
'CGX', nout, nerrs, ntestt, 0 )
903 9999
FORMAT(
' CDRGSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
904 $ i6,
', JTYPE=', i6,
')' )
906 9998
FORMAT(
' CDRGSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
907 $ i6,
', Input Example #', i2,
')' )
909 9997
FORMAT(
' CDRGSX: S not in Schur form at eigenvalue ', i6,
'.',
910 $ / 9x,
'N=', i6,
', JTYPE=', i6,
')' )
912 9996
FORMAT( / 1x, a3,
' -- Complex Expert Generalized Schur form',
913 $
' problem driver' )
915 9995
FORMAT(
'Input Example' )
917 9994
FORMAT(
' Matrix types: ', /
918 $
' 1: A is a block diagonal matrix of Jordan blocks ',
919 $
'and B is the identity ', /
' matrix, ',
920 $ /
' 2: A and B are upper triangular matrices, ',
921 $ /
' 3: A and B are as type 2, but each second diagonal ',
922 $
'block in A_11 and ', /
923 $
' each third diaongal block in A_22 are 2x2 blocks,',
924 $ /
' 4: A and B are block diagonal matrices, ',
925 $ /
' 5: (A,B) has potentially close or common ',
926 $
'eigenvalues.', / )
928 9993
FORMAT( /
' Tests performed: (S is Schur, T is triangular, ',
929 $
'Q and Z are ', a,
',', / 19x,
930 $
' a is alpha, b is beta, and ', a,
' means ', a,
'.)',
931 $ /
' 1 = | A - Q S Z', a,
932 $
' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
933 $
' | / ( |B| n ulp )', /
' 3 = | I - QQ', a,
934 $
' | / ( n ulp ) 4 = | I - ZZ', a,
935 $
' | / ( n ulp )', /
' 5 = 1/ULP if A is not in ',
936 $
'Schur form S', /
' 6 = difference between (alpha,beta)',
937 $
' and diagonals of (S,T)', /
938 $
' 7 = 1/ULP if SDIM is not the correct number of ',
939 $
'selected eigenvalues', /
940 $
' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
941 $
'DIFTRU/DIFEST > 10*THRESH',
942 $ /
' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
943 $
'when reordering fails', /
944 $
' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
945 $
'PLTRU/PLEST > THRESH', /
946 $
' ( Test 10 is only for input examples )', / )
947 9992
FORMAT(
' Matrix order=', i2,
', type=', i2,
', a=', e10.3,
948 $
', order(A_11)=', i2,
', result ', i2,
' is ', 0p, f8.2 )
949 9991
FORMAT(
' Matrix order=', i2,
', type=', i2,
', a=', e10.3,
950 $
', order(A_11)=', i2,
', result ', i2,
' is ', 0p, e10.3 )
951 9990
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
952 $
' result ', i2,
' is', 0p, f8.2 )
953 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
954 $
' result ', i2,
' is', 1p, e10.3 )
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine cgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, RWORK, INFO)
CGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine cdrgsx(NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK, LWORK, RWORK, IWORK, LIWORK, BWORK, INFO)
CDRGSX
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine cggesx(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)
CGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
subroutine clakf2(M, N, A, LDA, B, D, E, Z, LDZ)
CLAKF2
subroutine cget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
CGET51
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine slabad(SMALL, LARGE)
SLABAD
subroutine clatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
CLATM5