LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cget52.f
Go to the documentation of this file.
00001 *> \brief \b CGET52
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 CGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA,
00012 *                          WORK, RWORK, RESULT )
00013 * 
00014 *       .. Scalar Arguments ..
00015 *       LOGICAL            LEFT
00016 *       INTEGER            LDA, LDB, LDE, N
00017 *       ..
00018 *       .. Array Arguments ..
00019 *       REAL               RESULT( 2 ), RWORK( * )
00020 *       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
00021 *      $                   BETA( * ), E( LDE, * ), WORK( * )
00022 *       ..
00023 *  
00024 *
00025 *> \par Purpose:
00026 *  =============
00027 *>
00028 *> \verbatim
00029 *>
00030 *> CGET52  does an eigenvector check for the generalized eigenvalue
00031 *> problem.
00032 *>
00033 *> The basic test for right eigenvectors is:
00034 *>
00035 *>                           | b(i) A E(i) -  a(i) B E(i) |
00036 *>         RESULT(1) = max   -------------------------------
00037 *>                      i    n ulp max( |b(i) A|, |a(i) B| )
00038 *>
00039 *> using the 1-norm.  Here, a(i)/b(i) = w is the i-th generalized
00040 *> eigenvalue of A - w B, or, equivalently, b(i)/a(i) = m is the i-th
00041 *> generalized eigenvalue of m A - B.
00042 *>
00043 *>                         H   H  _      _
00044 *> For left eigenvectors, A , B , a, and b  are used.
00045 *>
00046 *> CGET52 also tests the normalization of E.  Each eigenvector is
00047 *> supposed to be normalized so that the maximum "absolute value"
00048 *> of its elements is 1, where in this case, "absolute value"
00049 *> of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
00050 *> maximum "absolute value" norm of a vector v  M(v).  
00051 *> if a(i)=b(i)=0, then the eigenvector is set to be the jth coordinate
00052 *> vector. The normalization test is:
00053 *>
00054 *>         RESULT(2) =      max       | M(v(i)) - 1 | / ( n ulp )
00055 *>                    eigenvectors v(i)
00056 *> \endverbatim
00057 *
00058 *  Arguments:
00059 *  ==========
00060 *
00061 *> \param[in] LEFT
00062 *> \verbatim
00063 *>          LEFT is LOGICAL
00064 *>          =.TRUE.:  The eigenvectors in the columns of E are assumed
00065 *>                    to be *left* eigenvectors.
00066 *>          =.FALSE.: The eigenvectors in the columns of E are assumed
00067 *>                    to be *right* eigenvectors.
00068 *> \endverbatim
00069 *>
00070 *> \param[in] N
00071 *> \verbatim
00072 *>          N is INTEGER
00073 *>          The size of the matrices.  If it is zero, CGET52 does
00074 *>          nothing.  It must be at least zero.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] A
00078 *> \verbatim
00079 *>          A is COMPLEX array, dimension (LDA, N)
00080 *>          The matrix A.
00081 *> \endverbatim
00082 *>
00083 *> \param[in] LDA
00084 *> \verbatim
00085 *>          LDA is INTEGER
00086 *>          The leading dimension of A.  It must be at least 1
00087 *>          and at least N.
00088 *> \endverbatim
00089 *>
00090 *> \param[in] B
00091 *> \verbatim
00092 *>          B is COMPLEX array, dimension (LDB, N)
00093 *>          The matrix B.
00094 *> \endverbatim
00095 *>
00096 *> \param[in] LDB
00097 *> \verbatim
00098 *>          LDB is INTEGER
00099 *>          The leading dimension of B.  It must be at least 1
00100 *>          and at least N.
00101 *> \endverbatim
00102 *>
00103 *> \param[in] E
00104 *> \verbatim
00105 *>          E is COMPLEX array, dimension (LDE, N)
00106 *>          The matrix of eigenvectors.  It must be O( 1 ).
00107 *> \endverbatim
00108 *>
00109 *> \param[in] LDE
00110 *> \verbatim
00111 *>          LDE is INTEGER
00112 *>          The leading dimension of E.  It must be at least 1 and at
00113 *>          least N.
00114 *> \endverbatim
00115 *>
00116 *> \param[in] ALPHA
00117 *> \verbatim
00118 *>          ALPHA is COMPLEX array, dimension (N)
00119 *>          The values a(i) as described above, which, along with b(i),
00120 *>          define the generalized eigenvalues.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] BETA
00124 *> \verbatim
00125 *>          BETA is COMPLEX array, dimension (N)
00126 *>          The values b(i) as described above, which, along with a(i),
00127 *>          define the generalized eigenvalues.
00128 *> \endverbatim
00129 *>
00130 *> \param[out] WORK
00131 *> \verbatim
00132 *>          WORK is COMPLEX array, dimension (N**2)
00133 *> \endverbatim
00134 *>
00135 *> \param[out] RWORK
00136 *> \verbatim
00137 *>          RWORK is REAL array, dimension (N)
00138 *> \endverbatim
00139 *>
00140 *> \param[out] RESULT
00141 *> \verbatim
00142 *>          RESULT is REAL array, dimension (2)
00143 *>          The values computed by the test described above.  If A E or
00144 *>          B E is likely to overflow, then RESULT(1:2) is set to
00145 *>          10 / ulp.
00146 *> \endverbatim
00147 *
00148 *  Authors:
00149 *  ========
00150 *
00151 *> \author Univ. of Tennessee 
00152 *> \author Univ. of California Berkeley 
00153 *> \author Univ. of Colorado Denver 
00154 *> \author NAG Ltd. 
00155 *
00156 *> \date November 2011
00157 *
00158 *> \ingroup complex_eig
00159 *
00160 *  =====================================================================
00161       SUBROUTINE CGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHA, BETA,
00162      $                   WORK, RWORK, RESULT )
00163 *
00164 *  -- LAPACK test routine (version 3.4.0) --
00165 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00166 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00167 *     November 2011
00168 *
00169 *     .. Scalar Arguments ..
00170       LOGICAL            LEFT
00171       INTEGER            LDA, LDB, LDE, N
00172 *     ..
00173 *     .. Array Arguments ..
00174       REAL               RESULT( 2 ), RWORK( * )
00175       COMPLEX            A( LDA, * ), ALPHA( * ), B( LDB, * ),
00176      $                   BETA( * ), E( LDE, * ), WORK( * )
00177 *     ..
00178 *
00179 *  =====================================================================
00180 *
00181 *     .. Parameters ..
00182       REAL               ZERO, ONE
00183       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
00184       COMPLEX            CZERO, CONE
00185       PARAMETER          ( CZERO = ( 0.0E+0, 0.0E+0 ),
00186      $                   CONE = ( 1.0E+0, 0.0E+0 ) )
00187 *     ..
00188 *     .. Local Scalars ..
00189       CHARACTER          NORMAB, TRANS
00190       INTEGER            J, JVEC
00191       REAL               ABMAX, ALFMAX, ANORM, BETMAX, BNORM, ENORM,
00192      $                   ENRMER, ERRNRM, SAFMAX, SAFMIN, SCALE, TEMP1,
00193      $                   ULP
00194       COMPLEX            ACOEFF, ALPHAI, BCOEFF, BETAI, X
00195 *     ..
00196 *     .. External Functions ..
00197       REAL               CLANGE, SLAMCH
00198       EXTERNAL           CLANGE, SLAMCH
00199 *     ..
00200 *     .. External Subroutines ..
00201       EXTERNAL           CGEMV
00202 *     ..
00203 *     .. Intrinsic Functions ..
00204       INTRINSIC          ABS, AIMAG, CONJG, MAX, REAL
00205 *     ..
00206 *     .. Statement Functions ..
00207       REAL               ABS1
00208 *     ..
00209 *     .. Statement Function definitions ..
00210       ABS1( X ) = ABS( REAL( X ) ) + ABS( AIMAG( X ) )
00211 *     ..
00212 *     .. Executable Statements ..
00213 *
00214       RESULT( 1 ) = ZERO
00215       RESULT( 2 ) = ZERO
00216       IF( N.LE.0 )
00217      $   RETURN
00218 *
00219       SAFMIN = SLAMCH( 'Safe minimum' )
00220       SAFMAX = ONE / SAFMIN
00221       ULP = SLAMCH( 'Epsilon' )*SLAMCH( 'Base' )
00222 *
00223       IF( LEFT ) THEN
00224          TRANS = 'C'
00225          NORMAB = 'I'
00226       ELSE
00227          TRANS = 'N'
00228          NORMAB = 'O'
00229       END IF
00230 *
00231 *     Norm of A, B, and E:
00232 *
00233       ANORM = MAX( CLANGE( NORMAB, N, N, A, LDA, RWORK ), SAFMIN )
00234       BNORM = MAX( CLANGE( NORMAB, N, N, B, LDB, RWORK ), SAFMIN )
00235       ENORM = MAX( CLANGE( 'O', N, N, E, LDE, RWORK ), ULP )
00236       ALFMAX = SAFMAX / MAX( ONE, BNORM )
00237       BETMAX = SAFMAX / MAX( ONE, ANORM )
00238 *
00239 *     Compute error matrix.
00240 *     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
00241 *
00242       DO 10 JVEC = 1, N
00243          ALPHAI = ALPHA( JVEC )
00244          BETAI = BETA( JVEC )
00245          ABMAX = MAX( ABS1( ALPHAI ), ABS1( BETAI ) )
00246          IF( ABS1( ALPHAI ).GT.ALFMAX .OR. ABS1( BETAI ).GT.BETMAX .OR.
00247      $       ABMAX.LT.ONE ) THEN
00248             SCALE = ONE / MAX( ABMAX, SAFMIN )
00249             ALPHAI = SCALE*ALPHAI
00250             BETAI = SCALE*BETAI
00251          END IF
00252          SCALE = ONE / MAX( ABS1( ALPHAI )*BNORM, ABS1( BETAI )*ANORM,
00253      $           SAFMIN )
00254          ACOEFF = SCALE*BETAI
00255          BCOEFF = SCALE*ALPHAI
00256          IF( LEFT ) THEN
00257             ACOEFF = CONJG( ACOEFF )
00258             BCOEFF = CONJG( BCOEFF )
00259          END IF
00260          CALL CGEMV( TRANS, N, N, ACOEFF, A, LDA, E( 1, JVEC ), 1,
00261      $               CZERO, WORK( N*( JVEC-1 )+1 ), 1 )
00262          CALL CGEMV( TRANS, N, N, -BCOEFF, B, LDA, E( 1, JVEC ), 1,
00263      $               CONE, WORK( N*( JVEC-1 )+1 ), 1 )
00264    10 CONTINUE
00265 *
00266       ERRNRM = CLANGE( 'One', N, N, WORK, N, RWORK ) / ENORM
00267 *
00268 *     Compute RESULT(1)
00269 *
00270       RESULT( 1 ) = ERRNRM / ULP
00271 *
00272 *     Normalization of E:
00273 *
00274       ENRMER = ZERO
00275       DO 30 JVEC = 1, N
00276          TEMP1 = ZERO
00277          DO 20 J = 1, N
00278             TEMP1 = MAX( TEMP1, ABS1( E( J, JVEC ) ) )
00279    20    CONTINUE
00280          ENRMER = MAX( ENRMER, TEMP1-ONE )
00281    30 CONTINUE
00282 *
00283 *     Compute RESULT(2) : the normalization error in E.
00284 *
00285       RESULT( 2 ) = ENRMER / ( REAL( N )*ULP )
00286 *
00287       RETURN
00288 *
00289 *     End of CGET52
00290 *
00291       END
 All Files Functions