LAPACK: TESTING/EIG/dget52.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
00002      $                   ALPHAI, BETA, WORK, RESULT )
00003 *
00004 *  -- LAPACK test routine (version 3.1) --
00005 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
00006 *     November 2006
00007 *
00008 *     .. Scalar Arguments ..
00009       LOGICAL            LEFT
00010       INTEGER            LDA, LDB, LDE, N
00011 *     ..
00012 *     .. Array Arguments ..
00013       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00014      $                   B( LDB, * ), BETA( * ), E( LDE, * ),
00015      $                   RESULT( 2 ), WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  DGET52  does an eigenvector check for the generalized eigenvalue
00022 *  problem.
00023 *
00024 *  The basic test for right eigenvectors is:
00025 *
00026 *                            | b(j) A E(j) -  a(j) B E(j) |
00027 *          RESULT(1) = max   -------------------------------
00028 *                       j    n ulp max( |b(j) A|, |a(j) B| )
00029 *
00030 *  using the 1-norm.  Here, a(j)/b(j) = w is the j-th generalized
00031 *  eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th
00032 *  generalized eigenvalue of m A - B.
00033 *
00034 *  For real eigenvalues, the test is straightforward.  For complex
00035 *  eigenvalues, E(j) and a(j) are complex, represented by
00036 *  Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that
00037 *  eigenvector becomes
00038 *
00039 *                  max( |Wr|, |Wi| )
00040 *      --------------------------------------------
00041 *      n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| )
00042 *
00043 *  where
00044 *
00045 *      Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j)
00046 *
00047 *      Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j)
00048 *
00049 *                          T   T  _
00050 *  For left eigenvectors, A , B , a, and b  are used.
00051 *
00052 *  DGET52 also tests the normalization of E.  Each eigenvector is
00053 *  supposed to be normalized so that the maximum "absolute value"
00054 *  of its elements is 1, where in this case, "absolute value"
00055 *  of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
00056 *  maximum "absolute value" norm of a vector v  M(v).
00057 *  if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate
00058 *  vector.  The normalization test is:
00059 *
00060 *          RESULT(2) =      max       | M(v(j)) - 1 | / ( n ulp )
00061 *                     eigenvectors v(j)
00062 *
00063 *  Arguments
00064 *  =========
00065 *
00066 *  LEFT    (input) LOGICAL
00067 *          =.TRUE.:  The eigenvectors in the columns of E are assumed
00068 *                    to be *left* eigenvectors.
00069 *          =.FALSE.: The eigenvectors in the columns of E are assumed
00070 *                    to be *right* eigenvectors.
00071 *
00072 *  N       (input) INTEGER
00073 *          The size of the matrices.  If it is zero, DGET52 does
00074 *          nothing.  It must be at least zero.
00075 *
00076 *  A       (input) DOUBLE PRECISION array, dimension (LDA, N)
00077 *          The matrix A.
00078 *
00079 *  LDA     (input) INTEGER
00080 *          The leading dimension of A.  It must be at least 1
00081 *          and at least N.
00082 *
00083 *  B       (input) DOUBLE PRECISION array, dimension (LDB, N)
00084 *          The matrix B.
00085 *
00086 *  LDB     (input) INTEGER
00087 *          The leading dimension of B.  It must be at least 1
00088 *          and at least N.
00089 *
00090 *  E       (input) DOUBLE PRECISION array, dimension (LDE, N)
00091 *          The matrix of eigenvectors.  It must be O( 1 ).  Complex
00092 *          eigenvalues and eigenvectors always come in pairs, the
00093 *          eigenvalue and its conjugate being stored in adjacent
00094 *          elements of ALPHAR, ALPHAI, and BETA.  Thus, if a(j)/b(j)
00095 *          and a(j+1)/b(j+1) are a complex conjugate pair of
00096 *          generalized eigenvalues, then E(,j) contains the real part
00097 *          of the eigenvector and E(,j+1) contains the imaginary part.
00098 *          Note that whether E(,j) is a real eigenvector or part of a
00099 *          complex one is specified by whether ALPHAI(j) is zero or not.
00100 *
00101 *  LDE     (input) INTEGER
00102 *          The leading dimension of E.  It must be at least 1 and at
00103 *          least N.
00104 *
00105 *  ALPHAR  (input) DOUBLE PRECISION array, dimension (N)
00106 *          The real parts of the values a(j) as described above, which,
00107 *          along with b(j), define the generalized eigenvalues.
00108 *          Complex eigenvalues always come in complex conjugate pairs
00109 *          a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent
00110 *          elements in ALPHAR, ALPHAI, and BETA.  Thus, if the j-th
00111 *          and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1)
00112 *          is assumed to be equal to ALPHAR(j)/BETA(j).
00113 *
00114 *  ALPHAI  (input) DOUBLE PRECISION array, dimension (N)
00115 *          The imaginary parts of the values a(j) as described above,
00116 *          which, along with b(j), define the generalized eigenvalues.
00117 *          If ALPHAI(j)=0, then the eigenvalue is real, otherwise it
00118 *          is part of a complex conjugate pair.  Complex eigenvalues
00119 *          always come in complex conjugate pairs a(j)/b(j) and
00120 *          a(j+1)/b(j+1), which are stored in adjacent elements in
00121 *          ALPHAR, ALPHAI, and BETA.  Thus, if the j-th and (j+1)-st
00122 *          eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to
00123 *          be equal to  -ALPHAI(j)/BETA(j).  Also, nonzero values in
00124 *          ALPHAI are assumed to always come in adjacent pairs.
00125 *
00126 *  BETA    (input) DOUBLE PRECISION array, dimension (N)
00127 *          The values b(j) as described above, which, along with a(j),
00128 *          define the generalized eigenvalues.
00129 *
00130 *  WORK    (workspace) DOUBLE PRECISION array, dimension (N**2+N)
00131 *
00132 *  RESULT  (output) DOUBLE PRECISION array, dimension (2)
00133 *          The values computed by the test described above.  If A E or
00134 *          B E is likely to overflow, then RESULT(1:2) is set to
00135 *          10 / ulp.
00136 *
00137 *  =====================================================================
00138 *
00139 *     .. Parameters ..
00140       DOUBLE PRECISION   ZERO, ONE, TEN
00141       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TEN = 10.0D0 )
00142 *     ..
00143 *     .. Local Scalars ..
00144       LOGICAL            ILCPLX
00145       CHARACTER          NORMAB, TRANS
00146       INTEGER            J, JVEC
00147       DOUBLE PRECISION   ABMAX, ACOEF, ALFMAX, ANORM, BCOEFI, BCOEFR,
00148      $                   BETMAX, BNORM, ENORM, ENRMER, ERRNRM, SAFMAX,
00149      $                   SAFMIN, SALFI, SALFR, SBETA, SCALE, TEMP1, ULP
00150 *     ..
00151 *     .. External Functions ..
00152       DOUBLE PRECISION   DLAMCH, DLANGE
00153       EXTERNAL           DLAMCH, DLANGE
00154 *     ..
00155 *     .. External Subroutines ..
00156       EXTERNAL           DGEMV
00157 *     ..
00158 *     .. Intrinsic Functions ..
00159       INTRINSIC          ABS, DBLE, MAX
00160 *     ..
00161 *     .. Executable Statements ..
00162 *
00163       RESULT( 1 ) = ZERO
00164       RESULT( 2 ) = ZERO
00165       IF( N.LE.0 )
00166      $   RETURN
00167 *
00168       SAFMIN = DLAMCH( 'Safe minimum' )
00169       SAFMAX = ONE / SAFMIN
00170       ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
00171 *
00172       IF( LEFT ) THEN
00173          TRANS = 'T'
00174          NORMAB = 'I'
00175       ELSE
00176          TRANS = 'N'
00177          NORMAB = 'O'
00178       END IF
00179 *
00180 *     Norm of A, B, and E:
00181 *
00182       ANORM = MAX( DLANGE( NORMAB, N, N, A, LDA, WORK ), SAFMIN )
00183       BNORM = MAX( DLANGE( NORMAB, N, N, B, LDB, WORK ), SAFMIN )
00184       ENORM = MAX( DLANGE( 'O', N, N, E, LDE, WORK ), ULP )
00185       ALFMAX = SAFMAX / MAX( ONE, BNORM )
00186       BETMAX = SAFMAX / MAX( ONE, ANORM )
00187 *
00188 *     Compute error matrix.
00189 *     Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
00190 *
00191       ILCPLX = .FALSE.
00192       DO 10 JVEC = 1, N
00193          IF( ILCPLX ) THEN
00194 *
00195 *           2nd Eigenvalue/-vector of pair -- do nothing
00196 *
00197             ILCPLX = .FALSE.
00198          ELSE
00199             SALFR = ALPHAR( JVEC )
00200             SALFI = ALPHAI( JVEC )
00201             SBETA = BETA( JVEC )
00202             IF( SALFI.EQ.ZERO ) THEN
00203 *
00204 *              Real eigenvalue and -vector
00205 *
00206                ABMAX = MAX( ABS( SALFR ), ABS( SBETA ) )
00207                IF( ABS( SALFR ).GT.ALFMAX .OR. ABS( SBETA ).GT.
00208      $             BETMAX .OR. ABMAX.LT.ONE ) THEN
00209                   SCALE = ONE / MAX( ABMAX, SAFMIN )
00210                   SALFR = SCALE*SALFR
00211                   SBETA = SCALE*SBETA
00212                END IF
00213                SCALE = ONE / MAX( ABS( SALFR )*BNORM,
00214      $                 ABS( SBETA )*ANORM, SAFMIN )
00215                ACOEF = SCALE*SBETA
00216                BCOEFR = SCALE*SALFR
00217                CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
00218      $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
00219                CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
00220      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
00221             ELSE
00222 *
00223 *              Complex conjugate pair
00224 *
00225                ILCPLX = .TRUE.
00226                IF( JVEC.EQ.N ) THEN
00227                   RESULT( 1 ) = TEN / ULP
00228                   RETURN
00229                END IF
00230                ABMAX = MAX( ABS( SALFR )+ABS( SALFI ), ABS( SBETA ) )
00231                IF( ABS( SALFR )+ABS( SALFI ).GT.ALFMAX .OR.
00232      $             ABS( SBETA ).GT.BETMAX .OR. ABMAX.LT.ONE ) THEN
00233                   SCALE = ONE / MAX( ABMAX, SAFMIN )
00234                   SALFR = SCALE*SALFR
00235                   SALFI = SCALE*SALFI
00236                   SBETA = SCALE*SBETA
00237                END IF
00238                SCALE = ONE / MAX( ( ABS( SALFR )+ABS( SALFI ) )*BNORM,
00239      $                 ABS( SBETA )*ANORM, SAFMIN )
00240                ACOEF = SCALE*SBETA
00241                BCOEFR = SCALE*SALFR
00242                BCOEFI = SCALE*SALFI
00243                IF( LEFT ) THEN
00244                   BCOEFI = -BCOEFI
00245                END IF
00246 *
00247                CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC ), 1,
00248      $                     ZERO, WORK( N*( JVEC-1 )+1 ), 1 )
00249                CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC ),
00250      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
00251                CALL DGEMV( TRANS, N, N, BCOEFI, B, LDA, E( 1, JVEC+1 ),
00252      $                     1, ONE, WORK( N*( JVEC-1 )+1 ), 1 )
00253 *
00254                CALL DGEMV( TRANS, N, N, ACOEF, A, LDA, E( 1, JVEC+1 ),
00255      $                     1, ZERO, WORK( N*JVEC+1 ), 1 )
00256                CALL DGEMV( TRANS, N, N, -BCOEFI, B, LDA, E( 1, JVEC ),
00257      $                     1, ONE, WORK( N*JVEC+1 ), 1 )
00258                CALL DGEMV( TRANS, N, N, -BCOEFR, B, LDA, E( 1, JVEC+1 ),
00259      $                     1, ONE, WORK( N*JVEC+1 ), 1 )
00260             END IF
00261          END IF
00262    10 CONTINUE
00263 *
00264       ERRNRM = DLANGE( 'One', N, N, WORK, N, WORK( N**2+1 ) ) / ENORM
00265 *
00266 *     Compute RESULT(1)
00267 *
00268       RESULT( 1 ) = ERRNRM / ULP
00269 *
00270 *     Normalization of E:
00271 *
00272       ENRMER = ZERO
00273       ILCPLX = .FALSE.
00274       DO 40 JVEC = 1, N
00275          IF( ILCPLX ) THEN
00276             ILCPLX = .FALSE.
00277          ELSE
00278             TEMP1 = ZERO
00279             IF( ALPHAI( JVEC ).EQ.ZERO ) THEN
00280                DO 20 J = 1, N
00281                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) ) )
00282    20          CONTINUE
00283                ENRMER = MAX( ENRMER, TEMP1-ONE )
00284             ELSE
00285                ILCPLX = .TRUE.
00286                DO 30 J = 1, N
00287                   TEMP1 = MAX( TEMP1, ABS( E( J, JVEC ) )+
00288      $                    ABS( E( J, JVEC+1 ) ) )
00289    30          CONTINUE
00290                ENRMER = MAX( ENRMER, TEMP1-ONE )
00291             END IF
00292          END IF
00293    40 CONTINUE
00294 *
00295 *     Compute RESULT(2) : the normalization error in E.
00296 *
00297       RESULT( 2 ) = ENRMER / ( DBLE( N )*ULP )
00298 *
00299       RETURN
00300 *
00301 *     End of DGET52
00302 *
00303       END