LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cdrgvx.f
Go to the documentation of this file.
00001 *> \brief \b CDRGVX
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 CDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
00012 *                          ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
00013 *                          S, STRU, DIF, DIFTRU, WORK, LWORK, RWORK,
00014 *                          IWORK, LIWORK, RESULT, BWORK, INFO )
00015 * 
00016 *       .. Scalar Arguments ..
00017 *       INTEGER            IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
00018 *      $                   NSIZE
00019 *       REAL               THRESH
00020 *       ..
00021 *       .. Array Arguments ..
00022 *       LOGICAL            BWORK( * )
00023 *       INTEGER            IWORK( * )
00024 *       REAL               DIF( * ), DIFTRU( * ), LSCALE( * ),
00025 *      $                   RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * ),
00026 *      $                   STRU( * )
00027 *       COMPLEX            A( LDA, * ), AI( LDA, * ), ALPHA( * ),
00028 *      $                   B( LDA, * ), BETA( * ), BI( LDA, * ),
00029 *      $                   VL( LDA, * ), VR( LDA, * ), WORK( * )
00030 *       ..
00031 *  
00032 *
00033 *> \par Purpose:
00034 *  =============
00035 *>
00036 *> \verbatim
00037 *>
00038 *> CDRGVX checks the nonsymmetric generalized eigenvalue problem
00039 *> expert driver CGGEVX.
00040 *>
00041 *> CGGEVX computes the generalized eigenvalues, (optionally) the left
00042 *> and/or right eigenvectors, (optionally) computes a balancing
00043 *> transformation to improve the conditioning, and (optionally)
00044 *> reciprocal condition numbers for the eigenvalues and eigenvectors.
00045 *>
00046 *> When CDRGVX is called with NSIZE > 0, two types of test matrix pairs
00047 *> are generated by the subroutine SLATM6 and test the driver CGGEVX.
00048 *> The test matrices have the known exact condition numbers for
00049 *> eigenvalues. For the condition numbers of the eigenvectors
00050 *> corresponding the first and last eigenvalues are also know
00051 *> ``exactly'' (see CLATM6).
00052 *> For each matrix pair, the following tests will be performed and
00053 *> compared with the threshhold THRESH.
00054 *>
00055 *> (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
00056 *>
00057 *>    | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
00058 *>
00059 *>     where l**H is the conjugate tranpose of l.
00060 *>
00061 *> (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
00062 *>
00063 *>       | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
00064 *>
00065 *> (3) The condition number S(i) of eigenvalues computed by CGGEVX
00066 *>     differs less than a factor THRESH from the exact S(i) (see
00067 *>     CLATM6).
00068 *>
00069 *> (4) DIF(i) computed by CTGSNA differs less than a factor 10*THRESH
00070 *>     from the exact value (for the 1st and 5th vectors only).
00071 *>
00072 *> Test Matrices
00073 *> =============
00074 *>
00075 *> Two kinds of test matrix pairs
00076 *>          (A, B) = inverse(YH) * (Da, Db) * inverse(X)
00077 *> are used in the tests:
00078 *>
00079 *> 1: Da = 1+a   0    0    0    0    Db = 1   0   0   0   0
00080 *>          0   2+a   0    0    0         0   1   0   0   0
00081 *>          0    0   3+a   0    0         0   0   1   0   0
00082 *>          0    0    0   4+a   0         0   0   0   1   0
00083 *>          0    0    0    0   5+a ,      0   0   0   0   1 , and
00084 *>
00085 *> 2: Da =  1   -1    0    0    0    Db = 1   0   0   0   0
00086 *>          1    1    0    0    0         0   1   0   0   0
00087 *>          0    0    1    0    0         0   0   1   0   0
00088 *>          0    0    0   1+a  1+b        0   0   0   1   0
00089 *>          0    0    0  -1-b  1+a ,      0   0   0   0   1 .
00090 *>
00091 *> In both cases the same inverse(YH) and inverse(X) are used to compute
00092 *> (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
00093 *>
00094 *> YH:  =  1    0   -y    y   -y    X =  1   0  -x  -x   x
00095 *>         0    1   -y    y   -y         0   1   x  -x  -x
00096 *>         0    0    1    0    0         0   0   1   0   0
00097 *>         0    0    0    1    0         0   0   0   1   0
00098 *>         0    0    0    0    1,        0   0   0   0   1 , where
00099 *>
00100 *> a, b, x and y will have all values independently of each other from
00101 *> { sqrt(sqrt(ULP)),  0.1,  1,  10,  1/sqrt(sqrt(ULP)) }.
00102 *> \endverbatim
00103 *
00104 *  Arguments:
00105 *  ==========
00106 *
00107 *> \param[in] NSIZE
00108 *> \verbatim
00109 *>          NSIZE is INTEGER
00110 *>          The number of sizes of matrices to use.  NSIZE must be at
00111 *>          least zero. If it is zero, no randomly generated matrices
00112 *>          are tested, but any test matrices read from NIN will be
00113 *>          tested.  If it is not zero, then N = 5.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] THRESH
00117 *> \verbatim
00118 *>          THRESH is REAL
00119 *>          A test will count as "failed" if the "error", computed as
00120 *>          described above, exceeds THRESH.  Note that the error
00121 *>          is scaled to be O(1), so THRESH should be a reasonably
00122 *>          small multiple of 1, e.g., 10 or 100.  In particular,
00123 *>          it should not depend on the precision (single vs. double)
00124 *>          or the size of the matrix.  It must be at least zero.
00125 *> \endverbatim
00126 *>
00127 *> \param[in] NIN
00128 *> \verbatim
00129 *>          NIN is INTEGER
00130 *>          The FORTRAN unit number for reading in the data file of
00131 *>          problems to solve.
00132 *> \endverbatim
00133 *>
00134 *> \param[in] NOUT
00135 *> \verbatim
00136 *>          NOUT is INTEGER
00137 *>          The FORTRAN unit number for printing out error messages
00138 *>          (e.g., if a routine returns IINFO not equal to 0.)
00139 *> \endverbatim
00140 *>
00141 *> \param[out] A
00142 *> \verbatim
00143 *>          A is COMPLEX array, dimension (LDA, NSIZE)
00144 *>          Used to hold the matrix whose eigenvalues are to be
00145 *>          computed.  On exit, A contains the last matrix actually used.
00146 *> \endverbatim
00147 *>
00148 *> \param[in] LDA
00149 *> \verbatim
00150 *>          LDA is INTEGER
00151 *>          The leading dimension of A, B, AI, BI, Ao, and Bo.
00152 *>          It must be at least 1 and at least NSIZE.
00153 *> \endverbatim
00154 *>
00155 *> \param[out] B
00156 *> \verbatim
00157 *>          B is COMPLEX array, dimension (LDA, NSIZE)
00158 *>          Used to hold the matrix whose eigenvalues are to be
00159 *>          computed.  On exit, B contains the last matrix actually used.
00160 *> \endverbatim
00161 *>
00162 *> \param[out] AI
00163 *> \verbatim
00164 *>          AI is COMPLEX array, dimension (LDA, NSIZE)
00165 *>          Copy of A, modified by CGGEVX.
00166 *> \endverbatim
00167 *>
00168 *> \param[out] BI
00169 *> \verbatim
00170 *>          BI is COMPLEX array, dimension (LDA, NSIZE)
00171 *>          Copy of B, modified by CGGEVX.
00172 *> \endverbatim
00173 *>
00174 *> \param[out] ALPHA
00175 *> \verbatim
00176 *>          ALPHA is COMPLEX array, dimension (NSIZE)
00177 *> \endverbatim
00178 *>
00179 *> \param[out] BETA
00180 *> \verbatim
00181 *>          BETA is COMPLEX array, dimension (NSIZE)
00182 *>
00183 *>          On exit, ALPHA/BETA are the eigenvalues.
00184 *> \endverbatim
00185 *>
00186 *> \param[out] VL
00187 *> \verbatim
00188 *>          VL is COMPLEX array, dimension (LDA, NSIZE)
00189 *>          VL holds the left eigenvectors computed by CGGEVX.
00190 *> \endverbatim
00191 *>
00192 *> \param[out] VR
00193 *> \verbatim
00194 *>          VR is COMPLEX array, dimension (LDA, NSIZE)
00195 *>          VR holds the right eigenvectors computed by CGGEVX.
00196 *> \endverbatim
00197 *>
00198 *> \param[out] ILO
00199 *> \verbatim
00200 *>              ILO is INTEGER
00201 *> \endverbatim
00202 *>
00203 *> \param[out] IHI
00204 *> \verbatim
00205 *>              IHI is INTEGER
00206 *> \endverbatim
00207 *>
00208 *> \param[out] LSCALE   
00209 *> \verbatim
00210 *>              LSCALE is REAL array, dimension (N)
00211 *> \endverbatim
00212 *>
00213 *> \param[out] RSCALE   
00214 *> \verbatim
00215 *>              RSCALE is REAL array, dimension (N)
00216 *> \endverbatim
00217 *>
00218 *> \param[out] S        
00219 *> \verbatim
00220 *>              S is REAL array, dimension (N)
00221 *> \endverbatim
00222 *>
00223 *> \param[out] STRU     
00224 *> \verbatim
00225 *>              STRU is REAL array, dimension (N)
00226 *> \endverbatim
00227 *>
00228 *> \param[out] DIF              
00229 *> \verbatim
00230 *>              DIF is REAL array, dimension (N)
00231 *> \endverbatim
00232 *>
00233 *> \param[out] DIFTRU           
00234 *> \verbatim
00235 *>              DIFTRU is REAL array, dimension (N)
00236 *> \endverbatim
00237 *>
00238 *> \param[out] WORK
00239 *> \verbatim
00240 *>          WORK is COMPLEX array, dimension (LWORK)
00241 *> \endverbatim
00242 *>
00243 *> \param[in] LWORK
00244 *> \verbatim
00245 *>          LWORK is INTEGER
00246 *>          Leading dimension of WORK.  LWORK >= 2*N*N + 2*N
00247 *> \endverbatim
00248 *>
00249 *> \param[out] RWORK
00250 *> \verbatim
00251 *>          RWORK is REAL array, dimension (6*N)
00252 *> \endverbatim
00253 *>
00254 *> \param[out] IWORK
00255 *> \verbatim
00256 *>          IWORK is INTEGER array, dimension (LIWORK)
00257 *> \endverbatim
00258 *>
00259 *> \param[in] LIWORK
00260 *> \verbatim
00261 *>          LIWORK is INTEGER
00262 *>          Leading dimension of IWORK.  LIWORK >= N+2.
00263 *> \endverbatim
00264 *>
00265 *> \param[out] RESULT
00266 *> \verbatim
00267 *>              RESULT is REAL array, dimension (4)
00268 *> \endverbatim
00269 *>
00270 *> \param[out] BWORK
00271 *> \verbatim
00272 *>          BWORK is LOGICAL array, dimension (N)
00273 *> \endverbatim
00274 *>
00275 *> \param[out] INFO
00276 *> \verbatim
00277 *>          INFO is INTEGER
00278 *>          = 0:  successful exit
00279 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00280 *>          > 0:  A routine returned an error code.
00281 *> \endverbatim
00282 *
00283 *  Authors:
00284 *  ========
00285 *
00286 *> \author Univ. of Tennessee 
00287 *> \author Univ. of California Berkeley 
00288 *> \author Univ. of Colorado Denver 
00289 *> \author NAG Ltd. 
00290 *
00291 *> \date November 2011
00292 *
00293 *> \ingroup complex_eig
00294 *
00295 *  =====================================================================
00296       SUBROUTINE CDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
00297      $                   ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
00298      $                   S, STRU, DIF, DIFTRU, WORK, LWORK, RWORK,
00299      $                   IWORK, LIWORK, RESULT, BWORK, INFO )
00300 *
00301 *  -- LAPACK test routine (version 3.4.0) --
00302 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00303 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00304 *     November 2011
00305 *
00306 *     .. Scalar Arguments ..
00307       INTEGER            IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
00308      $                   NSIZE
00309       REAL               THRESH
00310 *     ..
00311 *     .. Array Arguments ..
00312       LOGICAL            BWORK( * )
00313       INTEGER            IWORK( * )
00314       REAL               DIF( * ), DIFTRU( * ), LSCALE( * ),
00315      $                   RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * ),
00316      $                   STRU( * )
00317       COMPLEX            A( LDA, * ), AI( LDA, * ), ALPHA( * ),
00318      $                   B( LDA, * ), BETA( * ), BI( LDA, * ),
00319      $                   VL( LDA, * ), VR( LDA, * ), WORK( * )
00320 *     ..
00321 *
00322 *  =====================================================================
00323 *
00324 *     .. Parameters ..
00325       REAL               ZERO, ONE, TEN, TNTH, HALF
00326       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1,
00327      $                   TNTH = 1.0E-1, HALF = 0.5E+0 )
00328 *     ..
00329 *     .. Local Scalars ..
00330       INTEGER            I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
00331      $                   MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
00332       REAL               ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
00333      $                   ULP, ULPINV
00334 *     ..
00335 *     .. Local Arrays ..
00336       COMPLEX            WEIGHT( 5 )
00337 *     ..
00338 *     .. External Functions ..
00339       INTEGER            ILAENV
00340       REAL               CLANGE, SLAMCH
00341       EXTERNAL           ILAENV, CLANGE, SLAMCH
00342 *     ..
00343 *     .. External Subroutines ..
00344       EXTERNAL           ALASVM, CGET52, CGGEVX, CLACPY, CLATM6, XERBLA
00345 *     ..
00346 *     .. Intrinsic Functions ..
00347       INTRINSIC          ABS, CMPLX, MAX, SQRT
00348 *     ..
00349 *     .. Executable Statements ..
00350 *
00351 *     Check for errors
00352 *
00353       INFO = 0
00354 *
00355       NMAX = 5
00356 *
00357       IF( NSIZE.LT.0 ) THEN
00358          INFO = -1
00359       ELSE IF( THRESH.LT.ZERO ) THEN
00360          INFO = -2
00361       ELSE IF( NIN.LE.0 ) THEN
00362          INFO = -3
00363       ELSE IF( NOUT.LE.0 ) THEN
00364          INFO = -4
00365       ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
00366          INFO = -6
00367       ELSE IF( LIWORK.LT.NMAX+2 ) THEN
00368          INFO = -26
00369       END IF
00370 *
00371 *     Compute workspace
00372 *      (Note: Comments in the code beginning "Workspace:" describe the
00373 *       minimal amount of workspace needed at that point in the code,
00374 *       as well as the preferred amount for good performance.
00375 *       NB refers to the optimal block size for the immediately
00376 *       following subroutine, as returned by ILAENV.)
00377 *
00378       MINWRK = 1
00379       IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
00380          MINWRK = 2*NMAX*( NMAX+1 )
00381          MAXWRK = NMAX*( 1+ILAENV( 1, 'CGEQRF', ' ', NMAX, 1, NMAX,
00382      $            0 ) )
00383          MAXWRK = MAX( MAXWRK, 2*NMAX*( NMAX+1 ) )
00384          WORK( 1 ) = MAXWRK
00385       END IF
00386 *
00387       IF( LWORK.LT.MINWRK )
00388      $   INFO = -23
00389 *
00390       IF( INFO.NE.0 ) THEN
00391          CALL XERBLA( 'CDRGVX', -INFO )
00392          RETURN
00393       END IF
00394 *
00395       N = 5
00396       ULP = SLAMCH( 'P' )
00397       ULPINV = ONE / ULP
00398       THRSH2 = TEN*THRESH
00399       NERRS = 0
00400       NPTKNT = 0
00401       NTESTT = 0
00402 *
00403       IF( NSIZE.EQ.0 )
00404      $   GO TO 90
00405 *
00406 *     Parameters used for generating test matrices.
00407 *
00408       WEIGHT( 1 ) = CMPLX( TNTH, ZERO )
00409       WEIGHT( 2 ) = CMPLX( HALF, ZERO )
00410       WEIGHT( 3 ) = ONE
00411       WEIGHT( 4 ) = ONE / WEIGHT( 2 )
00412       WEIGHT( 5 ) = ONE / WEIGHT( 1 )
00413 *
00414       DO 80 IPTYPE = 1, 2
00415          DO 70 IWA = 1, 5
00416             DO 60 IWB = 1, 5
00417                DO 50 IWX = 1, 5
00418                   DO 40 IWY = 1, 5
00419 *
00420 *                    generated a pair of test matrix
00421 *
00422                      CALL CLATM6( IPTYPE, 5, A, LDA, B, VR, LDA, VL,
00423      $                            LDA, WEIGHT( IWA ), WEIGHT( IWB ),
00424      $                            WEIGHT( IWX ), WEIGHT( IWY ), STRU,
00425      $                            DIFTRU )
00426 *
00427 *                    Compute eigenvalues/eigenvectors of (A, B).
00428 *                    Compute eigenvalue/eigenvector condition numbers
00429 *                    using computed eigenvectors.
00430 *
00431                      CALL CLACPY( 'F', N, N, A, LDA, AI, LDA )
00432                      CALL CLACPY( 'F', N, N, B, LDA, BI, LDA )
00433 *
00434                      CALL CGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI,
00435      $                            LDA, ALPHA, BETA, VL, LDA, VR, LDA,
00436      $                            ILO, IHI, LSCALE, RSCALE, ANORM,
00437      $                            BNORM, S, DIF, WORK, LWORK, RWORK,
00438      $                            IWORK, BWORK, LINFO )
00439                      IF( LINFO.NE.0 ) THEN
00440                         WRITE( NOUT, FMT = 9999 )'CGGEVX', LINFO, N,
00441      $                     IPTYPE, IWA, IWB, IWX, IWY
00442                         GO TO 30
00443                      END IF
00444 *
00445 *                    Compute the norm(A, B)
00446 *
00447                      CALL CLACPY( 'Full', N, N, AI, LDA, WORK, N )
00448                      CALL CLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ),
00449      $                            N )
00450                      ABNORM = CLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
00451 *
00452 *                    Tests (1) and (2)
00453 *
00454                      RESULT( 1 ) = ZERO
00455                      CALL CGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA,
00456      $                            ALPHA, BETA, WORK, RWORK,
00457      $                            RESULT( 1 ) )
00458                      IF( RESULT( 2 ).GT.THRESH ) THEN
00459                         WRITE( NOUT, FMT = 9998 )'Left', 'CGGEVX',
00460      $                     RESULT( 2 ), N, IPTYPE, IWA, IWB, IWX, IWY
00461                      END IF
00462 *
00463                      RESULT( 2 ) = ZERO
00464                      CALL CGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA,
00465      $                            ALPHA, BETA, WORK, RWORK,
00466      $                            RESULT( 2 ) )
00467                      IF( RESULT( 3 ).GT.THRESH ) THEN
00468                         WRITE( NOUT, FMT = 9998 )'Right', 'CGGEVX',
00469      $                     RESULT( 3 ), N, IPTYPE, IWA, IWB, IWX, IWY
00470                      END IF
00471 *
00472 *                    Test (3)
00473 *
00474                      RESULT( 3 ) = ZERO
00475                      DO 10 I = 1, N
00476                         IF( S( I ).EQ.ZERO ) THEN
00477                            IF( STRU( I ).GT.ABNORM*ULP )
00478      $                        RESULT( 3 ) = ULPINV
00479                         ELSE IF( STRU( I ).EQ.ZERO ) THEN
00480                            IF( S( I ).GT.ABNORM*ULP )
00481      $                        RESULT( 3 ) = ULPINV
00482                         ELSE
00483                            RWORK( I ) = MAX( ABS( STRU( I ) / S( I ) ),
00484      $                                  ABS( S( I ) / STRU( I ) ) )
00485                            RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
00486                         END IF
00487    10                CONTINUE
00488 *
00489 *                    Test (4)
00490 *
00491                      RESULT( 4 ) = ZERO
00492                      IF( DIF( 1 ).EQ.ZERO ) THEN
00493                         IF( DIFTRU( 1 ).GT.ABNORM*ULP )
00494      $                     RESULT( 4 ) = ULPINV
00495                      ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
00496                         IF( DIF( 1 ).GT.ABNORM*ULP )
00497      $                     RESULT( 4 ) = ULPINV
00498                      ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
00499                         IF( DIFTRU( 5 ).GT.ABNORM*ULP )
00500      $                     RESULT( 4 ) = ULPINV
00501                      ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
00502                         IF( DIF( 5 ).GT.ABNORM*ULP )
00503      $                     RESULT( 4 ) = ULPINV
00504                      ELSE
00505                         RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
00506      $                           ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
00507                         RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
00508      $                           ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
00509                         RESULT( 4 ) = MAX( RATIO1, RATIO2 )
00510                      END IF
00511 *
00512                      NTESTT = NTESTT + 4
00513 *
00514 *                    Print out tests which fail.
00515 *
00516                      DO 20 J = 1, 4
00517                         IF( ( RESULT( J ).GE.THRSH2 .AND. J.GE.4 ) .OR.
00518      $                      ( RESULT( J ).GE.THRESH .AND. J.LE.3 ) )
00519      $                       THEN
00520 *
00521 *                       If this is the first test to fail,
00522 *                       print a header to the data file.
00523 *
00524                            IF( NERRS.EQ.0 ) THEN
00525                               WRITE( NOUT, FMT = 9997 )'CXV'
00526 *
00527 *                          Print out messages for built-in examples
00528 *
00529 *                          Matrix types
00530 *
00531                               WRITE( NOUT, FMT = 9995 )
00532                               WRITE( NOUT, FMT = 9994 )
00533                               WRITE( NOUT, FMT = 9993 )
00534 *
00535 *                          Tests performed
00536 *
00537                               WRITE( NOUT, FMT = 9992 )'''',
00538      $                           'transpose', ''''
00539 *
00540                            END IF
00541                            NERRS = NERRS + 1
00542                            IF( RESULT( J ).LT.10000.0 ) THEN
00543                               WRITE( NOUT, FMT = 9991 )IPTYPE, IWA,
00544      $                           IWB, IWX, IWY, J, RESULT( J )
00545                            ELSE
00546                               WRITE( NOUT, FMT = 9990 )IPTYPE, IWA,
00547      $                           IWB, IWX, IWY, J, RESULT( J )
00548                            END IF
00549                         END IF
00550    20                CONTINUE
00551 *
00552    30                CONTINUE
00553 *
00554    40             CONTINUE
00555    50          CONTINUE
00556    60       CONTINUE
00557    70    CONTINUE
00558    80 CONTINUE
00559 *
00560       GO TO 150
00561 *
00562    90 CONTINUE
00563 *
00564 *     Read in data from file to check accuracy of condition estimation
00565 *     Read input data until N=0
00566 *
00567       READ( NIN, FMT = *, END = 150 )N
00568       IF( N.EQ.0 )
00569      $   GO TO 150
00570       DO 100 I = 1, N
00571          READ( NIN, FMT = * )( A( I, J ), J = 1, N )
00572   100 CONTINUE
00573       DO 110 I = 1, N
00574          READ( NIN, FMT = * )( B( I, J ), J = 1, N )
00575   110 CONTINUE
00576       READ( NIN, FMT = * )( STRU( I ), I = 1, N )
00577       READ( NIN, FMT = * )( DIFTRU( I ), I = 1, N )
00578 *
00579       NPTKNT = NPTKNT + 1
00580 *
00581 *     Compute eigenvalues/eigenvectors of (A, B).
00582 *     Compute eigenvalue/eigenvector condition numbers
00583 *     using computed eigenvectors.
00584 *
00585       CALL CLACPY( 'F', N, N, A, LDA, AI, LDA )
00586       CALL CLACPY( 'F', N, N, B, LDA, BI, LDA )
00587 *
00588       CALL CGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI, LDA, ALPHA, BETA,
00589      $             VL, LDA, VR, LDA, ILO, IHI, LSCALE, RSCALE, ANORM,
00590      $             BNORM, S, DIF, WORK, LWORK, RWORK, IWORK, BWORK,
00591      $             LINFO )
00592 *
00593       IF( LINFO.NE.0 ) THEN
00594          WRITE( NOUT, FMT = 9987 )'CGGEVX', LINFO, N, NPTKNT
00595          GO TO 140
00596       END IF
00597 *
00598 *     Compute the norm(A, B)
00599 *
00600       CALL CLACPY( 'Full', N, N, AI, LDA, WORK, N )
00601       CALL CLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ), N )
00602       ABNORM = CLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
00603 *
00604 *     Tests (1) and (2)
00605 *
00606       RESULT( 1 ) = ZERO
00607       CALL CGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, ALPHA, BETA,
00608      $             WORK, RWORK, RESULT( 1 ) )
00609       IF( RESULT( 2 ).GT.THRESH ) THEN
00610          WRITE( NOUT, FMT = 9986 )'Left', 'CGGEVX', RESULT( 2 ), N,
00611      $      NPTKNT
00612       END IF
00613 *
00614       RESULT( 2 ) = ZERO
00615       CALL CGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, ALPHA, BETA,
00616      $             WORK, RWORK, RESULT( 2 ) )
00617       IF( RESULT( 3 ).GT.THRESH ) THEN
00618          WRITE( NOUT, FMT = 9986 )'Right', 'CGGEVX', RESULT( 3 ), N,
00619      $      NPTKNT
00620       END IF
00621 *
00622 *     Test (3)
00623 *
00624       RESULT( 3 ) = ZERO
00625       DO 120 I = 1, N
00626          IF( S( I ).EQ.ZERO ) THEN
00627             IF( STRU( I ).GT.ABNORM*ULP )
00628      $         RESULT( 3 ) = ULPINV
00629          ELSE IF( STRU( I ).EQ.ZERO ) THEN
00630             IF( S( I ).GT.ABNORM*ULP )
00631      $         RESULT( 3 ) = ULPINV
00632          ELSE
00633             RWORK( I ) = MAX( ABS( STRU( I ) / S( I ) ),
00634      $                   ABS( S( I ) / STRU( I ) ) )
00635             RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
00636          END IF
00637   120 CONTINUE
00638 *
00639 *     Test (4)
00640 *
00641       RESULT( 4 ) = ZERO
00642       IF( DIF( 1 ).EQ.ZERO ) THEN
00643          IF( DIFTRU( 1 ).GT.ABNORM*ULP )
00644      $      RESULT( 4 ) = ULPINV
00645       ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
00646          IF( DIF( 1 ).GT.ABNORM*ULP )
00647      $      RESULT( 4 ) = ULPINV
00648       ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
00649          IF( DIFTRU( 5 ).GT.ABNORM*ULP )
00650      $      RESULT( 4 ) = ULPINV
00651       ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
00652          IF( DIF( 5 ).GT.ABNORM*ULP )
00653      $      RESULT( 4 ) = ULPINV
00654       ELSE
00655          RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
00656      $            ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
00657          RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
00658      $            ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
00659          RESULT( 4 ) = MAX( RATIO1, RATIO2 )
00660       END IF
00661 *
00662       NTESTT = NTESTT + 4
00663 *
00664 *     Print out tests which fail.
00665 *
00666       DO 130 J = 1, 4
00667          IF( RESULT( J ).GE.THRSH2 ) THEN
00668 *
00669 *           If this is the first test to fail,
00670 *           print a header to the data file.
00671 *
00672             IF( NERRS.EQ.0 ) THEN
00673                WRITE( NOUT, FMT = 9997 )'CXV'
00674 *
00675 *              Print out messages for built-in examples
00676 *
00677 *              Matrix types
00678 *
00679                WRITE( NOUT, FMT = 9996 )
00680 *
00681 *              Tests performed
00682 *
00683                WRITE( NOUT, FMT = 9992 )'''', 'transpose', ''''
00684 *
00685             END IF
00686             NERRS = NERRS + 1
00687             IF( RESULT( J ).LT.10000.0 ) THEN
00688                WRITE( NOUT, FMT = 9989 )NPTKNT, N, J, RESULT( J )
00689             ELSE
00690                WRITE( NOUT, FMT = 9988 )NPTKNT, N, J, RESULT( J )
00691             END IF
00692          END IF
00693   130 CONTINUE
00694 *
00695   140 CONTINUE
00696 *
00697       GO TO 90
00698   150 CONTINUE
00699 *
00700 *     Summary
00701 *
00702       CALL ALASVM( 'CXV', NOUT, NERRS, NTESTT, 0 )
00703 *
00704       WORK( 1 ) = MAXWRK
00705 *
00706       RETURN
00707 *
00708  9999 FORMAT( ' CDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00709      $      I6, ', JTYPE=', I6, ')' )
00710 *
00711  9998 FORMAT( ' CDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
00712      $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
00713      $      'N=', I6, ', JTYPE=', I6, ', IWA=', I5, ', IWB=', I5,
00714      $      ', IWX=', I5, ', IWY=', I5 )
00715 *
00716  9997 FORMAT( / 1X, A3, ' -- Complex Expert Eigenvalue/vector',
00717      $      ' problem driver' )
00718 *
00719  9996 FORMAT( 'Input Example' )
00720 *
00721  9995 FORMAT( ' Matrix types: ', / )
00722 *
00723  9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ',
00724      $      / '     A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
00725      $      / '     YH and X are left and right eigenvectors. ', / )
00726 *
00727  9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ',
00728      $      / '     A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
00729      $      / '     YH and X are left and right eigenvectors. ', / )
00730 *
00731  9992 FORMAT( / ' Tests performed:  ', / 4X,
00732      $      ' a is alpha, b is beta, l is a left eigenvector, ', / 4X,
00733      $      ' r is a right eigenvector and ', A, ' means ', A, '.',
00734      $      / ' 1 = max | ( b A - a B )', A, ' l | / const.',
00735      $      / ' 2 = max | ( b A - a B ) r | / const.',
00736      $      / ' 3 = max ( Sest/Stru, Stru/Sest ) ',
00737      $      ' over all eigenvalues', /
00738      $      ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
00739      $      ' over the 1st and 5th eigenvectors', / )
00740 *
00741  9991 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
00742      $      I2, ', IWY=', I2, ', result ', I2, ' is', 0P, F8.2 )
00743 *
00744  9990 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
00745      $      I2, ', IWY=', I2, ', result ', I2, ' is', 1P, E10.3 )
00746 *
00747  9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00748      $      ' result ', I2, ' is', 0P, F8.2 )
00749 *
00750  9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
00751      $      ' result ', I2, ' is', 1P, E10.3 )
00752 *
00753  9987 FORMAT( ' CDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
00754      $      I6, ', Input example #', I2, ')' )
00755 *
00756  9986 FORMAT( ' CDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
00757      $      'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
00758      $      'N=', I6, ', Input Example #', I2, ')' )
00759 *
00760 *     End of CDRGVX
00761 *
00762       END
 All Files Functions