1022 lines
35 KiB
Fortran
1022 lines
35 KiB
Fortran
*> \brief \b DDRGSX
|
|
*
|
|
* =========== DOCUMENTATION ===========
|
|
*
|
|
* Online html documentation available at
|
|
* http://www.netlib.org/lapack/explore-html/
|
|
*
|
|
* Definition:
|
|
* ===========
|
|
*
|
|
* 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 )
|
|
*
|
|
* .. Scalar Arguments ..
|
|
* INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
|
|
* $ NOUT, NSIZE
|
|
* DOUBLE PRECISION THRESH
|
|
* ..
|
|
* .. Array Arguments ..
|
|
* LOGICAL BWORK( * )
|
|
* INTEGER IWORK( * )
|
|
* DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
|
|
* $ ALPHAR( * ), B( LDA, * ), BETA( * ),
|
|
* $ BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
|
|
* $ WORK( * ), Z( LDA, * )
|
|
* ..
|
|
*
|
|
*
|
|
*> \par Purpose:
|
|
* =============
|
|
*>
|
|
*> \verbatim
|
|
*>
|
|
*> DDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
|
|
*> problem expert driver DGGESX.
|
|
*>
|
|
*> DGGESX factors A and B as Q S Z' and Q T Z', where ' means
|
|
*> transpose, T is upper triangular, S is in generalized Schur form
|
|
*> (block upper triangular, with 1x1 and 2x2 blocks on the diagonal,
|
|
*> the 2x2 blocks corresponding to complex conjugate pairs of
|
|
*> generalized eigenvalues), and Q and Z are orthogonal. It also
|
|
*> computes the generalized eigenvalues (alpha(1),beta(1)), ...,
|
|
*> (alpha(n),beta(n)). Thus, w(j) = alpha(j)/beta(j) is a root of the
|
|
*> characteristic equation
|
|
*>
|
|
*> det( A - w(j) B ) = 0
|
|
*>
|
|
*> Optionally it also reorders the eigenvalues so that a selected
|
|
*> cluster of eigenvalues appears in the leading diagonal block of the
|
|
*> Schur forms; computes a reciprocal condition number for the average
|
|
*> of the selected eigenvalues; and computes a reciprocal condition
|
|
*> number for the right and left deflating subspaces corresponding to
|
|
*> the selected eigenvalues.
|
|
*>
|
|
*> When DDRGSX is called with NSIZE > 0, five (5) types of built-in
|
|
*> matrix pairs are used to test the routine DGGESX.
|
|
*>
|
|
*> When DDRGSX is called with NSIZE = 0, it reads in test matrix data
|
|
*> to test DGGESX.
|
|
*>
|
|
*> For each matrix pair, the following tests will be performed and
|
|
*> compared with the threshold THRESH except for the tests (7) and (9):
|
|
*>
|
|
*> (1) | A - Q S Z' | / ( |A| n ulp )
|
|
*>
|
|
*> (2) | B - Q T Z' | / ( |B| n ulp )
|
|
*>
|
|
*> (3) | I - QQ' | / ( n ulp )
|
|
*>
|
|
*> (4) | I - ZZ' | / ( n ulp )
|
|
*>
|
|
*> (5) if A is in Schur form (i.e. quasi-triangular form)
|
|
*>
|
|
*> (6) maximum over j of D(j) where:
|
|
*>
|
|
*> if alpha(j) is real:
|
|
*> |alpha(j) - S(j,j)| |beta(j) - T(j,j)|
|
|
*> D(j) = ------------------------ + -----------------------
|
|
*> max(|alpha(j)|,|S(j,j)|) max(|beta(j)|,|T(j,j)|)
|
|
*>
|
|
*> if alpha(j) is complex:
|
|
*> | det( s S - w T ) |
|
|
*> D(j) = ---------------------------------------------------
|
|
*> ulp max( s norm(S), |w| norm(T) )*norm( s S - w T )
|
|
*>
|
|
*> and S and T are here the 2 x 2 diagonal blocks of S and T
|
|
*> corresponding to the j-th and j+1-th eigenvalues.
|
|
*>
|
|
*> (7) if sorting worked and SDIM is the number of eigenvalues
|
|
*> which were selected.
|
|
*>
|
|
*> (8) the estimated value DIF does not differ from the true values of
|
|
*> Difu and Difl more than a factor 10*THRESH. If the estimate DIF
|
|
*> equals zero the corresponding true values of Difu and Difl
|
|
*> should be less than EPS*norm(A, B). If the true value of Difu
|
|
*> and Difl equal zero, the estimate DIF should be less than
|
|
*> EPS*norm(A, B).
|
|
*>
|
|
*> (9) If INFO = N+3 is returned by DGGESX, the reordering "failed"
|
|
*> and we check that DIF = PL = PR = 0 and that the true value of
|
|
*> Difu and Difl is < EPS*norm(A, B). We count the events when
|
|
*> INFO=N+3.
|
|
*>
|
|
*> For read-in test matrices, the above tests are run except that the
|
|
*> exact value for DIF (and PL) is input data. Additionally, there is
|
|
*> one more test run for read-in test matrices:
|
|
*>
|
|
*> (10) the estimated value PL does not differ from the true value of
|
|
*> PLTRU more than a factor THRESH. If the estimate PL equals
|
|
*> zero the corresponding true value of PLTRU should be less than
|
|
*> EPS*norm(A, B). If the true value of PLTRU equal zero, the
|
|
*> estimate PL should be less than EPS*norm(A, B).
|
|
*>
|
|
*> Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
|
|
*> matrix pairs are generated and tested. NSIZE should be kept small.
|
|
*>
|
|
*> SVD (routine DGESVD) is used for computing the true value of DIF_u
|
|
*> and DIF_l when testing the built-in test problems.
|
|
*>
|
|
*> Built-in Test Matrices
|
|
*> ======================
|
|
*>
|
|
*> All built-in test matrices are the 2 by 2 block of triangular
|
|
*> matrices
|
|
*>
|
|
*> A = [ A11 A12 ] and B = [ B11 B12 ]
|
|
*> [ A22 ] [ B22 ]
|
|
*>
|
|
*> where for different type of A11 and A22 are given as the following.
|
|
*> A12 and B12 are chosen so that the generalized Sylvester equation
|
|
*>
|
|
*> A11*R - L*A22 = -A12
|
|
*> B11*R - L*B22 = -B12
|
|
*>
|
|
*> have prescribed solution R and L.
|
|
*>
|
|
*> Type 1: A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
|
|
*> B11 = I_m, B22 = I_k
|
|
*> where J_k(a,b) is the k-by-k Jordan block with ``a'' on
|
|
*> diagonal and ``b'' on superdiagonal.
|
|
*>
|
|
*> Type 2: A11 = (a_ij) = ( 2(.5-sin(i)) ) and
|
|
*> B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
|
|
*> A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
|
|
*> B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
|
|
*>
|
|
*> Type 3: A11, A22 and B11, B22 are chosen as for Type 2, but each
|
|
*> second diagonal block in A_11 and each third diagonal block
|
|
*> in A_22 are made as 2 by 2 blocks.
|
|
*>
|
|
*> Type 4: A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
|
|
*> for i=1,...,m, j=1,...,m and
|
|
*> A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
|
|
*> for i=m+1,...,k, j=m+1,...,k
|
|
*>
|
|
*> Type 5: (A,B) and have potentially close or common eigenvalues and
|
|
*> very large departure from block diagonality A_11 is chosen
|
|
*> as the m x m leading submatrix of A_1:
|
|
*> | 1 b |
|
|
*> | -b 1 |
|
|
*> | 1+d b |
|
|
*> | -b 1+d |
|
|
*> A_1 = | d 1 |
|
|
*> | -1 d |
|
|
*> | -d 1 |
|
|
*> | -1 -d |
|
|
*> | 1 |
|
|
*> and A_22 is chosen as the k x k leading submatrix of A_2:
|
|
*> | -1 b |
|
|
*> | -b -1 |
|
|
*> | 1-d b |
|
|
*> | -b 1-d |
|
|
*> A_2 = | d 1+b |
|
|
*> | -1-b d |
|
|
*> | -d 1+b |
|
|
*> | -1+b -d |
|
|
*> | 1-d |
|
|
*> and matrix B are chosen as identity matrices (see DLATM5).
|
|
*>
|
|
*> \endverbatim
|
|
*
|
|
* Arguments:
|
|
* ==========
|
|
*
|
|
*> \param[in] NSIZE
|
|
*> \verbatim
|
|
*> NSIZE is INTEGER
|
|
*> The maximum size of the matrices to use. NSIZE >= 0.
|
|
*> If NSIZE = 0, no built-in tests matrices are used, but
|
|
*> read-in test matrices are used to test DGGESX.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NCMAX
|
|
*> \verbatim
|
|
*> NCMAX is INTEGER
|
|
*> Maximum allowable NMAX for generating Kroneker matrix
|
|
*> in call to DLAKF2
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] THRESH
|
|
*> \verbatim
|
|
*> THRESH is DOUBLE PRECISION
|
|
*> A test will count as "failed" if the "error", computed as
|
|
*> described above, exceeds THRESH. Note that the error
|
|
*> is scaled to be O(1), so THRESH should be a reasonably
|
|
*> small multiple of 1, e.g., 10 or 100. In particular,
|
|
*> it should not depend on the precision (single vs. double)
|
|
*> or the size of the matrix. THRESH >= 0.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NIN
|
|
*> \verbatim
|
|
*> NIN is INTEGER
|
|
*> The FORTRAN unit number for reading in the data file of
|
|
*> problems to solve.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] NOUT
|
|
*> \verbatim
|
|
*> NOUT is INTEGER
|
|
*> The FORTRAN unit number for printing out error messages
|
|
*> (e.g., if a routine returns IINFO not equal to 0.)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] A
|
|
*> \verbatim
|
|
*> A is DOUBLE PRECISION array, dimension (LDA, NSIZE)
|
|
*> Used to store the matrix whose eigenvalues are to be
|
|
*> computed. On exit, A contains the last matrix actually used.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDA
|
|
*> \verbatim
|
|
*> LDA is INTEGER
|
|
*> The leading dimension of A, B, AI, BI, Z and Q,
|
|
*> LDA >= max( 1, NSIZE ). For the read-in test,
|
|
*> LDA >= max( 1, N ), N is the size of the test matrices.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] B
|
|
*> \verbatim
|
|
*> B is DOUBLE PRECISION array, dimension (LDA, NSIZE)
|
|
*> Used to store the matrix whose eigenvalues are to be
|
|
*> computed. On exit, B contains the last matrix actually used.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] AI
|
|
*> \verbatim
|
|
*> AI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
|
|
*> Copy of A, modified by DGGESX.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] BI
|
|
*> \verbatim
|
|
*> BI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
|
|
*> Copy of B, modified by DGGESX.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] Z
|
|
*> \verbatim
|
|
*> Z is DOUBLE PRECISION array, dimension (LDA, NSIZE)
|
|
*> Z holds the left Schur vectors computed by DGGESX.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] Q
|
|
*> \verbatim
|
|
*> Q is DOUBLE PRECISION array, dimension (LDA, NSIZE)
|
|
*> Q holds the right Schur vectors computed by DGGESX.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] ALPHAR
|
|
*> \verbatim
|
|
*> ALPHAR is DOUBLE PRECISION array, dimension (NSIZE)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] ALPHAI
|
|
*> \verbatim
|
|
*> ALPHAI is DOUBLE PRECISION array, dimension (NSIZE)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] BETA
|
|
*> \verbatim
|
|
*> BETA is DOUBLE PRECISION array, dimension (NSIZE)
|
|
*>
|
|
*> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] C
|
|
*> \verbatim
|
|
*> C is DOUBLE PRECISION array, dimension (LDC, LDC)
|
|
*> Store the matrix generated by subroutine DLAKF2, this is the
|
|
*> matrix formed by Kronecker products used for estimating
|
|
*> DIF.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LDC
|
|
*> \verbatim
|
|
*> LDC is INTEGER
|
|
*> The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] S
|
|
*> \verbatim
|
|
*> S is DOUBLE PRECISION array, dimension (LDC)
|
|
*> Singular values of C
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] WORK
|
|
*> \verbatim
|
|
*> WORK is DOUBLE PRECISION array, dimension (LWORK)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LWORK
|
|
*> \verbatim
|
|
*> LWORK is INTEGER
|
|
*> The dimension of the array WORK.
|
|
*> LWORK >= MAX( 5*NSIZE*NSIZE/2 - 2, 10*(NSIZE+1) )
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] IWORK
|
|
*> \verbatim
|
|
*> IWORK is INTEGER array, dimension (LIWORK)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[in] LIWORK
|
|
*> \verbatim
|
|
*> LIWORK is INTEGER
|
|
*> The dimension of the array IWORK. LIWORK >= NSIZE + 6.
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] BWORK
|
|
*> \verbatim
|
|
*> BWORK is LOGICAL array, dimension (LDA)
|
|
*> \endverbatim
|
|
*>
|
|
*> \param[out] INFO
|
|
*> \verbatim
|
|
*> INFO is INTEGER
|
|
*> = 0: successful exit
|
|
*> < 0: if INFO = -i, the i-th argument had an illegal value.
|
|
*> > 0: A routine returned an error code.
|
|
*> \endverbatim
|
|
*
|
|
* Authors:
|
|
* ========
|
|
*
|
|
*> \author Univ. of Tennessee
|
|
*> \author Univ. of California Berkeley
|
|
*> \author Univ. of Colorado Denver
|
|
*> \author NAG Ltd.
|
|
*
|
|
*> \ingroup double_eig
|
|
*
|
|
* =====================================================================
|
|
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 )
|
|
*
|
|
* -- LAPACK test routine --
|
|
* -- LAPACK is a software package provided by Univ. of Tennessee, --
|
|
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
|
|
*
|
|
* .. Scalar Arguments ..
|
|
INTEGER INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
|
|
$ NOUT, NSIZE
|
|
DOUBLE PRECISION THRESH
|
|
* ..
|
|
* .. Array Arguments ..
|
|
LOGICAL BWORK( * )
|
|
INTEGER IWORK( * )
|
|
DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
|
|
$ ALPHAR( * ), B( LDA, * ), BETA( * ),
|
|
$ BI( LDA, * ), C( LDC, * ), Q( LDA, * ), S( * ),
|
|
$ WORK( * ), Z( LDA, * )
|
|
* ..
|
|
*
|
|
* =====================================================================
|
|
*
|
|
* .. Parameters ..
|
|
DOUBLE PRECISION ZERO, ONE, TEN
|
|
PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, TEN = 1.0D+1 )
|
|
* ..
|
|
* .. Local Scalars ..
|
|
LOGICAL ILABAD
|
|
CHARACTER SENSE
|
|
INTEGER BDSPAC, I, I1, IFUNC, IINFO, J, LINFO, MAXWRK,
|
|
$ MINWRK, MM, MN2, NERRS, NPTKNT, NTEST, NTESTT,
|
|
$ PRTYPE, QBA, QBB
|
|
DOUBLE PRECISION ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
|
|
$ TEMP2, THRSH2, ULP, ULPINV, WEIGHT
|
|
* ..
|
|
* .. Local Arrays ..
|
|
DOUBLE PRECISION DIFEST( 2 ), PL( 2 ), RESULT( 10 )
|
|
* ..
|
|
* .. External Functions ..
|
|
LOGICAL DLCTSX
|
|
INTEGER ILAENV
|
|
DOUBLE PRECISION DLAMCH, DLANGE
|
|
EXTERNAL DLCTSX, ILAENV, DLAMCH, DLANGE
|
|
* ..
|
|
* .. External Subroutines ..
|
|
EXTERNAL ALASVM, DGESVD, DGET51, DGET53, DGGESX,
|
|
$ DLACPY, DLAKF2, DLASET, DLATM5, XERBLA
|
|
* ..
|
|
* .. Intrinsic Functions ..
|
|
INTRINSIC ABS, MAX, SQRT
|
|
* ..
|
|
* .. Scalars in Common ..
|
|
LOGICAL FS
|
|
INTEGER K, M, MPLUSN, N
|
|
* ..
|
|
* .. Common blocks ..
|
|
COMMON / MN / M, N, MPLUSN, K, FS
|
|
* ..
|
|
* .. Executable Statements ..
|
|
*
|
|
* Check for errors
|
|
*
|
|
IF( NSIZE.LT.0 ) THEN
|
|
INFO = -1
|
|
ELSE IF( THRESH.LT.ZERO ) THEN
|
|
INFO = -2
|
|
ELSE IF( NIN.LE.0 ) THEN
|
|
INFO = -3
|
|
ELSE IF( NOUT.LE.0 ) THEN
|
|
INFO = -4
|
|
ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN
|
|
INFO = -6
|
|
ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN
|
|
INFO = -17
|
|
ELSE IF( LIWORK.LT.NSIZE+6 ) THEN
|
|
INFO = -21
|
|
END IF
|
|
*
|
|
* Compute workspace
|
|
* (Note: Comments in the code beginning "Workspace:" describe the
|
|
* minimal amount of workspace needed at that point in the code,
|
|
* as well as the preferred amount for good performance.
|
|
* NB refers to the optimal block size for the immediately
|
|
* following subroutine, as returned by ILAENV.)
|
|
*
|
|
MINWRK = 1
|
|
IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
|
|
MINWRK = MAX( 10*( NSIZE+1 ), 5*NSIZE*NSIZE / 2 )
|
|
*
|
|
* workspace for sggesx
|
|
*
|
|
MAXWRK = 9*( NSIZE+1 ) + NSIZE*
|
|
$ ILAENV( 1, 'DGEQRF', ' ', NSIZE, 1, NSIZE, 0 )
|
|
MAXWRK = MAX( MAXWRK, 9*( NSIZE+1 )+NSIZE*
|
|
$ ILAENV( 1, 'DORGQR', ' ', NSIZE, 1, NSIZE, -1 ) )
|
|
*
|
|
* workspace for dgesvd
|
|
*
|
|
BDSPAC = 5*NSIZE*NSIZE / 2
|
|
MAXWRK = MAX( MAXWRK, 3*NSIZE*NSIZE / 2+NSIZE*NSIZE*
|
|
$ ILAENV( 1, 'DGEBRD', ' ', NSIZE*NSIZE / 2,
|
|
$ NSIZE*NSIZE / 2, -1, -1 ) )
|
|
MAXWRK = MAX( MAXWRK, BDSPAC )
|
|
*
|
|
MAXWRK = MAX( MAXWRK, MINWRK )
|
|
*
|
|
WORK( 1 ) = MAXWRK
|
|
END IF
|
|
*
|
|
IF( LWORK.LT.MINWRK )
|
|
$ INFO = -19
|
|
*
|
|
IF( INFO.NE.0 ) THEN
|
|
CALL XERBLA( 'DDRGSX', -INFO )
|
|
RETURN
|
|
END IF
|
|
*
|
|
* Important constants
|
|
*
|
|
ULP = DLAMCH( 'P' )
|
|
ULPINV = ONE / ULP
|
|
SMLNUM = DLAMCH( 'S' ) / ULP
|
|
BIGNUM = ONE / SMLNUM
|
|
THRSH2 = TEN*THRESH
|
|
NTESTT = 0
|
|
NERRS = 0
|
|
*
|
|
* Go to the tests for read-in matrix pairs
|
|
*
|
|
IFUNC = 0
|
|
IF( NSIZE.EQ.0 )
|
|
$ GO TO 70
|
|
*
|
|
* Test the built-in matrix pairs.
|
|
* Loop over different functions (IFUNC) of DGGESX, types (PRTYPE)
|
|
* of test matrices, different size (M+N)
|
|
*
|
|
PRTYPE = 0
|
|
QBA = 3
|
|
QBB = 4
|
|
WEIGHT = SQRT( ULP )
|
|
*
|
|
DO 60 IFUNC = 0, 3
|
|
DO 50 PRTYPE = 1, 5
|
|
DO 40 M = 1, NSIZE - 1
|
|
DO 30 N = 1, NSIZE - M
|
|
*
|
|
WEIGHT = ONE / WEIGHT
|
|
MPLUSN = M + N
|
|
*
|
|
* Generate test matrices
|
|
*
|
|
FS = .TRUE.
|
|
K = 0
|
|
*
|
|
CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, AI,
|
|
$ LDA )
|
|
CALL DLASET( 'Full', MPLUSN, MPLUSN, ZERO, ZERO, BI,
|
|
$ LDA )
|
|
*
|
|
CALL DLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
|
|
$ LDA, AI( 1, M+1 ), LDA, BI, LDA,
|
|
$ BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
|
|
$ Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
|
|
*
|
|
* Compute the Schur factorization and swapping the
|
|
* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
|
|
* Swapping is accomplished via the function DLCTSX
|
|
* which is supplied below.
|
|
*
|
|
IF( IFUNC.EQ.0 ) THEN
|
|
SENSE = 'N'
|
|
ELSE IF( IFUNC.EQ.1 ) THEN
|
|
SENSE = 'E'
|
|
ELSE IF( IFUNC.EQ.2 ) THEN
|
|
SENSE = 'V'
|
|
ELSE IF( IFUNC.EQ.3 ) THEN
|
|
SENSE = 'B'
|
|
END IF
|
|
*
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
|
|
*
|
|
CALL DGGESX( 'V', 'V', 'S', DLCTSX, SENSE, MPLUSN, AI,
|
|
$ LDA, BI, LDA, MM, ALPHAR, ALPHAI, BETA,
|
|
$ Q, LDA, Z, LDA, PL, DIFEST, WORK, LWORK,
|
|
$ IWORK, LIWORK, BWORK, LINFO )
|
|
*
|
|
IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
|
|
RESULT( 1 ) = ULPINV
|
|
WRITE( NOUT, FMT = 9999 )'DGGESX', LINFO, MPLUSN,
|
|
$ PRTYPE
|
|
INFO = LINFO
|
|
GO TO 30
|
|
END IF
|
|
*
|
|
* Compute the norm(A, B)
|
|
*
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
|
|
$ MPLUSN )
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
|
|
$ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
|
|
ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
|
|
$ WORK )
|
|
*
|
|
* Do tests (1) to (4)
|
|
*
|
|
CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
|
|
$ LDA, WORK, RESULT( 1 ) )
|
|
CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
|
|
$ LDA, WORK, RESULT( 2 ) )
|
|
CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
|
|
$ LDA, WORK, RESULT( 3 ) )
|
|
CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
|
|
$ LDA, WORK, RESULT( 4 ) )
|
|
NTEST = 4
|
|
*
|
|
* Do tests (5) and (6): check Schur form of A and
|
|
* compare eigenvalues with diagonals.
|
|
*
|
|
TEMP1 = ZERO
|
|
RESULT( 5 ) = ZERO
|
|
RESULT( 6 ) = ZERO
|
|
*
|
|
DO 10 J = 1, MPLUSN
|
|
ILABAD = .FALSE.
|
|
IF( ALPHAI( J ).EQ.ZERO ) THEN
|
|
TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
|
|
$ MAX( SMLNUM, ABS( ALPHAR( J ) ),
|
|
$ ABS( AI( J, J ) ) )+
|
|
$ ABS( BETA( J )-BI( J, J ) ) /
|
|
$ MAX( SMLNUM, ABS( BETA( J ) ),
|
|
$ ABS( BI( J, J ) ) ) ) / ULP
|
|
IF( J.LT.MPLUSN ) THEN
|
|
IF( AI( J+1, J ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
END IF
|
|
IF( J.GT.1 ) THEN
|
|
IF( AI( J, J-1 ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
END IF
|
|
ELSE
|
|
IF( ALPHAI( J ).GT.ZERO ) THEN
|
|
I1 = J
|
|
ELSE
|
|
I1 = J - 1
|
|
END IF
|
|
IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
|
|
ILABAD = .TRUE.
|
|
ELSE IF( I1.LT.MPLUSN-1 ) THEN
|
|
IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
ELSE IF( I1.GT.1 ) THEN
|
|
IF( AI( I1, I1-1 ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
END IF
|
|
IF( .NOT.ILABAD ) THEN
|
|
CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ),
|
|
$ LDA, BETA( J ), ALPHAR( J ),
|
|
$ ALPHAI( J ), TEMP2, IINFO )
|
|
IF( IINFO.GE.3 ) THEN
|
|
WRITE( NOUT, FMT = 9997 )IINFO, J,
|
|
$ MPLUSN, PRTYPE
|
|
INFO = ABS( IINFO )
|
|
END IF
|
|
ELSE
|
|
TEMP2 = ULPINV
|
|
END IF
|
|
END IF
|
|
TEMP1 = MAX( TEMP1, TEMP2 )
|
|
IF( ILABAD ) THEN
|
|
WRITE( NOUT, FMT = 9996 )J, MPLUSN, PRTYPE
|
|
END IF
|
|
10 CONTINUE
|
|
RESULT( 6 ) = TEMP1
|
|
NTEST = NTEST + 2
|
|
*
|
|
* Test (7) (if sorting worked)
|
|
*
|
|
RESULT( 7 ) = ZERO
|
|
IF( LINFO.EQ.MPLUSN+3 ) THEN
|
|
RESULT( 7 ) = ULPINV
|
|
ELSE IF( MM.NE.N ) THEN
|
|
RESULT( 7 ) = ULPINV
|
|
END IF
|
|
NTEST = NTEST + 1
|
|
*
|
|
* Test (8): compare the estimated value DIF and its
|
|
* value. first, compute the exact DIF.
|
|
*
|
|
RESULT( 8 ) = ZERO
|
|
MN2 = MM*( MPLUSN-MM )*2
|
|
IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
|
|
*
|
|
* Note: for either following two causes, there are
|
|
* almost same number of test cases fail the test.
|
|
*
|
|
CALL DLAKF2( MM, MPLUSN-MM, AI, LDA,
|
|
$ AI( MM+1, MM+1 ), BI,
|
|
$ BI( MM+1, MM+1 ), C, LDC )
|
|
*
|
|
CALL DGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
|
|
$ 1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
|
|
$ INFO )
|
|
DIFTRU = S( MN2 )
|
|
*
|
|
IF( DIFEST( 2 ).EQ.ZERO ) THEN
|
|
IF( DIFTRU.GT.ABNRM*ULP )
|
|
$ RESULT( 8 ) = ULPINV
|
|
ELSE IF( DIFTRU.EQ.ZERO ) THEN
|
|
IF( DIFEST( 2 ).GT.ABNRM*ULP )
|
|
$ RESULT( 8 ) = ULPINV
|
|
ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
|
|
$ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
|
|
RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
|
|
$ DIFEST( 2 ) / DIFTRU )
|
|
END IF
|
|
NTEST = NTEST + 1
|
|
END IF
|
|
*
|
|
* Test (9)
|
|
*
|
|
RESULT( 9 ) = ZERO
|
|
IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
|
|
IF( DIFTRU.GT.ABNRM*ULP )
|
|
$ RESULT( 9 ) = ULPINV
|
|
IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
|
|
$ RESULT( 9 ) = ULPINV
|
|
IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
|
|
$ RESULT( 9 ) = ULPINV
|
|
NTEST = NTEST + 1
|
|
END IF
|
|
*
|
|
NTESTT = NTESTT + NTEST
|
|
*
|
|
* Print out tests which fail.
|
|
*
|
|
DO 20 J = 1, 9
|
|
IF( RESULT( J ).GE.THRESH ) THEN
|
|
*
|
|
* If this is the first test to fail,
|
|
* print a header to the data file.
|
|
*
|
|
IF( NERRS.EQ.0 ) THEN
|
|
WRITE( NOUT, FMT = 9995 )'DGX'
|
|
*
|
|
* Matrix types
|
|
*
|
|
WRITE( NOUT, FMT = 9993 )
|
|
*
|
|
* Tests performed
|
|
*
|
|
WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
|
|
$ 'transpose', ( '''', I = 1, 4 )
|
|
*
|
|
END IF
|
|
NERRS = NERRS + 1
|
|
IF( RESULT( J ).LT.10000.0D0 ) THEN
|
|
WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
|
|
$ WEIGHT, M, J, RESULT( J )
|
|
ELSE
|
|
WRITE( NOUT, FMT = 9990 )MPLUSN, PRTYPE,
|
|
$ WEIGHT, M, J, RESULT( J )
|
|
END IF
|
|
END IF
|
|
20 CONTINUE
|
|
*
|
|
30 CONTINUE
|
|
40 CONTINUE
|
|
50 CONTINUE
|
|
60 CONTINUE
|
|
*
|
|
GO TO 150
|
|
*
|
|
70 CONTINUE
|
|
*
|
|
* Read in data from file to check accuracy of condition estimation
|
|
* Read input data until N=0
|
|
*
|
|
NPTKNT = 0
|
|
*
|
|
80 CONTINUE
|
|
READ( NIN, FMT = *, END = 140 )MPLUSN
|
|
IF( MPLUSN.EQ.0 )
|
|
$ GO TO 140
|
|
READ( NIN, FMT = *, END = 140 )N
|
|
DO 90 I = 1, MPLUSN
|
|
READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
|
|
90 CONTINUE
|
|
DO 100 I = 1, MPLUSN
|
|
READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
|
|
100 CONTINUE
|
|
READ( NIN, FMT = * )PLTRU, DIFTRU
|
|
*
|
|
NPTKNT = NPTKNT + 1
|
|
FS = .TRUE.
|
|
K = 0
|
|
M = MPLUSN - N
|
|
*
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
|
|
*
|
|
* Compute the Schur factorization while swapping the
|
|
* m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
|
|
*
|
|
CALL DGGESX( 'V', 'V', 'S', DLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
|
|
$ MM, ALPHAR, ALPHAI, BETA, Q, LDA, Z, LDA, PL, DIFEST,
|
|
$ WORK, LWORK, IWORK, LIWORK, BWORK, LINFO )
|
|
*
|
|
IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
|
|
RESULT( 1 ) = ULPINV
|
|
WRITE( NOUT, FMT = 9998 )'DGGESX', LINFO, MPLUSN, NPTKNT
|
|
GO TO 130
|
|
END IF
|
|
*
|
|
* Compute the norm(A, B)
|
|
* (should this be norm of (A,B) or (AI,BI)?)
|
|
*
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
|
|
CALL DLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
|
|
$ WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
|
|
ABNRM = DLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, WORK )
|
|
*
|
|
* Do tests (1) to (4)
|
|
*
|
|
CALL DGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
|
|
$ RESULT( 1 ) )
|
|
CALL DGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
|
|
$ RESULT( 2 ) )
|
|
CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
|
|
$ RESULT( 3 ) )
|
|
CALL DGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
|
|
$ RESULT( 4 ) )
|
|
*
|
|
* Do tests (5) and (6): check Schur form of A and compare
|
|
* eigenvalues with diagonals.
|
|
*
|
|
NTEST = 6
|
|
TEMP1 = ZERO
|
|
RESULT( 5 ) = ZERO
|
|
RESULT( 6 ) = ZERO
|
|
*
|
|
DO 110 J = 1, MPLUSN
|
|
ILABAD = .FALSE.
|
|
IF( ALPHAI( J ).EQ.ZERO ) THEN
|
|
TEMP2 = ( ABS( ALPHAR( J )-AI( J, J ) ) /
|
|
$ MAX( SMLNUM, ABS( ALPHAR( J ) ), ABS( AI( J,
|
|
$ J ) ) )+ABS( BETA( J )-BI( J, J ) ) /
|
|
$ MAX( SMLNUM, ABS( BETA( J ) ), ABS( BI( J, J ) ) ) )
|
|
$ / ULP
|
|
IF( J.LT.MPLUSN ) THEN
|
|
IF( AI( J+1, J ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
END IF
|
|
IF( J.GT.1 ) THEN
|
|
IF( AI( J, J-1 ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
END IF
|
|
ELSE
|
|
IF( ALPHAI( J ).GT.ZERO ) THEN
|
|
I1 = J
|
|
ELSE
|
|
I1 = J - 1
|
|
END IF
|
|
IF( I1.LE.0 .OR. I1.GE.MPLUSN ) THEN
|
|
ILABAD = .TRUE.
|
|
ELSE IF( I1.LT.MPLUSN-1 ) THEN
|
|
IF( AI( I1+2, I1+1 ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
ELSE IF( I1.GT.1 ) THEN
|
|
IF( AI( I1, I1-1 ).NE.ZERO ) THEN
|
|
ILABAD = .TRUE.
|
|
RESULT( 5 ) = ULPINV
|
|
END IF
|
|
END IF
|
|
IF( .NOT.ILABAD ) THEN
|
|
CALL DGET53( AI( I1, I1 ), LDA, BI( I1, I1 ), LDA,
|
|
$ BETA( J ), ALPHAR( J ), ALPHAI( J ), TEMP2,
|
|
$ IINFO )
|
|
IF( IINFO.GE.3 ) THEN
|
|
WRITE( NOUT, FMT = 9997 )IINFO, J, MPLUSN, NPTKNT
|
|
INFO = ABS( IINFO )
|
|
END IF
|
|
ELSE
|
|
TEMP2 = ULPINV
|
|
END IF
|
|
END IF
|
|
TEMP1 = MAX( TEMP1, TEMP2 )
|
|
IF( ILABAD ) THEN
|
|
WRITE( NOUT, FMT = 9996 )J, MPLUSN, NPTKNT
|
|
END IF
|
|
110 CONTINUE
|
|
RESULT( 6 ) = TEMP1
|
|
*
|
|
* Test (7) (if sorting worked) <--------- need to be checked.
|
|
*
|
|
NTEST = 7
|
|
RESULT( 7 ) = ZERO
|
|
IF( LINFO.EQ.MPLUSN+3 )
|
|
$ RESULT( 7 ) = ULPINV
|
|
*
|
|
* Test (8): compare the estimated value of DIF and its true value.
|
|
*
|
|
NTEST = 8
|
|
RESULT( 8 ) = ZERO
|
|
IF( DIFEST( 2 ).EQ.ZERO ) THEN
|
|
IF( DIFTRU.GT.ABNRM*ULP )
|
|
$ RESULT( 8 ) = ULPINV
|
|
ELSE IF( DIFTRU.EQ.ZERO ) THEN
|
|
IF( DIFEST( 2 ).GT.ABNRM*ULP )
|
|
$ RESULT( 8 ) = ULPINV
|
|
ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
|
|
$ ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
|
|
RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
|
|
END IF
|
|
*
|
|
* Test (9)
|
|
*
|
|
NTEST = 9
|
|
RESULT( 9 ) = ZERO
|
|
IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
|
|
IF( DIFTRU.GT.ABNRM*ULP )
|
|
$ RESULT( 9 ) = ULPINV
|
|
IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
|
|
$ RESULT( 9 ) = ULPINV
|
|
IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
|
|
$ RESULT( 9 ) = ULPINV
|
|
END IF
|
|
*
|
|
* Test (10): compare the estimated value of PL and it true value.
|
|
*
|
|
NTEST = 10
|
|
RESULT( 10 ) = ZERO
|
|
IF( PL( 1 ).EQ.ZERO ) THEN
|
|
IF( PLTRU.GT.ABNRM*ULP )
|
|
$ RESULT( 10 ) = ULPINV
|
|
ELSE IF( PLTRU.EQ.ZERO ) THEN
|
|
IF( PL( 1 ).GT.ABNRM*ULP )
|
|
$ RESULT( 10 ) = ULPINV
|
|
ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
|
|
$ ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
|
|
RESULT( 10 ) = ULPINV
|
|
END IF
|
|
*
|
|
NTESTT = NTESTT + NTEST
|
|
*
|
|
* Print out tests which fail.
|
|
*
|
|
DO 120 J = 1, NTEST
|
|
IF( RESULT( J ).GE.THRESH ) THEN
|
|
*
|
|
* If this is the first test to fail,
|
|
* print a header to the data file.
|
|
*
|
|
IF( NERRS.EQ.0 ) THEN
|
|
WRITE( NOUT, FMT = 9995 )'DGX'
|
|
*
|
|
* Matrix types
|
|
*
|
|
WRITE( NOUT, FMT = 9994 )
|
|
*
|
|
* Tests performed
|
|
*
|
|
WRITE( NOUT, FMT = 9992 )'orthogonal', '''',
|
|
$ 'transpose', ( '''', I = 1, 4 )
|
|
*
|
|
END IF
|
|
NERRS = NERRS + 1
|
|
IF( RESULT( J ).LT.10000.0D0 ) THEN
|
|
WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
|
|
ELSE
|
|
WRITE( NOUT, FMT = 9988 )NPTKNT, MPLUSN, J, RESULT( J )
|
|
END IF
|
|
END IF
|
|
*
|
|
120 CONTINUE
|
|
*
|
|
130 CONTINUE
|
|
GO TO 80
|
|
140 CONTINUE
|
|
*
|
|
150 CONTINUE
|
|
*
|
|
* Summary
|
|
*
|
|
CALL ALASVM( 'DGX', NOUT, NERRS, NTESTT, 0 )
|
|
*
|
|
WORK( 1 ) = MAXWRK
|
|
*
|
|
RETURN
|
|
*
|
|
9999 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
|
|
$ I6, ', JTYPE=', I6, ')' )
|
|
*
|
|
9998 FORMAT( ' DDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
|
|
$ I6, ', Input Example #', I2, ')' )
|
|
*
|
|
9997 FORMAT( ' DDRGSX: DGET53 returned INFO=', I1, ' for eigenvalue ',
|
|
$ I6, '.', / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
|
|
*
|
|
9996 FORMAT( ' DDRGSX: S not in Schur form at eigenvalue ', I6, '.',
|
|
$ / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
|
|
*
|
|
9995 FORMAT( / 1X, A3, ' -- Real Expert Generalized Schur form',
|
|
$ ' problem driver' )
|
|
*
|
|
9994 FORMAT( 'Input Example' )
|
|
*
|
|
9993 FORMAT( ' Matrix types: ', /
|
|
$ ' 1: A is a block diagonal matrix of Jordan blocks ',
|
|
$ 'and B is the identity ', / ' matrix, ',
|
|
$ / ' 2: A and B are upper triangular matrices, ',
|
|
$ / ' 3: A and B are as type 2, but each second diagonal ',
|
|
$ 'block in A_11 and ', /
|
|
$ ' each third diagonal block in A_22 are 2x2 blocks,',
|
|
$ / ' 4: A and B are block diagonal matrices, ',
|
|
$ / ' 5: (A,B) has potentially close or common ',
|
|
$ 'eigenvalues.', / )
|
|
*
|
|
9992 FORMAT( / ' Tests performed: (S is Schur, T is triangular, ',
|
|
$ 'Q and Z are ', A, ',', / 19X,
|
|
$ ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
|
|
$ / ' 1 = | A - Q S Z', A,
|
|
$ ' | / ( |A| n ulp ) 2 = | B - Q T Z', A,
|
|
$ ' | / ( |B| n ulp )', / ' 3 = | I - QQ', A,
|
|
$ ' | / ( n ulp ) 4 = | I - ZZ', A,
|
|
$ ' | / ( n ulp )', / ' 5 = 1/ULP if A is not in ',
|
|
$ 'Schur form S', / ' 6 = difference between (alpha,beta)',
|
|
$ ' and diagonals of (S,T)', /
|
|
$ ' 7 = 1/ULP if SDIM is not the correct number of ',
|
|
$ 'selected eigenvalues', /
|
|
$ ' 8 = 1/ULP if DIFEST/DIFTRU > 10*THRESH or ',
|
|
$ 'DIFTRU/DIFEST > 10*THRESH',
|
|
$ / ' 9 = 1/ULP if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
|
|
$ 'when reordering fails', /
|
|
$ ' 10 = 1/ULP if PLEST/PLTRU > THRESH or ',
|
|
$ 'PLTRU/PLEST > THRESH', /
|
|
$ ' ( Test 10 is only for input examples )', / )
|
|
9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3,
|
|
$ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
|
|
9990 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', D10.3,
|
|
$ ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, D10.3 )
|
|
9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
|
|
$ ' result ', I2, ' is', 0P, F8.2 )
|
|
9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
|
|
$ ' result ', I2, ' is', 1P, D10.3 )
|
|
*
|
|
* End of DDRGSX
|
|
*
|
|
END
|