LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine dget52 ( logical  LEFT,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( ldb, * )  B,
integer  LDB,
double precision, dimension( lde, * )  E,
integer  LDE,
double precision, dimension( * )  ALPHAR,
double precision, dimension( * )  ALPHAI,
double precision, dimension( * )  BETA,
double precision, dimension( * )  WORK,
double precision, dimension( 2 )  RESULT 
)

DGET52

Purpose: 
 DGET52  does an eigenvector check for the generalized eigenvalue
 problem.

 The basic test for right eigenvectors is:

                           | b(j) A E(j) -  a(j) B E(j) |
         RESULT(1) = max   -------------------------------
                      j    n ulp max( |b(j) A|, |a(j) B| )

 using the 1-norm.  Here, a(j)/b(j) = w is the j-th generalized
 eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th
 generalized eigenvalue of m A - B.

 For real eigenvalues, the test is straightforward.  For complex
 eigenvalues, E(j) and a(j) are complex, represented by
 Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that
 eigenvector becomes

                 max( |Wr|, |Wi| )
     --------------------------------------------
     n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| )

 where

     Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j)

     Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j)

                         T   T  _
 For left eigenvectors, A , B , a, and b  are used.

 DGET52 also tests the normalization of E.  Each eigenvector is
 supposed to be normalized so that the maximum "absolute value"
 of its elements is 1, where in this case, "absolute value"
 of a complex value x is  |Re(x)| + |Im(x)| ; let us call this
 maximum "absolute value" norm of a vector v  M(v).
 if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate
 vector.  The normalization test is:

         RESULT(2) =      max       | M(v(j)) - 1 | / ( n ulp )
                    eigenvectors v(j)
Parameters
[in]LEFT
          LEFT is LOGICAL
          =.TRUE.:  The eigenvectors in the columns of E are assumed
                    to be *left* eigenvectors.
          =.FALSE.: The eigenvectors in the columns of E are assumed
                    to be *right* eigenvectors.
[in]N
          N is INTEGER
          The size of the matrices.  If it is zero, DGET52 does
          nothing.  It must be at least zero.
[in]A
          A is DOUBLE PRECISION array, dimension (LDA, N)
          The matrix A.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  It must be at least 1
          and at least N.
[in]B
          B is DOUBLE PRECISION array, dimension (LDB, N)
          The matrix B.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  It must be at least 1
          and at least N.
[in]E
          E is DOUBLE PRECISION array, dimension (LDE, N)
          The matrix of eigenvectors.  It must be O( 1 ).  Complex
          eigenvalues and eigenvectors always come in pairs, the
          eigenvalue and its conjugate being stored in adjacent
          elements of ALPHAR, ALPHAI, and BETA.  Thus, if a(j)/b(j)
          and a(j+1)/b(j+1) are a complex conjugate pair of
          generalized eigenvalues, then E(,j) contains the real part
          of the eigenvector and E(,j+1) contains the imaginary part.
          Note that whether E(,j) is a real eigenvector or part of a
          complex one is specified by whether ALPHAI(j) is zero or not.
[in]LDE
          LDE is INTEGER
          The leading dimension of E.  It must be at least 1 and at
          least N.
[in]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (N)
          The real parts of the values a(j) as described above, which,
          along with b(j), define the generalized eigenvalues.
          Complex eigenvalues always come in complex conjugate pairs
          a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent
          elements in ALPHAR, ALPHAI, and BETA.  Thus, if the j-th
          and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1)
          is assumed to be equal to ALPHAR(j)/BETA(j).
[in]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (N)
          The imaginary parts of the values a(j) as described above,
          which, along with b(j), define the generalized eigenvalues.
          If ALPHAI(j)=0, then the eigenvalue is real, otherwise it
          is part of a complex conjugate pair.  Complex eigenvalues
          always come in complex conjugate pairs a(j)/b(j) and
          a(j+1)/b(j+1), which are stored in adjacent elements in
          ALPHAR, ALPHAI, and BETA.  Thus, if the j-th and (j+1)-st
          eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to
          be equal to  -ALPHAI(j)/BETA(j).  Also, nonzero values in
          ALPHAI are assumed to always come in adjacent pairs.
[in]BETA
          BETA is DOUBLE PRECISION array, dimension (N)
          The values b(j) as described above, which, along with a(j),
          define the generalized eigenvalues.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (N**2+N)
[out]RESULT
          RESULT is DOUBLE PRECISION array, dimension (2)
          The values computed by the test described above.  If A E or
          B E is likely to overflow, then RESULT(1:2) is set to
          10 / ulp.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 201 of file dget52.f.

201 *
202 * -- LAPACK test routine (version 3.4.0) --
203 * -- LAPACK is a software package provided by Univ. of Tennessee, --
204 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
205 * November 2011
206 *
207 * .. Scalar Arguments ..
208  LOGICAL left
209  INTEGER lda, ldb, lde, n
210 * ..
211 * .. Array Arguments ..
212  DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
213  $ b( ldb, * ), beta( * ), e( lde, * ),
214  $ result( 2 ), work( * )
215 * ..
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220  DOUBLE PRECISION zero, one, ten
221  parameter ( zero = 0.0d0, one = 1.0d0, ten = 10.0d0 )
222 * ..
223 * .. Local Scalars ..
224  LOGICAL ilcplx
225  CHARACTER normab, trans
226  INTEGER j, jvec
227  DOUBLE PRECISION abmax, acoef, alfmax, anorm, bcoefi, bcoefr,
228  $ betmax, bnorm, enorm, enrmer, errnrm, safmax,
229  $ safmin, salfi, salfr, sbeta, scale, temp1, ulp
230 * ..
231 * .. External Functions ..
232  DOUBLE PRECISION dlamch, dlange
233  EXTERNAL dlamch, dlange
234 * ..
235 * .. External Subroutines ..
236  EXTERNAL dgemv
237 * ..
238 * .. Intrinsic Functions ..
239  INTRINSIC abs, dble, max
240 * ..
241 * .. Executable Statements ..
242 *
243  result( 1 ) = zero
244  result( 2 ) = zero
245  IF( n.LE.0 )
246  $ RETURN
247 *
248  safmin = dlamch( 'Safe minimum' )
249  safmax = one / safmin
250  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
251 *
252  IF( left ) THEN
253  trans = 'T'
254  normab = 'I'
255  ELSE
256  trans = 'N'
257  normab = 'O'
258  END IF
259 *
260 * Norm of A, B, and E:
261 *
262  anorm = max( dlange( normab, n, n, a, lda, work ), safmin )
263  bnorm = max( dlange( normab, n, n, b, ldb, work ), safmin )
264  enorm = max( dlange( 'O', n, n, e, lde, work ), ulp )
265  alfmax = safmax / max( one, bnorm )
266  betmax = safmax / max( one, anorm )
267 *
268 * Compute error matrix.
269 * Column i = ( b(i) A - a(i) B ) E(i) / max( |a(i) B| |b(i) A| )
270 *
271  ilcplx = .false.
272  DO 10 jvec = 1, n
273  IF( ilcplx ) THEN
274 *
275 * 2nd Eigenvalue/-vector of pair -- do nothing
276 *
277  ilcplx = .false.
278  ELSE
279  salfr = alphar( jvec )
280  salfi = alphai( jvec )
281  sbeta = beta( jvec )
282  IF( salfi.EQ.zero ) THEN
283 *
284 * Real eigenvalue and -vector
285 *
286  abmax = max( abs( salfr ), abs( sbeta ) )
287  IF( abs( salfr ).GT.alfmax .OR. abs( sbeta ).GT.
288  $ betmax .OR. abmax.LT.one ) THEN
289  scale = one / max( abmax, safmin )
290  salfr = scale*salfr
291  sbeta = scale*sbeta
292  END IF
293  scale = one / max( abs( salfr )*bnorm,
294  $ abs( sbeta )*anorm, safmin )
295  acoef = scale*sbeta
296  bcoefr = scale*salfr
297  CALL dgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
298  $ zero, work( n*( jvec-1 )+1 ), 1 )
299  CALL dgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec ),
300  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
301  ELSE
302 *
303 * Complex conjugate pair
304 *
305  ilcplx = .true.
306  IF( jvec.EQ.n ) THEN
307  result( 1 ) = ten / ulp
308  RETURN
309  END IF
310  abmax = max( abs( salfr )+abs( salfi ), abs( sbeta ) )
311  IF( abs( salfr )+abs( salfi ).GT.alfmax .OR.
312  $ abs( sbeta ).GT.betmax .OR. abmax.LT.one ) THEN
313  scale = one / max( abmax, safmin )
314  salfr = scale*salfr
315  salfi = scale*salfi
316  sbeta = scale*sbeta
317  END IF
318  scale = one / max( ( abs( salfr )+abs( salfi ) )*bnorm,
319  $ abs( sbeta )*anorm, safmin )
320  acoef = scale*sbeta
321  bcoefr = scale*salfr
322  bcoefi = scale*salfi
323  IF( left ) THEN
324  bcoefi = -bcoefi
325  END IF
326 *
327  CALL dgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
328  $ zero, work( n*( jvec-1 )+1 ), 1 )
329  CALL dgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec ),
330  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
331  CALL dgemv( trans, n, n, bcoefi, b, lda, e( 1, jvec+1 ),
332  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
333 *
334  CALL dgemv( trans, n, n, acoef, a, lda, e( 1, jvec+1 ),
335  $ 1, zero, work( n*jvec+1 ), 1 )
336  CALL dgemv( trans, n, n, -bcoefi, b, lda, e( 1, jvec ),
337  $ 1, one, work( n*jvec+1 ), 1 )
338  CALL dgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec+1 ),
339  $ 1, one, work( n*jvec+1 ), 1 )
340  END IF
341  END IF
342  10 CONTINUE
343 *
344  errnrm = dlange( 'One', n, n, work, n, work( n**2+1 ) ) / enorm
345 *
346 * Compute RESULT(1)
347 *
348  result( 1 ) = errnrm / ulp
349 *
350 * Normalization of E:
351 *
352  enrmer = zero
353  ilcplx = .false.
354  DO 40 jvec = 1, n
355  IF( ilcplx ) THEN
356  ilcplx = .false.
357  ELSE
358  temp1 = zero
359  IF( alphai( jvec ).EQ.zero ) THEN
360  DO 20 j = 1, n
361  temp1 = max( temp1, abs( e( j, jvec ) ) )
362  20 CONTINUE
363  enrmer = max( enrmer, temp1-one )
364  ELSE
365  ilcplx = .true.
366  DO 30 j = 1, n
367  temp1 = max( temp1, abs( e( j, jvec ) )+
368  $ abs( e( j, jvec+1 ) ) )
369  30 CONTINUE
370  enrmer = max( enrmer, temp1-one )
371  END IF
372  END IF
373  40 CONTINUE
374 *
375 * Compute RESULT(2) : the normalization error in E.
376 *
377  result( 2 ) = enrmer / ( dble( n )*ulp )
378 *
379  RETURN
380 *
381 * End of DGET52
382 *
dlamch
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
lde
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
dlange
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:116

Here is the call graph for this function:

Here is the caller graph for this function: