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