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 )
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 )
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 claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
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
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
logical function clctsx(ALPHA, BETA)
CLCTSX
subroutine cget51(ITYPE, N, A, LDA, B, LDB, U, LDU, V, LDV, WORK, RWORK, RESULT)
CGET51
subroutine xerbla(SRNAME, INFO)
XERBLA
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
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine clakf2(M, N, A, LDA, B, D, E, Z, LDZ)
CLAKF2
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine slabad(SMALL, LARGE)
SLABAD
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
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 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