LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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.
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:2970
subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
DGEMV
Definition dgemv.f:158
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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: