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

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:2942
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: