LAPACK  3.8.0 LAPACK: Linear Algebra PACKage

## ◆ dget52()

 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.```
Date
December 2016

Definition at line 201 of file dget52.f.

201 *
202 * -- LAPACK test routine (version 3.7.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 * December 2016
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 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:65
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:158
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
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: