LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sget52()

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.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 197 of file sget52.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  REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
210  $ B( LDB, * ), BETA( * ), E( LDE, * ),
211  $ RESULT( 2 ), WORK( * )
212 * ..
213 *
214 * =====================================================================
215 *
216 * .. Parameters ..
217  REAL ZERO, ONE, TEN
218  parameter( zero = 0.0, one = 1.0, ten = 10.0 )
219 * ..
220 * .. Local Scalars ..
221  LOGICAL ILCPLX
222  CHARACTER NORMAB, TRANS
223  INTEGER J, JVEC
224  REAL 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  REAL SLAMCH, SLANGE
230  EXTERNAL slamch, slange
231 * ..
232 * .. External Subroutines ..
233  EXTERNAL sgemv
234 * ..
235 * .. Intrinsic Functions ..
236  INTRINSIC abs, max, real
237 * ..
238 * .. Executable Statements ..
239 *
240  result( 1 ) = zero
241  result( 2 ) = zero
242  IF( n.LE.0 )
243  $ RETURN
244 *
245  safmin = slamch( 'Safe minimum' )
246  safmax = one / safmin
247  ulp = slamch( 'Epsilon' )*slamch( '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( slange( normab, n, n, a, lda, work ), safmin )
260  bnorm = max( slange( normab, n, n, b, ldb, work ), safmin )
261  enorm = max( slange( '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 sgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
295  $ zero, work( n*( jvec-1 )+1 ), 1 )
296  CALL sgemv( 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 sgemv( trans, n, n, acoef, a, lda, e( 1, jvec ), 1,
325  $ zero, work( n*( jvec-1 )+1 ), 1 )
326  CALL sgemv( trans, n, n, -bcoefr, b, lda, e( 1, jvec ),
327  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
328  CALL sgemv( trans, n, n, bcoefi, b, lda, e( 1, jvec+1 ),
329  $ 1, one, work( n*( jvec-1 )+1 ), 1 )
330 *
331  CALL sgemv( trans, n, n, acoef, a, lda, e( 1, jvec+1 ),
332  $ 1, zero, work( n*jvec+1 ), 1 )
333  CALL sgemv( trans, n, n, -bcoefi, b, lda, e( 1, jvec ),
334  $ 1, one, work( n*jvec+1 ), 1 )
335  CALL sgemv( 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 = slange( '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 / ( real( n )*ulp )
375 *
376  RETURN
377 *
378 * End of SGET52
379 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2942
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:114
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: