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

SGET52

Purpose:
``` SGET52  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.

SGET52 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, SGET52 does nothing. It must be at least zero.``` [in] A ``` A is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (N) The values b(j) as described above, which, along with a(j), define the generalized eigenvalues.``` [out] WORK ` WORK is REAL array, dimension (N**2+N)` [out] RESULT ``` RESULT is REAL 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
November 2011

Definition at line 201 of file sget52.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  REAL a( lda, * ), alphai( * ), alphar( * ),
213  \$ b( ldb, * ), beta( * ), e( lde, * ),
214  \$ result( 2 ), work( * )
215 * ..
216 *
217 * =====================================================================
218 *
219 * .. Parameters ..
220  REAL zero, one, ten
221  parameter ( zero = 0.0, one = 1.0, ten = 10.0 )
222 * ..
223 * .. Local Scalars ..
224  LOGICAL ilcplx
225  CHARACTER normab, trans
226  INTEGER j, jvec
227  REAL 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  REAL slamch, slange
233  EXTERNAL slamch, slange
234 * ..
235 * .. External Subroutines ..
236  EXTERNAL sgemv
237 * ..
238 * .. Intrinsic Functions ..
239  INTRINSIC abs, max, real
240 * ..
241 * .. Executable Statements ..
242 *
243  result( 1 ) = zero
244  result( 2 ) = zero
245  IF( n.LE.0 )
246  \$ RETURN
247 *
248  safmin = slamch( 'Safe minimum' )
249  safmax = one / safmin
250  ulp = slamch( 'Epsilon' )*slamch( '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( slange( normab, n, n, a, lda, work ), safmin )
263  bnorm = max( slange( normab, n, n, b, ldb, work ), safmin )
264  enorm = max( slange( '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 sgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
298  \$ zero, work( n*( jvec-1 )+1 ), 1 )
299  CALL sgemv( 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 sgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
328  \$ zero, work( n*( jvec-1 )+1 ), 1 )
329  CALL sgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec ),
330  \$ 1, one, work( n*( jvec-1 )+1 ), 1 )
331  CALL sgemv( trans, n, n, bcoefi, b, lda, e( 1, jvec+1 ),
332  \$ 1, one, work( n*( jvec-1 )+1 ), 1 )
333 *
334  CALL sgemv( trans, n, n, acoef, a, lda, e( 1, jvec+1 ),
335  \$ 1, zero, work( n*jvec+1 ), 1 )
336  CALL sgemv( trans, n, n, -bcoefi, b, lda, e( 1, jvec ),
337  \$ 1, one, work( n*jvec+1 ), 1 )
338  CALL sgemv( 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 = slange( '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 / ( REAL( n )*ulp )
378 *
379  RETURN
380 *
381 * End of SGET52
382 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:158
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:116
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69

Here is the call graph for this function:

Here is the caller graph for this function: