358 SUBROUTINE ddrgsx( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI,
359 $ bi, z, q, alphar, alphai, beta, c, ldc, s,
360 $ work, lwork, iwork, liwork, bwork, info )
368 INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
370 DOUBLE PRECISION THRESH
375 DOUBLE PRECISION A( lda, * ), AI( lda, * ), ALPHAI( * ),
376 $ alphar( * ), b( lda, * ), beta( * ),
377 $ bi( lda, * ), c( ldc, * ), q( lda, * ), s( * ),
378 $ work( * ), z( lda, * )
384 DOUBLE PRECISION ZERO, ONE, TEN
385 parameter ( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1 )
390 INTEGER BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
391 $ minwrk, mm, mn2, nerrs, nptknt, ntest, ntestt,
393 DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
394 $ temp2, thrsh2, ulp, ulpinv, weight
397 DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 )
402 DOUBLE PRECISION DLAMCH, DLANGE
403 EXTERNAL dlctsx, ilaenv, dlamch, dlange
410 INTRINSIC abs, max, sqrt
414 INTEGER K, M, MPLUSN, N
417 COMMON / mn / m, n, mplusn, k, fs
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+6 )
THEN
447 IF( info.EQ.0 .AND. lwork.GE.1 )
THEN
448 minwrk = max( 10*( nsize+1 ), 5*nsize*nsize / 2 )
452 maxwrk = 9*( nsize+1 ) + nsize*
453 $ ilaenv( 1,
'DGEQRF',
' ', nsize, 1, nsize, 0 )
454 maxwrk = max( maxwrk, 9*( nsize+1 )+nsize*
455 $ ilaenv( 1,
'DORGQR',
' ', nsize, 1, nsize, -1 ) )
459 bdspac = 5*nsize*nsize / 2
460 maxwrk = max( maxwrk, 3*nsize*nsize / 2+nsize*nsize*
461 $ ilaenv( 1,
'DGEBRD',
' ', 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(
'DDRGSX', -info )
482 smlnum = dlamch(
'S' ) / ulp
483 bignum = one / smlnum
484 CALL dlabad( smlnum, bignum )
506 DO 40 m = 1, nsize - 1
507 DO 30 n = 1, nsize - m
509 weight = one / weight
517 CALL dlaset(
'Full', mplusn, mplusn, zero, zero, ai,
519 CALL dlaset(
'Full', mplusn, mplusn, zero, zero, bi,
522 CALL dlatm5( 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 dlacpy(
'Full', mplusn, mplusn, ai, lda, a, lda )
543 CALL dlacpy(
'Full', mplusn, mplusn, bi, lda, b, lda )
545 CALL dggesx(
'V',
'V',
'S', dlctsx, sense, mplusn, ai,
546 $ lda, bi, lda, mm, alphar, alphai, beta,
547 $ q, lda, z, lda, pl, difest, work, lwork,
548 $ iwork, liwork, bwork, linfo )
550 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 )
THEN
552 WRITE( nout, fmt = 9999 )
'DGGESX', linfo, mplusn,
560 CALL dlacpy(
'Full', mplusn, mplusn, ai, lda, work,
562 CALL dlacpy(
'Full', mplusn, mplusn, bi, lda,
563 $ work( mplusn*mplusn+1 ), mplusn )
564 abnrm = dlange(
'Fro', mplusn, 2*mplusn, work, mplusn,
569 CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z,
570 $ lda, work, result( 1 ) )
571 CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z,
572 $ lda, work, result( 2 ) )
573 CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q,
574 $ lda, work, result( 3 ) )
575 CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z,
576 $ lda, work, result( 4 ) )
588 IF( alphai( j ).EQ.zero )
THEN
589 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
590 $ max( smlnum, abs( alphar( j ) ),
591 $ abs( ai( j, j ) ) )+
592 $ abs( beta( j )-bi( j, j ) ) /
593 $ max( smlnum, abs( beta( j ) ),
594 $ abs( 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
608 IF( alphai( j ).GT.zero )
THEN
613 IF( i1.LE.0 .OR. i1.GE.mplusn )
THEN
615 ELSE IF( i1.LT.mplusn-1 )
THEN
616 IF( ai( i1+2, i1+1 ).NE.zero )
THEN
620 ELSE IF( i1.GT.1 )
THEN
621 IF( ai( i1, i1-1 ).NE.zero )
THEN
626 IF( .NOT.ilabad )
THEN
627 CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ),
628 $ lda, beta( j ), alphar( j ),
629 $ alphai( j ), temp2, iinfo )
630 IF( iinfo.GE.3 )
THEN
631 WRITE( nout, fmt = 9997 )iinfo, j,
639 temp1 = max( temp1, temp2 )
641 WRITE( nout, fmt = 9996 )j, mplusn, prtype
650 IF( linfo.EQ.mplusn+3 )
THEN
652 ELSE IF( mm.NE.n )
THEN
661 mn2 = mm*( mplusn-mm )*2
662 IF( ifunc.GE.2 .AND. mn2.LE.ncmax*ncmax )
THEN
667 CALL dlakf2( mm, mplusn-mm, ai, lda,
668 $ ai( mm+1, mm+1 ), bi,
669 $ bi( mm+1, mm+1 ), c, ldc )
671 CALL dgesvd(
'N',
'N', mn2, mn2, c, ldc, s, work,
672 $ 1, work( 2 ), 1, work( 3 ), lwork-2,
676 IF( difest( 2 ).EQ.zero )
THEN
677 IF( diftru.GT.abnrm*ulp )
678 $ result( 8 ) = ulpinv
679 ELSE IF( diftru.EQ.zero )
THEN
680 IF( difest( 2 ).GT.abnrm*ulp )
681 $ result( 8 ) = ulpinv
682 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
683 $ ( diftru*thrsh2.LT.difest( 2 ) ) )
THEN
684 result( 8 ) = max( diftru / difest( 2 ),
685 $ difest( 2 ) / diftru )
693 IF( linfo.EQ.( mplusn+2 ) )
THEN
694 IF( diftru.GT.abnrm*ulp )
695 $ result( 9 ) = ulpinv
696 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
697 $ result( 9 ) = ulpinv
698 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
699 $ result( 9 ) = ulpinv
703 ntestt = ntestt + ntest
708 IF( result( j ).GE.thresh )
THEN
713 IF( nerrs.EQ.0 )
THEN
714 WRITE( nout, fmt = 9995 )
'DGX'
718 WRITE( nout, fmt = 9993 )
722 WRITE( nout, fmt = 9992 )
'orthogonal',
'''',
723 $
'transpose', (
'''', i = 1, 4 )
727 IF( result( j ).LT.10000.0d0 )
THEN
728 WRITE( nout, fmt = 9991 )mplusn, prtype,
729 $ weight, m, j, result( j )
731 WRITE( nout, fmt = 9990 )mplusn, prtype,
732 $ weight, m, j, result( j )
752 READ( nin, fmt = *, end = 140 )mplusn
755 READ( nin, fmt = *, end = 140 )n
757 READ( nin, fmt = * )( ai( i, j ), j = 1, mplusn )
760 READ( nin, fmt = * )( bi( i, j ), j = 1, mplusn )
762 READ( nin, fmt = * )pltru, diftru
769 CALL dlacpy(
'Full', mplusn, mplusn, ai, lda, a, lda )
770 CALL dlacpy(
'Full', mplusn, mplusn, bi, lda, b, lda )
775 CALL dggesx(
'V',
'V',
'S', dlctsx,
'B', mplusn, ai, lda, bi, lda,
776 $ mm, alphar, alphai, beta, q, lda, z, lda, pl, difest,
777 $ work, lwork, iwork, liwork, bwork, linfo )
779 IF( linfo.NE.0 .AND. linfo.NE.mplusn+2 )
THEN
781 WRITE( nout, fmt = 9998 )
'DGGESX', linfo, mplusn, nptknt
788 CALL dlacpy(
'Full', mplusn, mplusn, ai, lda, work, mplusn )
789 CALL dlacpy(
'Full', mplusn, mplusn, bi, lda,
790 $ work( mplusn*mplusn+1 ), mplusn )
791 abnrm = dlange(
'Fro', mplusn, 2*mplusn, work, mplusn, work )
795 CALL dget51( 1, mplusn, a, lda, ai, lda, q, lda, z, lda, work,
797 CALL dget51( 1, mplusn, b, lda, bi, lda, q, lda, z, lda, work,
799 CALL dget51( 3, mplusn, b, lda, bi, lda, q, lda, q, lda, work,
801 CALL dget51( 3, mplusn, b, lda, bi, lda, z, lda, z, lda, work,
814 IF( alphai( j ).EQ.zero )
THEN
815 temp2 = ( abs( alphar( j )-ai( j, j ) ) /
816 $ max( smlnum, abs( alphar( j ) ), abs( ai( j,
817 $ j ) ) )+abs( beta( j )-bi( j, j ) ) /
818 $ max( smlnum, abs( beta( j ) ), abs( bi( j, j ) ) ) )
820 IF( j.LT.mplusn )
THEN
821 IF( ai( j+1, j ).NE.zero )
THEN
827 IF( ai( j, j-1 ).NE.zero )
THEN
833 IF( alphai( j ).GT.zero )
THEN
838 IF( i1.LE.0 .OR. i1.GE.mplusn )
THEN
840 ELSE IF( i1.LT.mplusn-1 )
THEN
841 IF( ai( i1+2, i1+1 ).NE.zero )
THEN
845 ELSE IF( i1.GT.1 )
THEN
846 IF( ai( i1, i1-1 ).NE.zero )
THEN
851 IF( .NOT.ilabad )
THEN
852 CALL dget53( ai( i1, i1 ), lda, bi( i1, i1 ), lda,
853 $ beta( j ), alphar( j ), alphai( j ), temp2,
855 IF( iinfo.GE.3 )
THEN
856 WRITE( nout, fmt = 9997 )iinfo, j, mplusn, nptknt
863 temp1 = max( temp1, temp2 )
865 WRITE( nout, fmt = 9996 )j, mplusn, nptknt
874 IF( linfo.EQ.mplusn+3 )
875 $ result( 7 ) = ulpinv
881 IF( difest( 2 ).EQ.zero )
THEN
882 IF( diftru.GT.abnrm*ulp )
883 $ result( 8 ) = ulpinv
884 ELSE IF( diftru.EQ.zero )
THEN
885 IF( difest( 2 ).GT.abnrm*ulp )
886 $ result( 8 ) = ulpinv
887 ELSE IF( ( diftru.GT.thrsh2*difest( 2 ) ) .OR.
888 $ ( diftru*thrsh2.LT.difest( 2 ) ) )
THEN
889 result( 8 ) = max( diftru / difest( 2 ), difest( 2 ) / diftru )
896 IF( linfo.EQ.( mplusn+2 ) )
THEN
897 IF( diftru.GT.abnrm*ulp )
898 $ result( 9 ) = ulpinv
899 IF( ( ifunc.GT.1 ) .AND. ( difest( 2 ).NE.zero ) )
900 $ result( 9 ) = ulpinv
901 IF( ( ifunc.EQ.1 ) .AND. ( pl( 1 ).NE.zero ) )
902 $ result( 9 ) = ulpinv
909 IF( pl( 1 ).EQ.zero )
THEN
910 IF( pltru.GT.abnrm*ulp )
911 $ result( 10 ) = ulpinv
912 ELSE IF( pltru.EQ.zero )
THEN
913 IF( pl( 1 ).GT.abnrm*ulp )
914 $ result( 10 ) = ulpinv
915 ELSE IF( ( pltru.GT.thresh*pl( 1 ) ) .OR.
916 $ ( pltru*thresh.LT.pl( 1 ) ) )
THEN
917 result( 10 ) = ulpinv
920 ntestt = ntestt + ntest
925 IF( result( j ).GE.thresh )
THEN
930 IF( nerrs.EQ.0 )
THEN
931 WRITE( nout, fmt = 9995 )
'DGX'
935 WRITE( nout, fmt = 9994 )
939 WRITE( nout, fmt = 9992 )
'orthogonal',
'''',
940 $
'transpose', (
'''', i = 1, 4 )
944 IF( result( j ).LT.10000.0d0 )
THEN
945 WRITE( nout, fmt = 9989 )nptknt, mplusn, j, result( j )
947 WRITE( nout, fmt = 9988 )nptknt, mplusn, j, result( j )
961 CALL alasvm(
'DGX', nout, nerrs, ntestt, 0 )
967 9999
FORMAT(
' DDRGSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
968 $ i6,
', JTYPE=', i6,
')' )
970 9998
FORMAT(
' DDRGSX: ', a,
' returned INFO=', i6,
'.', / 9x,
'N=',
971 $ i6,
', Input Example #', i2,
')' )
973 9997
FORMAT(
' DDRGSX: DGET53 returned INFO=', i1,
' for eigenvalue ',
974 $ i6,
'.', / 9x,
'N=', i6,
', JTYPE=', i6,
')' )
976 9996
FORMAT(
' DDRGSX: S not in Schur form at eigenvalue ', i6,
'.',
977 $ / 9x,
'N=', i6,
', JTYPE=', i6,
')' )
979 9995
FORMAT( / 1x, a3,
' -- Real Expert Generalized Schur form',
980 $
' problem driver' )
982 9994
FORMAT(
'Input Example' )
984 9993
FORMAT(
' Matrix types: ', /
985 $
' 1: A is a block diagonal matrix of Jordan blocks ',
986 $
'and B is the identity ', /
' matrix, ',
987 $ /
' 2: A and B are upper triangular matrices, ',
988 $ /
' 3: A and B are as type 2, but each second diagonal ',
989 $
'block in A_11 and ', /
990 $
' each third diaongal block in A_22 are 2x2 blocks,',
991 $ /
' 4: A and B are block diagonal matrices, ',
992 $ /
' 5: (A,B) has potentially close or common ',
993 $
'eigenvalues.', / )
995 9992
FORMAT( /
' Tests performed: (S is Schur, T is triangular, ',
996 $
'Q and Z are ', a,
',', / 19x,
997 $
' a is alpha, b is beta, and ', a,
' means ', a,
'.)',
998 $ /
' 1 = | A - Q S Z', a,
999 $
' | / ( |A| n ulp ) 2 = | B - Q T Z', a,
1000 $
' | / ( |B| n ulp )', /
' 3 = | I - QQ', a,
1001 $
' | / ( n ulp ) 4 = | I - ZZ', a,
1002 $
' | / ( n ulp )', /
' 5 = 1/ULP if A is not in ',
1003 $
'Schur form S', /
' 6 = difference between (alpha,beta)',
1004 $
' and diagonals of (S,T)', /
1005 $
' 7 = 1/ULP if SDIM is not the correct number of ',
1006 $
'selected eigenvalues', /
1007 $
' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
1008 $
'DIFTRU/DIFEST > 10*THRESH',
1009 $ /
' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
1010 $
'when reordering fails', /
1011 $
' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
1012 $
'PLTRU/PLEST > THRESH', /
1013 $
' ( Test 10 is only for input examples )', / )
1014 9991
FORMAT(
' Matrix order=', i2,
', type=', i2,
', a=', d10.3,
1015 $
', order(A_11)=', i2,
', result ', i2,
' is ', 0p, f8.2 )
1016 9990
FORMAT(
' Matrix order=', i2,
', type=', i2,
', a=', d10.3,
1017 $
', order(A_11)=', i2,
', result ', i2,
' is ', 0p, d10.3 )
1018 9989
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
1019 $
' result ', i2,
' is', 0p, f8.2 )
1020 9988
FORMAT(
' Input example #', i2,
', matrix order=', i4,
',',
1021 $
' result ', i2,
' is', 1p, d10.3 )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
subroutine dlakf2(M, N, A, LDA, B, D, E, Z, LDZ)
DLAKF2
subroutine dlatm5(PRTYPE, M, N, A, LDA, B, LDB, C, LDC, D, LDD, E, LDE, F, LDF, R, LDR, L, LDL, ALPHA, QBLCKA, QBLCKB)
DLATM5
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine xerbla(SRNAME, INFO)
XERBLA
subroutine dget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RESULT)
DGET51
subroutine dlabad(SMALL, LARGE)
DLABAD
subroutine dget53(A, LDA, B, LDB, SCALE, WR, WI, RESULT, INFO)
DGET53
subroutine dggesx(JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA, B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
DGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
subroutine dgesvd(JOBU, JOBVT, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK, LWORK, INFO)
DGESVD computes the singular value decomposition (SVD) for GE matrices
subroutine ddrgsx(NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B, AI, BI, Z, Q, ALPHAR, ALPHAI, BETA, C, LDC, S, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
DDRGSX