LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
sget52.f
Go to the documentation of this file.
00001 *> \brief \b SGET52
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 SGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
00012 *                          ALPHAI, BETA, WORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       LOGICAL            LEFT
00016 *       INTEGER            LDA, LDB, LDE, N
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00020 *      $                   B( LDB, * ), BETA( * ), E( LDE, * ),
00021 *      $                   RESULT( 2 ), WORK( * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> SGET52  does an eigenvector check for the generalized eigenvalue
00031 *> problem.
00032 *>
00033 *> The basic test for right eigenvectors is:
00034 *>
00035 *>                           | b(j) A E(j) -  a(j) B E(j) |
00036 *>         RESULT(1) = max   -------------------------------
00037 *>                      j    n ulp max( |b(j) A|, |a(j) B| )
00038 *>
00039 *> using the 1-norm.  Here, a(j)/b(j) = w is the j-th generalized
00040 *> eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th
00041 *> generalized eigenvalue of m A - B.
00042 *>
00043 *> For real eigenvalues, the test is straightforward.  For complex
00044 *> eigenvalues, E(j) and a(j) are complex, represented by
00045 *> Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that
00046 *> eigenvector becomes
00047 *>
00048 *>                 max( |Wr|, |Wi| )
00049 *>     --------------------------------------------
00050 *>     n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| )
00051 *>
00052 *> where
00053 *>
00054 *>     Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j)
00055 *>
00056 *>     Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j)
00057 *>
00058 *>                         T   T  _
00059 *> For left eigenvectors, A , B , a, and b  are used.
00060 *>
00061 *> SGET52 also tests the normalization of E.  Each eigenvector is
00062 *> supposed to be normalized so that the maximum "absolute value"
00063 *> of its elements is 1, where in this case, "absolute value"
00064 *> of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
00065 *> maximum "absolute value" norm of a vector v  M(v). 
00066 *> if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate
00067 *> vector.  The normalization test is:
00068 *>
00069 *>         RESULT(2) =      max       | M(v(j)) - 1 | / ( n ulp )
00070 *>                    eigenvectors v(j)
00071 *> \endverbatim
00072 *
00073 *  Arguments:
00074 *  ==========
00075 *
00076 *> \param[in] LEFT
00077 *> \verbatim
00078 *>          LEFT is LOGICAL
00079 *>          =.TRUE.:  The eigenvectors in the columns of E are assumed
00080 *>                    to be *left* eigenvectors.
00081 *>          =.FALSE.: The eigenvectors in the columns of E are assumed
00082 *>                    to be *right* eigenvectors.
00083 *> \endverbatim
00084 *>
00085 *> \param[in] N
00086 *> \verbatim
00087 *>          N is INTEGER
00088 *>          The size of the matrices.  If it is zero, SGET52 does
00089 *>          nothing.  It must be at least zero.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] A
00093 *> \verbatim
00094 *>          A is REAL array, dimension (LDA, N)
00095 *>          The matrix A.
00096 *> \endverbatim
00097 *>
00098 *> \param[in] LDA
00099 *> \verbatim
00100 *>          LDA is INTEGER
00101 *>          The leading dimension of A.  It must be at least 1
00102 *>          and at least N.
00103 *> \endverbatim
00104 *>
00105 *> \param[in] B
00106 *> \verbatim
00107 *>          B is REAL array, dimension (LDB, N)
00108 *>          The matrix B.
00109 *> \endverbatim
00110 *>
00111 *> \param[in] LDB
00112 *> \verbatim
00113 *>          LDB is INTEGER
00114 *>          The leading dimension of B.  It must be at least 1
00115 *>          and at least N.
00116 *> \endverbatim
00117 *>
00118 *> \param[in] E
00119 *> \verbatim
00120 *>          E is REAL array, dimension (LDE, N)
00121 *>          The matrix of eigenvectors.  It must be O( 1 ).  Complex
00122 *>          eigenvalues and eigenvectors always come in pairs, the
00123 *>          eigenvalue and its conjugate being stored in adjacent
00124 *>          elements of ALPHAR, ALPHAI, and BETA.  Thus, if a(j)/b(j)
00125 *>          and a(j+1)/b(j+1) are a complex conjugate pair of
00126 *>          generalized eigenvalues, then E(,j) contains the real part
00127 *>          of the eigenvector and E(,j+1) contains the imaginary part.
00128 *>          Note that whether E(,j) is a real eigenvector or part of a
00129 *>          complex one is specified by whether ALPHAI(j) is zero or not.
00130 *> \endverbatim
00131 *>
00132 *> \param[in] LDE
00133 *> \verbatim
00134 *>          LDE is INTEGER
00135 *>          The leading dimension of E.  It must be at least 1 and at
00136 *>          least N.
00137 *> \endverbatim
00138 *>
00139 *> \param[in] ALPHAR
00140 *> \verbatim
00141 *>          ALPHAR is REAL array, dimension (N)
00142 *>          The real parts of the values a(j) as described above, which,
00143 *>          along with b(j), define the generalized eigenvalues.
00144 *>          Complex eigenvalues always come in complex conjugate pairs
00145 *>          a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent
00146 *>          elements in ALPHAR, ALPHAI, and BETA.  Thus, if the j-th
00147 *>          and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1)
00148 *>          is assumed to be equal to ALPHAR(j)/BETA(j).
00149 *> \endverbatim
00150 *>
00151 *> \param[in] ALPHAI
00152 *> \verbatim
00153 *>          ALPHAI is REAL array, dimension (N)
00154 *>          The imaginary parts of the values a(j) as described above,
00155 *>          which, along with b(j), define the generalized eigenvalues.
00156 *>          If ALPHAI(j)=0, then the eigenvalue is real, otherwise it
00157 *>          is part of a complex conjugate pair.  Complex eigenvalues
00158 *>          always come in complex conjugate pairs a(j)/b(j) and
00159 *>          a(j+1)/b(j+1), which are stored in adjacent elements in
00160 *>          ALPHAR, ALPHAI, and BETA.  Thus, if the j-th and (j+1)-st
00161 *>          eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to
00162 *>          be equal to  -ALPHAI(j)/BETA(j).  Also, nonzero values in
00163 *>          ALPHAI are assumed to always come in adjacent pairs.
00164 *> \endverbatim
00165 *>
00166 *> \param[in] BETA
00167 *> \verbatim
00168 *>          BETA is REAL array, dimension (N)
00169 *>          The values b(j) as described above, which, along with a(j),
00170 *>          define the generalized eigenvalues.
00171 *> \endverbatim
00172 *>
00173 *> \param[out] WORK
00174 *> \verbatim
00175 *>          WORK is REAL array, dimension (N**2+N)
00176 *> \endverbatim
00177 *>
00178 *> \param[out] RESULT
00179 *> \verbatim
00180 *>          RESULT is REAL array, dimension (2)
00181 *>          The values computed by the test described above.  If A E or
00182 *>          B E is likely to overflow, then RESULT(1:2) is set to
00183 *>          10 / ulp.
00184 *> \endverbatim
00185 *
00186 *  Authors:
00187 *  ========
00188 *
00189 *> \author Univ. of Tennessee 
00190 *> \author Univ. of California Berkeley 
00191 *> \author Univ. of Colorado Denver 
00192 *> \author NAG Ltd. 
00193 *
00194 *> \date November 2011
00195 *
00196 *> \ingroup single_eig
00197 *
00198 *  =====================================================================
00199       SUBROUTINE SGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
00200      $                   ALPHAI, BETA, WORK, RESULT )
00201 *
00202 *  -- LAPACK test routine (version 3.4.0) --
00203 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00204 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00205 *     November 2011
00206 *
00207 *     .. Scalar Arguments ..
00208       LOGICAL            LEFT
00209       INTEGER            LDA, LDB, LDE, N
00210 *     ..
00211 *     .. Array Arguments ..
00212       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00213      $                   B( LDB, * ), BETA( * ), E( LDE, * ),
00214      $                   RESULT( 2 ), WORK( * )
00215 *     ..
00216 *
00217 *  =====================================================================
00218 *
00219 *     .. Parameters ..
00220       REAL               ZERO, ONE, TEN
00221       PARAMETER          ( ZERO = 0.0, ONE = 1.0, TEN = 10.0 )
00222 *     ..
00223 *     .. Local Scalars ..
00224       LOGICAL            ILCPLX
00225       CHARACTER          NORMAB, TRANS
00226       INTEGER            J, JVEC
00227       REAL               ABMAX, ACOEF, ALFMAX, ANORM, BCOEFI, BCOEFR,
00228      $                   BETMAX, BNORM, ENORM, ENRMER, ERRNRM, SAFMAX,
00229      $                   SAFMIN, SALFI, SALFR, SBETA, SCALE, TEMP1, ULP
00230 *     ..
00231 *     .. External Functions ..
00232       REAL               SLAMCH, SLANGE
00233       EXTERNAL           SLAMCH, SLANGE
00234 *     ..
00235 *     .. External Subroutines ..
00236       EXTERNAL           SGEMV
00237 *     ..
00238 *     .. Intrinsic Functions ..
00239       INTRINSIC          ABS, MAX, REAL
00240 *     ..
00241 *     .. Executable Statements ..
00242 *
00243       RESULT( 1 ) = ZERO
00244       RESULT( 2 ) = ZERO
00245       IF( N.LE.0 )
00246      $   RETURN
00247 *
00248       SAFMIN = SLAMCH( 'Safe minimum' )
00249       SAFMAX = ONE / SAFMIN
00250       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00251 *
00252       IF( LEFT ) THEN
00253          TRANS = 'T'
00254          NORMAB = 'I'
00255       ELSE
00256          TRANS = 'N'
00257          NORMAB = 'O'
00258       END IF
00259 *
00260 *     Norm of A, B, and E:
00261 *
00262       ANORM = MAX( SLANGE( NORMAB, N, N, A, LDA, WORK ), SAFMIN )
00263       BNORM = MAX( SLANGE( NORMAB, N, N, B, LDB, WORK ), SAFMIN )
00264       ENORM = MAX( SLANGE( 'O', N, N, E, LDE, WORK ), ULP )
00265       ALFMAX = SAFMAX / MAX( ONE, BNORM )
00266       BETMAX = SAFMAX / MAX( ONE, ANORM )
00267 *
00268 *     Compute error matrix.
00269 *     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
00270 *
00271       ILCPLX = .FALSE.
00272       DO 10 JVEC = 1, N
00273          IF( ILCPLX ) THEN
00274 *
00275 *           2nd Eigenvalue/-vector of pair -- do nothing
00276 *
00277             ILCPLX = .FALSE.
00278          ELSE
00279             SALFR = ALPHAR( JVEC )
00280             SALFI = ALPHAI( JVEC )
00281             SBETA = BETA( JVEC )
00282             IF( SALFI.EQ.ZERO ) THEN
00283 *
00284 *              Real eigenvalue and -vector
00285 *
00286                ABMAX = MAX( ABS( SALFR ), ABS( SBETA ) )
00287                IF( ABS( SALFR ).GT.ALFMAX .OR. ABS( SBETA ).GT.
00288      $             BETMAX .OR. ABMAX.LT.ONE ) THEN
00289                   SCALE = ONE / MAX( ABMAX, SAFMIN )
00290                   SALFR = SCALE*SALFR
00291                   SBETA = SCALE*SBETA
00292                END IF
00293                SCALE = ONE / MAX( ABS( SALFR )*BNORM,
00294      $                 ABS( SBETA )*ANORM, SAFMIN )
00295                ACOEF = SCALE*SBETA
00296                BCOEFR = SCALE*SALFR
00297                CALL SGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
00298      $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
00299                CALL SGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
00300      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
00301             ELSE
00302 *
00303 *              Complex conjugate pair
00304 *
00305                ILCPLX = .TRUE.
00306                IF( JVEC.EQ.N ) THEN
00307                   RESULT( 1 ) = TEN / ULP
00308                   RETURN
00309                END IF
00310                ABMAX = MAX( ABS( SALFR )+ABS( SALFI ), ABS( SBETA ) )
00311                IF( ABS( SALFR )+ABS( SALFI ).GT.ALFMAX .OR.
00312      $             ABS( SBETA ).GT.BETMAX .OR. ABMAX.LT.ONE ) THEN
00313                   SCALE = ONE / MAX( ABMAX, SAFMIN )
00314                   SALFR = SCALE*SALFR
00315                   SALFI = SCALE*SALFI
00316                   SBETA = SCALE*SBETA
00317                END IF
00318                SCALE = ONE / MAX( ( ABS( SALFR )+ABS( SALFI ) )*BNORM,
00319      $                 ABS( SBETA )*ANORM, SAFMIN )
00320                ACOEF = SCALE*SBETA
00321                BCOEFR = SCALE*SALFR
00322                BCOEFI = SCALE*SALFI
00323                IF( LEFT ) THEN
00324                   BCOEFI = -BCOEFI
00325                END IF
00326 *
00327                CALL SGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
00328      $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
00329                CALL SGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
00330      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
00331                CALL SGEMV( TRANS, N, N, BCOEFI, B, LDA, E( 1, JVEC+1 ),
00332      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
00333 *
00334                CALL SGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC+1 ),
00335      $                     1, ZERO, WORK( N*JVEC+1 ), 1 )
00336                CALL SGEMV( TRANS, N, N, -BCOEFI, B, LDA, E( 1, JVEC ),
00337      $                     1, ONE, WORK( N*JVEC+1 ), 1 )
00338                CALL SGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC+1 ),
00339      $                     1, ONE, WORK( N*JVEC+1 ), 1 )
00340             END IF
00341          END IF
00342    10 CONTINUE
00343 *
00344       ERRNRM = SLANGE( 'One', N, N, WORK, N, WORK( N**2+1 ) ) / ENORM
00345 *
00346 *     Compute RESULT(1)
00347 *
00348       RESULT( 1 ) = ERRNRM / ULP
00349 *
00350 *     Normalization of E:
00351 *
00352       ENRMER = ZERO
00353       ILCPLX = .FALSE.
00354       DO 40 JVEC = 1, N
00355          IF( ILCPLX ) THEN
00356             ILCPLX = .FALSE.
00357          ELSE
00358             TEMP1 = ZERO
00359             IF( ALPHAI( JVEC ).EQ.ZERO ) THEN
00360                DO 20 J = 1, N
00361                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
00362    20          CONTINUE
00363                ENRMER = MAX( ENRMER, TEMP1-ONE )
00364             ELSE
00365                ILCPLX = .TRUE.
00366                DO 30 J = 1, N
00367                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
00368      $                    ABS( E( J, JVEC+1 ) ) )
00369    30          CONTINUE
00370                ENRMER = MAX( ENRMER, TEMP1-ONE )
00371             END IF
00372          END IF
00373    40 CONTINUE
00374 *
00375 *     Compute RESULT(2) : the normalization error in E.
00376 *
00377       RESULT( 2 ) = ENRMER / ( REAL( N )*ULP )
00378 *
00379       RETURN
00380 *
00381 *     End of SGET52
00382 *
00383       END
 All Files Functions