 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ 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.```

Definition at line 197 of file dget52.f.

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