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