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