LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cdrgsx.f
Go to the documentation of this file.
00001 *> \brief \b CDRGSX
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *  Definition:
00009 *  ===========
00010 *
00011 *       SUBROUTINE CDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
00012 *                          AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK,
00013 *                          LWORK, RWORK, IWORK, LIWORK, BWORK, INFO )
00014 * 
00015 *       .. Scalar Arguments ..
00016 *       INTEGER            INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
00017 *      $                   NOUT, NSIZE
00018 *       REAL               THRESH
00019 *       ..
00020 *       .. Array Arguments ..
00021 *       LOGICAL            BWORK( * )
00022 *       INTEGER            IWORK( * )
00023 *       REAL               RWORK( * ), S( * )
00024 *       COMPLEX            A( LDA, * ), AI( LDA, * ), ALPHA( * ),
00025 *      $                   B( LDA, * ), BETA( * ), BI( LDA, * ),
00026 *      $                   C( LDC, * ), Q( LDA, * ), WORK( * ),
00027 *      $                   Z( LDA, * )
00028 *       ..
00029 *  
00030 *
00031 *> \par Purpose:
00032 *  =============
00033 *>
00034 *> \verbatim
00035 *>
00036 *> CDRGSX checks the nonsymmetric generalized eigenvalue (Schur form)
00037 *> problem expert driver CGGESX.
00038 *>
00039 *> CGGES factors A and B as Q*S*Z'  and Q*T*Z' , where ' means conjugate
00040 *> transpose, S and T are  upper triangular (i.e., in generalized Schur
00041 *> form), and Q and Z are unitary. It also computes the generalized
00042 *> eigenvalues (alpha(j),beta(j)), j=1,...,n.  Thus,
00043 *> w(j) = alpha(j)/beta(j) is a root of the characteristic equation
00044 *>
00045 *>                 det( A - w(j) B ) = 0
00046 *>
00047 *> Optionally it also reorders the eigenvalues so that a selected
00048 *> cluster of eigenvalues appears in the leading diagonal block of the
00049 *> Schur forms; computes a reciprocal condition number for the average
00050 *> of the selected eigenvalues; and computes a reciprocal condition
00051 *> number for the right and left deflating subspaces corresponding to
00052 *> the selected eigenvalues.
00053 *>
00054 *> When CDRGSX is called with NSIZE > 0, five (5) types of built-in
00055 *> matrix pairs are used to test the routine CGGESX.
00056 *>
00057 *> When CDRGSX is called with NSIZE = 0, it reads in test matrix data
00058 *> to test CGGESX.
00059 *> (need more details on what kind of read-in data are needed).
00060 *>
00061 *> For each matrix pair, the following tests will be performed and
00062 *> compared with the threshhold THRESH except for the tests (7) and (9):
00063 *>
00064 *> (1)   | A - Q S Z' | / ( |A| n ulp )
00065 *>
00066 *> (2)   | B - Q T Z' | / ( |B| n ulp )
00067 *>
00068 *> (3)   | I - QQ' | / ( n ulp )
00069 *>
00070 *> (4)   | I - ZZ' | / ( n ulp )
00071 *>
00072 *> (5)   if A is in Schur form (i.e. triangular form)
00073 *>
00074 *> (6)   maximum over j of D(j)  where:
00075 *>
00076 *>                     |alpha(j) - S(j,j)|        |beta(j) - T(j,j)|
00077 *>           D(j) = ------------------------ + -----------------------
00078 *>                  max(|alpha(j)|,|S(j,j)|)   max(|beta(j)|,|T(j,j)|)
00079 *>
00080 *> (7)   if sorting worked and SDIM is the number of eigenvalues
00081 *>       which were selected.
00082 *>
00083 *> (8)   the estimated value DIF does not differ from the true values of
00084 *>       Difu and Difl more than a factor 10*THRESH. If the estimate DIF
00085 *>       equals zero the corresponding true values of Difu and Difl
00086 *>       should be less than EPS*norm(A, B). If the true value of Difu
00087 *>       and Difl equal zero, the estimate DIF should be less than
00088 *>       EPS*norm(A, B).
00089 *>
00090 *> (9)   If INFO = N+3 is returned by CGGESX, the reordering "failed"
00091 *>       and we check that DIF = PL = PR = 0 and that the true value of
00092 *>       Difu and Difl is < EPS*norm(A, B). We count the events when
00093 *>       INFO=N+3.
00094 *>
00095 *> For read-in test matrices, the same tests are run except that the
00096 *> exact value for DIF (and PL) is input data.  Additionally, there is
00097 *> one more test run for read-in test matrices:
00098 *>
00099 *> (10)  the estimated value PL does not differ from the true value of
00100 *>       PLTRU more than a factor THRESH. If the estimate PL equals
00101 *>       zero the corresponding true value of PLTRU should be less than
00102 *>       EPS*norm(A, B). If the true value of PLTRU equal zero, the
00103 *>       estimate PL should be less than EPS*norm(A, B).
00104 *>
00105 *> Note that for the built-in tests, a total of 10*NSIZE*(NSIZE-1)
00106 *> matrix pairs are generated and tested. NSIZE should be kept small.
00107 *>
00108 *> SVD (routine CGESVD) is used for computing the true value of DIF_u
00109 *> and DIF_l when testing the built-in test problems.
00110 *>
00111 *> Built-in Test Matrices
00112 *> ======================
00113 *>
00114 *> All built-in test matrices are the 2 by 2 block of triangular
00115 *> matrices
00116 *>
00117 *>          A = [ A11 A12 ]    and      B = [ B11 B12 ]
00118 *>              [     A22 ]                 [     B22 ]
00119 *>
00120 *> where for different type of A11 and A22 are given as the following.
00121 *> A12 and B12 are chosen so that the generalized Sylvester equation
00122 *>
00123 *>          A11*R - L*A22 = -A12
00124 *>          B11*R - L*B22 = -B12
00125 *>
00126 *> have prescribed solution R and L.
00127 *>
00128 *> Type 1:  A11 = J_m(1,-1) and A_22 = J_k(1-a,1).
00129 *>          B11 = I_m, B22 = I_k
00130 *>          where J_k(a,b) is the k-by-k Jordan block with ``a'' on
00131 *>          diagonal and ``b'' on superdiagonal.
00132 *>
00133 *> Type 2:  A11 = (a_ij) = ( 2(.5-sin(i)) ) and
00134 *>          B11 = (b_ij) = ( 2(.5-sin(ij)) ) for i=1,...,m, j=i,...,m
00135 *>          A22 = (a_ij) = ( 2(.5-sin(i+j)) ) and
00136 *>          B22 = (b_ij) = ( 2(.5-sin(ij)) ) for i=m+1,...,k, j=i,...,k
00137 *>
00138 *> Type 3:  A11, A22 and B11, B22 are chosen as for Type 2, but each
00139 *>          second diagonal block in A_11 and each third diagonal block
00140 *>          in A_22 are made as 2 by 2 blocks.
00141 *>
00142 *> Type 4:  A11 = ( 20(.5 - sin(ij)) ) and B22 = ( 2(.5 - sin(i+j)) )
00143 *>             for i=1,...,m,  j=1,...,m and
00144 *>          A22 = ( 20(.5 - sin(i+j)) ) and B22 = ( 2(.5 - sin(ij)) )
00145 *>             for i=m+1,...,k,  j=m+1,...,k
00146 *>
00147 *> Type 5:  (A,B) and have potentially close or common eigenvalues and
00148 *>          very large departure from block diagonality A_11 is chosen
00149 *>          as the m x m leading submatrix of A_1:
00150 *>                  |  1  b                            |
00151 *>                  | -b  1                            |
00152 *>                  |        1+d  b                    |
00153 *>                  |         -b 1+d                   |
00154 *>           A_1 =  |                  d  1            |
00155 *>                  |                 -1  d            |
00156 *>                  |                        -d  1     |
00157 *>                  |                        -1 -d     |
00158 *>                  |                               1  |
00159 *>          and A_22 is chosen as the k x k leading submatrix of A_2:
00160 *>                  | -1  b                            |
00161 *>                  | -b -1                            |
00162 *>                  |       1-d  b                     |
00163 *>                  |       -b  1-d                    |
00164 *>           A_2 =  |                 d 1+b            |
00165 *>                  |               -1-b d             |
00166 *>                  |                       -d  1+b    |
00167 *>                  |                      -1+b  -d    |
00168 *>                  |                              1-d |
00169 *>          and matrix B are chosen as identity matrices (see SLATM5).
00170 *>
00171 *> \endverbatim
00172 *
00173 *  Arguments:
00174 *  ==========
00175 *
00176 *> \param[in] NSIZE
00177 *> \verbatim
00178 *>          NSIZE is INTEGER
00179 *>          The maximum size of the matrices to use. NSIZE >= 0.
00180 *>          If NSIZE = 0, no built-in tests matrices are used, but
00181 *>          read-in test matrices are used to test SGGESX.
00182 *> \endverbatim
00183 *>
00184 *> \param[in] NCMAX
00185 *> \verbatim
00186 *>          NCMAX is INTEGER
00187 *>          Maximum allowable NMAX for generating Kroneker matrix
00188 *>          in call to CLAKF2
00189 *> \endverbatim
00190 *>
00191 *> \param[in] THRESH
00192 *> \verbatim
00193 *>          THRESH is REAL
00194 *>          A test will count as "failed" if the "error", computed as
00195 *>          described above, exceeds THRESH.  Note that the error
00196 *>          is scaled to be O(1), so THRESH should be a reasonably
00197 *>          small multiple of 1, e.g., 10 or 100.  In particular,
00198 *>          it should not depend on the precision (single vs. double)
00199 *>          or the size of the matrix.  THRESH >= 0.
00200 *> \endverbatim
00201 *>
00202 *> \param[in] NIN
00203 *> \verbatim
00204 *>          NIN is INTEGER
00205 *>          The FORTRAN unit number for reading in the data file of
00206 *>          problems to solve.
00207 *> \endverbatim
00208 *>
00209 *> \param[in] NOUT
00210 *> \verbatim
00211 *>          NOUT is INTEGER
00212 *>          The FORTRAN unit number for printing out error messages
00213 *>          (e.g., if a routine returns INFO not equal to 0.)
00214 *> \endverbatim
00215 *>
00216 *> \param[out] A
00217 *> \verbatim
00218 *>          A is COMPLEX array, dimension (LDA, NSIZE)
00219 *>          Used to store the matrix whose eigenvalues are to be
00220 *>          computed.  On exit, A contains the last matrix actually used.
00221 *> \endverbatim
00222 *>
00223 *> \param[in] LDA
00224 *> \verbatim
00225 *>          LDA is INTEGER
00226 *>          The leading dimension of A, B, AI, BI, Z and Q,
00227 *>          LDA >= max( 1, NSIZE ). For the read-in test,
00228 *>          LDA >= max( 1, N ), N is the size of the test matrices.
00229 *> \endverbatim
00230 *>
00231 *> \param[out] B
00232 *> \verbatim
00233 *>          B is COMPLEX array, dimension (LDA, NSIZE)
00234 *>          Used to store the matrix whose eigenvalues are to be
00235 *>          computed.  On exit, B contains the last matrix actually used.
00236 *> \endverbatim
00237 *>
00238 *> \param[out] AI
00239 *> \verbatim
00240 *>          AI is COMPLEX array, dimension (LDA, NSIZE)
00241 *>          Copy of A, modified by CGGESX.
00242 *> \endverbatim
00243 *>
00244 *> \param[out] BI
00245 *> \verbatim
00246 *>          BI is COMPLEX array, dimension (LDA, NSIZE)
00247 *>          Copy of B, modified by CGGESX.
00248 *> \endverbatim
00249 *>
00250 *> \param[out] Z
00251 *> \verbatim
00252 *>          Z is COMPLEX array, dimension (LDA, NSIZE)
00253 *>          Z holds the left Schur vectors computed by CGGESX.
00254 *> \endverbatim
00255 *>
00256 *> \param[out] Q
00257 *> \verbatim
00258 *>          Q is COMPLEX array, dimension (LDA, NSIZE)
00259 *>          Q holds the right Schur vectors computed by CGGESX.
00260 *> \endverbatim
00261 *>
00262 *> \param[out] ALPHA
00263 *> \verbatim
00264 *>          ALPHA is COMPLEX array, dimension (NSIZE)
00265 *> \endverbatim
00266 *>
00267 *> \param[out] BETA
00268 *> \verbatim
00269 *>          BETA is COMPLEX array, dimension (NSIZE)
00270 *>
00271 *>          On exit, ALPHA/BETA are the eigenvalues.
00272 *> \endverbatim
00273 *>
00274 *> \param[out] C
00275 *> \verbatim
00276 *>          C is COMPLEX array, dimension (LDC, LDC)
00277 *>          Store the matrix generated by subroutine CLAKF2, this is the
00278 *>          matrix formed by Kronecker products used for estimating
00279 *>          DIF.
00280 *> \endverbatim
00281 *>
00282 *> \param[in] LDC
00283 *> \verbatim
00284 *>          LDC is INTEGER
00285 *>          The leading dimension of C. LDC >= max(1, LDA*LDA/2 ).
00286 *> \endverbatim
00287 *>
00288 *> \param[out] S
00289 *> \verbatim
00290 *>          S is REAL array, dimension (LDC)
00291 *>          Singular values of C
00292 *> \endverbatim
00293 *>
00294 *> \param[out] WORK
00295 *> \verbatim
00296 *>          WORK is COMPLEX array, dimension (LWORK)
00297 *> \endverbatim
00298 *>
00299 *> \param[in] LWORK
00300 *> \verbatim
00301 *>          LWORK is INTEGER
00302 *>          The dimension of the array WORK.  LWORK >= 3*NSIZE*NSIZE/2
00303 *> \endverbatim
00304 *>
00305 *> \param[out] RWORK
00306 *> \verbatim
00307 *>          RWORK is REAL array,
00308 *>                                 dimension (5*NSIZE*NSIZE/2 - 4)
00309 *> \endverbatim
00310 *>
00311 *> \param[out] IWORK
00312 *> \verbatim
00313 *>          IWORK is INTEGER array, dimension (LIWORK)
00314 *> \endverbatim
00315 *>
00316 *> \param[in] LIWORK
00317 *> \verbatim
00318 *>          LIWORK is INTEGER
00319 *>          The dimension of the array IWORK. LIWORK >= NSIZE + 2.
00320 *> \endverbatim
00321 *>
00322 *> \param[out] BWORK
00323 *> \verbatim
00324 *>          BWORK is LOGICAL array, dimension (NSIZE)
00325 *> \endverbatim
00326 *>
00327 *> \param[out] INFO
00328 *> \verbatim
00329 *>          INFO is INTEGER
00330 *>          = 0:  successful exit
00331 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00332 *>          > 0:  A routine returned an error code.
00333 *> \endverbatim
00334 *
00335 *  Authors:
00336 *  ========
00337 *
00338 *> \author Univ. of Tennessee 
00339 *> \author Univ. of California Berkeley 
00340 *> \author Univ. of Colorado Denver 
00341 *> \author NAG Ltd. 
00342 *
00343 *> \date November 2011
00344 *
00345 *> \ingroup complex_eig
00346 *
00347 *  =====================================================================
00348       SUBROUTINE CDRGSX( NSIZE, NCMAX, THRESH, NIN, NOUT, A, LDA, B,
00349      $                   AI, BI, Z, Q, ALPHA, BETA, C, LDC, S, WORK,
00350      $                   LWORK, RWORK, IWORK, LIWORK, BWORK, INFO )
00351 *
00352 *  -- LAPACK test routine (version 3.4.0) --
00353 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00354 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00355 *     November 2011
00356 *
00357 *     .. Scalar Arguments ..
00358       INTEGER            INFO, LDA, LDC, LIWORK, LWORK, NCMAX, NIN,
00359      $                   NOUT, NSIZE
00360       REAL               THRESH
00361 *     ..
00362 *     .. Array Arguments ..
00363       LOGICAL            BWORK( * )
00364       INTEGER            IWORK( * )
00365       REAL               RWORK( * ), S( * )
00366       COMPLEX            A( LDA, * ), AI( LDA, * ), ALPHA( * ),
00367      $                   B( LDA, * ), BETA( * ), BI( LDA, * ),
00368      $                   C( LDC, * ), Q( LDA, * ), WORK( * ),
00369      $                   Z( LDA, * )
00370 *     ..
00371 *
00372 *  =====================================================================
00373 *
00374 *     .. Parameters ..
00375       REAL               ZERO, ONE, TEN
00376       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1 )
00377       COMPLEX            CZERO
00378       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ) )
00379 *     ..
00380 *     .. Local Scalars ..
00381       LOGICAL            ILABAD
00382       CHARACTER          SENSE
00383       INTEGER            BDSPAC, I, IFUNC, J, LINFO, MAXWRK, MINWRK, MM,
00384      $                   MN2, NERRS, NPTKNT, NTEST, NTESTT, PRTYPE, QBA,
00385      $                   QBB
00386       REAL               ABNRM, BIGNUM, DIFTRU, PLTRU, SMLNUM, TEMP1,
00387      $                   TEMP2, THRSH2, ULP, ULPINV, WEIGHT
00388       COMPLEX            X
00389 *     ..
00390 *     .. Local Arrays ..
00391       REAL               DIFEST( 2 ), PL( 2 ), RESULT( 10 )
00392 *     ..
00393 *     .. External Functions ..
00394       LOGICAL            CLCTSX
00395       INTEGER            ILAENV
00396       REAL               CLANGE, SLAMCH
00397       EXTERNAL           CLCTSX, ILAENV, CLANGE, SLAMCH
00398 *     ..
00399 *     .. External Subroutines ..
00400       EXTERNAL           ALASVM, CGESVD, CGET51, CGGESX, CLACPY, CLAKF2,
00401      $                   CLASET, CLATM5, SLABAD, XERBLA
00402 *     ..
00403 *     .. Scalars in Common ..
00404       LOGICAL            FS
00405       INTEGER            K, M, MPLUSN, N
00406 *     ..
00407 *     .. Common blocks ..
00408       COMMON             / MN / M, N, MPLUSN, K, FS
00409 *     ..
00410 *     .. Intrinsic Functions ..
00411       INTRINSIC          ABS, AIMAG, MAX, REAL, SQRT
00412 *     ..
00413 *     .. Statement Functions ..
00414       REAL               ABS1
00415 *     ..
00416 *     .. Statement Function definitions ..
00417       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
00418 *     ..
00419 *     .. Executable Statements ..
00420 *
00421 *     Check for errors
00422 *
00423       IF( NSIZE.LT.0 ) THEN
00424          INFO = -1
00425       ELSE IF( THRESH.LT.ZERO ) THEN
00426          INFO = -2
00427       ELSE IF( NIN.LE.0 ) THEN
00428          INFO = -3
00429       ELSE IF( NOUT.LE.0 ) THEN
00430          INFO = -4
00431       ELSE IF( LDA.LT.1 .OR. LDA.LT.NSIZE ) THEN
00432          INFO = -6
00433       ELSE IF( LDC.LT.1 .OR. LDC.LT.NSIZE*NSIZE / 2 ) THEN
00434          INFO = -15
00435       ELSE IF( LIWORK.LT.NSIZE+2 ) THEN
00436          INFO = -21
00437       END IF
00438 *
00439 *     Compute workspace
00440 *      (Note: Comments in the code beginning "Workspace:" describe the
00441 *       minimal amount of workspace needed at that point in the code,
00442 *       as well as the preferred amount for good performance.
00443 *       NB refers to the optimal block size for the immediately
00444 *       following subroutine, as returned by ILAENV.)
00445 *
00446       MINWRK = 1
00447       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
00448          MINWRK = 3*NSIZE*NSIZE / 2
00449 *
00450 *        workspace for cggesx
00451 *
00452          MAXWRK = NSIZE*( 1+ILAENV( 1, 'CGEQRF', ' ', NSIZE, 1, NSIZE,
00453      $            0 ) )
00454          MAXWRK = MAX( MAXWRK, NSIZE*( 1+ILAENV( 1, 'CUNGQR', ' ',
00455      $            NSIZE, 1, NSIZE, -1 ) ) )
00456 *
00457 *        workspace for cgesvd
00458 *
00459          BDSPAC = 3*NSIZE*NSIZE / 2
00460          MAXWRK = MAX( MAXWRK, NSIZE*NSIZE*
00461      $            ( 1+ILAENV( 1, 'CGEBRD', ' ', NSIZE*NSIZE / 2,
00462      $            NSIZE*NSIZE / 2, -1, -1 ) ) )
00463          MAXWRK = MAX( MAXWRK, BDSPAC )
00464 *
00465          MAXWRK = MAX( MAXWRK, MINWRK )
00466 *
00467          WORK( 1 ) = MAXWRK
00468       END IF
00469 *
00470       IF( LWORK.LT.MINWRK )
00471      $   INFO = -18
00472 *
00473       IF( INFO.NE.0 ) THEN
00474          CALL XERBLA( 'CDRGSX', -INFO )
00475          RETURN
00476       END IF
00477 *
00478 *     Important constants
00479 *
00480       ULP = SLAMCH( 'P' )
00481       ULPINV = ONE / ULP
00482       SMLNUM = SLAMCH( 'S' ) / ULP
00483       BIGNUM = ONE / SMLNUM
00484       CALL SLABAD( SMLNUM, BIGNUM )
00485       THRSH2 = TEN*THRESH
00486       NTESTT = 0
00487       NERRS = 0
00488 *
00489 *     Go to the tests for read-in matrix pairs
00490 *
00491       IFUNC = 0
00492       IF( NSIZE.EQ.0 )
00493      $   GO TO 70
00494 *
00495 *     Test the built-in matrix pairs.
00496 *     Loop over different functions (IFUNC) of CGGESX, types (PRTYPE)
00497 *     of test matrices, different size (M+N)
00498 *
00499       PRTYPE = 0
00500       QBA = 3
00501       QBB = 4
00502       WEIGHT = SQRT( ULP )
00503 *
00504       DO 60 IFUNC = 0, 3
00505          DO 50 PRTYPE = 1, 5
00506             DO 40 M = 1, NSIZE - 1
00507                DO 30 N = 1, NSIZE - M
00508 *
00509                   WEIGHT = ONE / WEIGHT
00510                   MPLUSN = M + N
00511 *
00512 *                 Generate test matrices
00513 *
00514                   FS = .TRUE.
00515                   K = 0
00516 *
00517                   CALL CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, AI,
00518      $                         LDA )
00519                   CALL CLASET( 'Full', MPLUSN, MPLUSN, CZERO, CZERO, BI,
00520      $                         LDA )
00521 *
00522                   CALL CLATM5( PRTYPE, M, N, AI, LDA, AI( M+1, M+1 ),
00523      $                         LDA, AI( 1, M+1 ), LDA, BI, LDA,
00524      $                         BI( M+1, M+1 ), LDA, BI( 1, M+1 ), LDA,
00525      $                         Q, LDA, Z, LDA, WEIGHT, QBA, QBB )
00526 *
00527 *                 Compute the Schur factorization and swapping the
00528 *                 m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00529 *                 Swapping is accomplished via the function CLCTSX
00530 *                 which is supplied below.
00531 *
00532                   IF( IFUNC.EQ.0 ) THEN
00533                      SENSE = 'N'
00534                   ELSE IF( IFUNC.EQ.1 ) THEN
00535                      SENSE = 'E'
00536                   ELSE IF( IFUNC.EQ.2 ) THEN
00537                      SENSE = 'V'
00538                   ELSE IF( IFUNC.EQ.3 ) THEN
00539                      SENSE = 'B'
00540                   END IF
00541 *
00542                   CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00543                   CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00544 *
00545                   CALL CGGESX( 'V', 'V', 'S', CLCTSX, SENSE, MPLUSN, AI,
00546      $                         LDA, BI, LDA, MM, ALPHA, BETA, Q, LDA, Z,
00547      $                         LDA, PL, DIFEST, WORK, LWORK, RWORK,
00548      $                         IWORK, LIWORK, BWORK, LINFO )
00549 *
00550                   IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00551                      RESULT( 1 ) = ULPINV
00552                      WRITE( NOUT, FMT = 9999 )'CGGESX', LINFO, MPLUSN,
00553      $                  PRTYPE
00554                      INFO = LINFO
00555                      GO TO 30
00556                   END IF
00557 *
00558 *                 Compute the norm(A, B)
00559 *
00560                   CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK,
00561      $                         MPLUSN )
00562                   CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00563      $                         WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00564                   ABNRM = CLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN,
00565      $                    RWORK )
00566 *
00567 *                 Do tests (1) to (4)
00568 *
00569                   RESULT( 2 ) = ZERO
00570                   CALL CGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z,
00571      $                         LDA, WORK, RWORK, RESULT( 1 ) )
00572                   CALL CGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z,
00573      $                         LDA, WORK, RWORK, RESULT( 2 ) )
00574                   CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q,
00575      $                         LDA, WORK, RWORK, RESULT( 3 ) )
00576                   CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z,
00577      $                         LDA, WORK, RWORK, RESULT( 4 ) )
00578                   NTEST = 4
00579 *
00580 *                 Do tests (5) and (6): check Schur form of A and
00581 *                 compare eigenvalues with diagonals.
00582 *
00583                   TEMP1 = ZERO
00584                   RESULT( 5 ) = ZERO
00585                   RESULT( 6 ) = ZERO
00586 *
00587                   DO 10 J = 1, MPLUSN
00588                      ILABAD = .FALSE.
00589                      TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) /
00590      $                       MAX( SMLNUM, ABS1( ALPHA( J ) ),
00591      $                       ABS1( AI( J, J ) ) )+
00592      $                       ABS1( BETA( J )-BI( J, J ) ) /
00593      $                       MAX( SMLNUM, ABS1( BETA( J ) ),
00594      $                       ABS1( BI( J, J ) ) ) ) / ULP
00595                      IF( J.LT.MPLUSN ) THEN
00596                         IF( AI( J+1, J ).NE.ZERO ) THEN
00597                            ILABAD = .TRUE.
00598                            RESULT( 5 ) = ULPINV
00599                         END IF
00600                      END IF
00601                      IF( J.GT.1 ) THEN
00602                         IF( AI( J, J-1 ).NE.ZERO ) THEN
00603                            ILABAD = .TRUE.
00604                            RESULT( 5 ) = ULPINV
00605                         END IF
00606                      END IF
00607                      TEMP1 = MAX( TEMP1, TEMP2 )
00608                      IF( ILABAD ) THEN
00609                         WRITE( NOUT, FMT = 9997 )J, MPLUSN, PRTYPE
00610                      END IF
00611    10             CONTINUE
00612                   RESULT( 6 ) = TEMP1
00613                   NTEST = NTEST + 2
00614 *
00615 *                 Test (7) (if sorting worked)
00616 *
00617                   RESULT( 7 ) = ZERO
00618                   IF( LINFO.EQ.MPLUSN+3 ) THEN
00619                      RESULT( 7 ) = ULPINV
00620                   ELSE IF( MM.NE.N ) THEN
00621                      RESULT( 7 ) = ULPINV
00622                   END IF
00623                   NTEST = NTEST + 1
00624 *
00625 *                 Test (8): compare the estimated value DIF and its
00626 *                 value. first, compute the exact DIF.
00627 *
00628                   RESULT( 8 ) = ZERO
00629                   MN2 = MM*( MPLUSN-MM )*2
00630                   IF( IFUNC.GE.2 .AND. MN2.LE.NCMAX*NCMAX ) THEN
00631 *
00632 *                    Note: for either following two cases, there are
00633 *                    almost same number of test cases fail the test.
00634 *
00635                      CALL CLAKF2( MM, MPLUSN-MM, AI, LDA,
00636      $                            AI( MM+1, MM+1 ), BI,
00637      $                            BI( MM+1, MM+1 ), C, LDC )
00638 *
00639                      CALL CGESVD( 'N', 'N', MN2, MN2, C, LDC, S, WORK,
00640      $                            1, WORK( 2 ), 1, WORK( 3 ), LWORK-2,
00641      $                            RWORK, INFO )
00642                      DIFTRU = S( MN2 )
00643 *
00644                      IF( DIFEST( 2 ).EQ.ZERO ) THEN
00645                         IF( DIFTRU.GT.ABNRM*ULP )
00646      $                     RESULT( 8 ) = ULPINV
00647                      ELSE IF( DIFTRU.EQ.ZERO ) THEN
00648                         IF( DIFEST( 2 ).GT.ABNRM*ULP )
00649      $                     RESULT( 8 ) = ULPINV
00650                      ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00651      $                        ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00652                         RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ),
00653      $                                DIFEST( 2 ) / DIFTRU )
00654                      END IF
00655                      NTEST = NTEST + 1
00656                   END IF
00657 *
00658 *                 Test (9)
00659 *
00660                   RESULT( 9 ) = ZERO
00661                   IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00662                      IF( DIFTRU.GT.ABNRM*ULP )
00663      $                  RESULT( 9 ) = ULPINV
00664                      IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00665      $                  RESULT( 9 ) = ULPINV
00666                      IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00667      $                  RESULT( 9 ) = ULPINV
00668                      NTEST = NTEST + 1
00669                   END IF
00670 *
00671                   NTESTT = NTESTT + NTEST
00672 *
00673 *                 Print out tests which fail.
00674 *
00675                   DO 20 J = 1, 9
00676                      IF( RESULT( J ).GE.THRESH ) THEN
00677 *
00678 *                       If this is the first test to fail,
00679 *                       print a header to the data file.
00680 *
00681                         IF( NERRS.EQ.0 ) THEN
00682                            WRITE( NOUT, FMT = 9996 )'CGX'
00683 *
00684 *                          Matrix types
00685 *
00686                            WRITE( NOUT, FMT = 9994 )
00687 *
00688 *                          Tests performed
00689 *
00690                            WRITE( NOUT, FMT = 9993 )'unitary', '''',
00691      $                        'transpose', ( '''', I = 1, 4 )
00692 *
00693                         END IF
00694                         NERRS = NERRS + 1
00695                         IF( RESULT( J ).LT.10000.0 ) THEN
00696                            WRITE( NOUT, FMT = 9992 )MPLUSN, PRTYPE,
00697      $                        WEIGHT, M, J, RESULT( J )
00698                         ELSE
00699                            WRITE( NOUT, FMT = 9991 )MPLUSN, PRTYPE,
00700      $                        WEIGHT, M, J, RESULT( J )
00701                         END IF
00702                      END IF
00703    20             CONTINUE
00704 *
00705    30          CONTINUE
00706    40       CONTINUE
00707    50    CONTINUE
00708    60 CONTINUE
00709 *
00710       GO TO 150
00711 *
00712    70 CONTINUE
00713 *
00714 *     Read in data from file to check accuracy of condition estimation
00715 *     Read input data until N=0
00716 *
00717       NPTKNT = 0
00718 *
00719    80 CONTINUE
00720       READ( NIN, FMT = *, END = 140 )MPLUSN
00721       IF( MPLUSN.EQ.0 )
00722      $   GO TO 140
00723       READ( NIN, FMT = *, END = 140 )N
00724       DO 90 I = 1, MPLUSN
00725          READ( NIN, FMT = * )( AI( I, J ), J = 1, MPLUSN )
00726    90 CONTINUE
00727       DO 100 I = 1, MPLUSN
00728          READ( NIN, FMT = * )( BI( I, J ), J = 1, MPLUSN )
00729   100 CONTINUE
00730       READ( NIN, FMT = * )PLTRU, DIFTRU
00731 *
00732       NPTKNT = NPTKNT + 1
00733       FS = .TRUE.
00734       K = 0
00735       M = MPLUSN - N
00736 *
00737       CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, A, LDA )
00738       CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA, B, LDA )
00739 *
00740 *     Compute the Schur factorization while swaping the
00741 *     m-by-m (1,1)-blocks with n-by-n (2,2)-blocks.
00742 *
00743       CALL CGGESX( 'V', 'V', 'S', CLCTSX, 'B', MPLUSN, AI, LDA, BI, LDA,
00744      $             MM, ALPHA, BETA, Q, LDA, Z, LDA, PL, DIFEST, WORK,
00745      $             LWORK, RWORK, IWORK, LIWORK, BWORK, LINFO )
00746 *
00747       IF( LINFO.NE.0 .AND. LINFO.NE.MPLUSN+2 ) THEN
00748          RESULT( 1 ) = ULPINV
00749          WRITE( NOUT, FMT = 9998 )'CGGESX', LINFO, MPLUSN, NPTKNT
00750          GO TO 130
00751       END IF
00752 *
00753 *     Compute the norm(A, B)
00754 *        (should this be norm of (A,B) or (AI,BI)?)
00755 *
00756       CALL CLACPY( 'Full', MPLUSN, MPLUSN, AI, LDA, WORK, MPLUSN )
00757       CALL CLACPY( 'Full', MPLUSN, MPLUSN, BI, LDA,
00758      $             WORK( MPLUSN*MPLUSN+1 ), MPLUSN )
00759       ABNRM = CLANGE( 'Fro', MPLUSN, 2*MPLUSN, WORK, MPLUSN, RWORK )
00760 *
00761 *     Do tests (1) to (4)
00762 *
00763       CALL CGET51( 1, MPLUSN, A, LDA, AI, LDA, Q, LDA, Z, LDA, WORK,
00764      $             RWORK, RESULT( 1 ) )
00765       CALL CGET51( 1, MPLUSN, B, LDA, BI, LDA, Q, LDA, Z, LDA, WORK,
00766      $             RWORK, RESULT( 2 ) )
00767       CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Q, LDA, Q, LDA, WORK,
00768      $             RWORK, RESULT( 3 ) )
00769       CALL CGET51( 3, MPLUSN, B, LDA, BI, LDA, Z, LDA, Z, LDA, WORK,
00770      $             RWORK, RESULT( 4 ) )
00771 *
00772 *     Do tests (5) and (6): check Schur form of A and compare
00773 *     eigenvalues with diagonals.
00774 *
00775       NTEST = 6
00776       TEMP1 = ZERO
00777       RESULT( 5 ) = ZERO
00778       RESULT( 6 ) = ZERO
00779 *
00780       DO 110 J = 1, MPLUSN
00781          ILABAD = .FALSE.
00782          TEMP2 = ( ABS1( ALPHA( J )-AI( J, J ) ) /
00783      $           MAX( SMLNUM, ABS1( ALPHA( J ) ), ABS1( AI( J, J ) ) )+
00784      $           ABS1( BETA( J )-BI( J, J ) ) /
00785      $           MAX( SMLNUM, ABS1( BETA( J ) ), ABS1( BI( J, J ) ) ) )
00786      $            / ULP
00787          IF( J.LT.MPLUSN ) THEN
00788             IF( AI( J+1, J ).NE.ZERO ) THEN
00789                ILABAD = .TRUE.
00790                RESULT( 5 ) = ULPINV
00791             END IF
00792          END IF
00793          IF( J.GT.1 ) THEN
00794             IF( AI( J, J-1 ).NE.ZERO ) THEN
00795                ILABAD = .TRUE.
00796                RESULT( 5 ) = ULPINV
00797             END IF
00798          END IF
00799          TEMP1 = MAX( TEMP1, TEMP2 )
00800          IF( ILABAD ) THEN
00801             WRITE( NOUT, FMT = 9997 )J, MPLUSN, NPTKNT
00802          END IF
00803   110 CONTINUE
00804       RESULT( 6 ) = TEMP1
00805 *
00806 *     Test (7) (if sorting worked)  <--------- need to be checked.
00807 *
00808       NTEST = 7
00809       RESULT( 7 ) = ZERO
00810       IF( LINFO.EQ.MPLUSN+3 )
00811      $   RESULT( 7 ) = ULPINV
00812 *
00813 *     Test (8): compare the estimated value of DIF and its true value.
00814 *
00815       NTEST = 8
00816       RESULT( 8 ) = ZERO
00817       IF( DIFEST( 2 ).EQ.ZERO ) THEN
00818          IF( DIFTRU.GT.ABNRM*ULP )
00819      $      RESULT( 8 ) = ULPINV
00820       ELSE IF( DIFTRU.EQ.ZERO ) THEN
00821          IF( DIFEST( 2 ).GT.ABNRM*ULP )
00822      $      RESULT( 8 ) = ULPINV
00823       ELSE IF( ( DIFTRU.GT.THRSH2*DIFEST( 2 ) ) .OR.
00824      $         ( DIFTRU*THRSH2.LT.DIFEST( 2 ) ) ) THEN
00825          RESULT( 8 ) = MAX( DIFTRU / DIFEST( 2 ), DIFEST( 2 ) / DIFTRU )
00826       END IF
00827 *
00828 *     Test (9)
00829 *
00830       NTEST = 9
00831       RESULT( 9 ) = ZERO
00832       IF( LINFO.EQ.( MPLUSN+2 ) ) THEN
00833          IF( DIFTRU.GT.ABNRM*ULP )
00834      $      RESULT( 9 ) = ULPINV
00835          IF( ( IFUNC.GT.1 ) .AND. ( DIFEST( 2 ).NE.ZERO ) )
00836      $      RESULT( 9 ) = ULPINV
00837          IF( ( IFUNC.EQ.1 ) .AND. ( PL( 1 ).NE.ZERO ) )
00838      $      RESULT( 9 ) = ULPINV
00839       END IF
00840 *
00841 *     Test (10): compare the estimated value of PL and it true value.
00842 *
00843       NTEST = 10
00844       RESULT( 10 ) = ZERO
00845       IF( PL( 1 ).EQ.ZERO ) THEN
00846          IF( PLTRU.GT.ABNRM*ULP )
00847      $      RESULT( 10 ) = ULPINV
00848       ELSE IF( PLTRU.EQ.ZERO ) THEN
00849          IF( PL( 1 ).GT.ABNRM*ULP )
00850      $      RESULT( 10 ) = ULPINV
00851       ELSE IF( ( PLTRU.GT.THRESH*PL( 1 ) ) .OR.
00852      $         ( PLTRU*THRESH.LT.PL( 1 ) ) ) THEN
00853          RESULT( 10 ) = ULPINV
00854       END IF
00855 *
00856       NTESTT = NTESTT + NTEST
00857 *
00858 *     Print out tests which fail.
00859 *
00860       DO 120 J = 1, NTEST
00861          IF( RESULT( J ).GE.THRESH ) THEN
00862 *
00863 *           If this is the first test to fail,
00864 *           print a header to the data file.
00865 *
00866             IF( NERRS.EQ.0 ) THEN
00867                WRITE( NOUT, FMT = 9996 )'CGX'
00868 *
00869 *              Matrix types
00870 *
00871                WRITE( NOUT, FMT = 9995 )
00872 *
00873 *              Tests performed
00874 *
00875                WRITE( NOUT, FMT = 9993 )'unitary', '''', 'transpose',
00876      $            ( '''', I = 1, 4 )
00877 *
00878             END IF
00879             NERRS = NERRS + 1
00880             IF( RESULT( J ).LT.10000.0 ) THEN
00881                WRITE( NOUT, FMT = 9990 )NPTKNT, MPLUSN, J, RESULT( J )
00882             ELSE
00883                WRITE( NOUT, FMT = 9989 )NPTKNT, MPLUSN, J, RESULT( J )
00884             END IF
00885          END IF
00886 *
00887   120 CONTINUE
00888 *
00889   130 CONTINUE
00890       GO TO 80
00891   140 CONTINUE
00892 *
00893   150 CONTINUE
00894 *
00895 *     Summary
00896 *
00897       CALL ALASVM( 'CGX', NOUT, NERRS, NTESTT, 0 )
00898 *
00899       WORK( 1 ) = MAXWRK
00900 *
00901       RETURN
00902 *
00903  9999 FORMAT( ' CDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00904      $      I6, ', JTYPE=', I6, ')' )
00905 *
00906  9998 FORMAT( ' CDRGSX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00907      $      I6, ', Input Example #', I2, ')' )
00908 *
00909  9997 FORMAT( ' CDRGSX: S not in Schur form at eigenvalue ', I6, '.',
00910      $      / 9X, 'N=', I6, ', JTYPE=', I6, ')' )
00911 *
00912  9996 FORMAT( / 1X, A3, ' -- Complex Expert Generalized Schur form',
00913      $      ' problem driver' )
00914 *
00915  9995 FORMAT( 'Input Example' )
00916 *
00917  9994 FORMAT( ' Matrix types: ', /
00918      $      '  1:  A is a block diagonal matrix of Jordan blocks ',
00919      $      'and B is the identity ', / '      matrix, ',
00920      $      / '  2:  A and B are upper triangular matrices, ',
00921      $      / '  3:  A and B are as type 2, but each second diagonal ',
00922      $      'block in A_11 and ', /
00923      $      '      each third diaongal block in A_22 are 2x2 blocks,',
00924      $      / '  4:  A and B are block diagonal matrices, ',
00925      $      / '  5:  (A,B) has potentially close or common ',
00926      $      'eigenvalues.', / )
00927 *
00928  9993 FORMAT( / ' Tests performed:  (S is Schur, T is triangular, ',
00929      $      'Q and Z are ', A, ',', / 19X,
00930      $      ' a is alpha, b is beta, and ', A, ' means ', A, '.)',
00931      $      / '  1 = | A - Q S Z', A,
00932      $      ' | / ( |A| n ulp )      2 = | B - Q T Z', A,
00933      $      ' | / ( |B| n ulp )', / '  3 = | I - QQ', A,
00934      $      ' | / ( n ulp )             4 = | I - ZZ', A,
00935      $      ' | / ( n ulp )', / '  5 = 1/ULP  if A is not in ',
00936      $      'Schur form S', / '  6 = difference between (alpha,beta)',
00937      $      ' and diagonals of (S,T)', /
00938      $      '  7 = 1/ULP  if SDIM is not the correct number of ',
00939      $      'selected eigenvalues', /
00940      $      '  8 = 1/ULP  if DIFEST/DIFTRU > 10*THRESH or ',
00941      $      'DIFTRU/DIFEST > 10*THRESH',
00942      $      / '  9 = 1/ULP  if DIFEST <> 0 or DIFTRU > ULP*norm(A,B) ',
00943      $      'when reordering fails', /
00944      $      ' 10 = 1/ULP  if PLEST/PLTRU > THRESH or ',
00945      $      'PLTRU/PLEST > THRESH', /
00946      $      '    ( Test 10 is only for input examples )', / )
00947  9992 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.3,
00948      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, F8.2 )
00949  9991 FORMAT( ' Matrix order=', I2, ', type=', I2, ', a=', E10.3,
00950      $      ', order(A_11)=', I2, ', result ', I2, ' is ', 0P, E10.3 )
00951  9990 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00952      $      ' result ', I2, ' is', 0P, F8.2 )
00953  9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00954      $      ' result ', I2, ' is', 1P, E10.3 )
00955 *
00956 *     End of CDRGSX
00957 *
00958       END
 All Files Functions