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