LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
sget52.f
Go to the documentation of this file.
1 *> \brief \b SGET52
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE SGET52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
12 * ALPHAI, BETA, WORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * LOGICAL LEFT
16 * INTEGER LDA, LDB, LDE, N
17 * ..
18 * .. Array Arguments ..
19 * REAL A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
20 * $ B( LDB, * ), BETA( * ), E( LDE, * ),
21 * $ RESULT( 2 ), WORK( * )
22 * ..
23 *
24 *
25 *> \par Purpose:
26 * =============
27 *>
28 *> \verbatim
29 *>
30 *> SGET52 does an eigenvector check for the generalized eigenvalue
31 *> problem.
32 *>
33 *> The basic test for right eigenvectors is:
34 *>
35 *> | b(j) A E(j) - a(j) B E(j) |
36 *> RESULT(1) = max -------------------------------
37 *> j n ulp max( |b(j) A|, |a(j) B| )
38 *>
39 *> using the 1-norm. Here, a(j)/b(j) = w is the j-th generalized
40 *> eigenvalue of A - w B, or, equivalently, b(j)/a(j) = m is the j-th
41 *> generalized eigenvalue of m A - B.
42 *>
43 *> For real eigenvalues, the test is straightforward. For complex
44 *> eigenvalues, E(j) and a(j) are complex, represented by
45 *> Er(j) + i*Ei(j) and ar(j) + i*ai(j), resp., so the test for that
46 *> eigenvector becomes
47 *>
48 *> max( |Wr|, |Wi| )
49 *> --------------------------------------------
50 *> n ulp max( |b(j) A|, (|ar(j)|+|ai(j)|) |B| )
51 *>
52 *> where
53 *>
54 *> Wr = b(j) A Er(j) - ar(j) B Er(j) + ai(j) B Ei(j)
55 *>
56 *> Wi = b(j) A Ei(j) - ai(j) B Er(j) - ar(j) B Ei(j)
57 *>
58 *> T T _
59 *> For left eigenvectors, A , B , a, and b are used.
60 *>
61 *> SGET52 also tests the normalization of E. Each eigenvector is
62 *> supposed to be normalized so that the maximum "absolute value"
63 *> of its elements is 1, where in this case, "absolute value"
64 *> of a complex value x is |Re(x)| + |Im(x)| ; let us call this
65 *> maximum "absolute value" norm of a vector v M(v).
66 *> if a(j)=b(j)=0, then the eigenvector is set to be the jth coordinate
67 *> vector. The normalization test is:
68 *>
69 *> RESULT(2) = max | M(v(j)) - 1 | / ( n ulp )
70 *> eigenvectors v(j)
71 *> \endverbatim
72 *
73 * Arguments:
74 * ==========
75 *
76 *> \param[in] LEFT
77 *> \verbatim
78 *> LEFT is LOGICAL
79 *> =.TRUE.: The eigenvectors in the columns of E are assumed
80 *> to be *left* eigenvectors.
81 *> =.FALSE.: The eigenvectors in the columns of E are assumed
82 *> to be *right* eigenvectors.
83 *> \endverbatim
84 *>
85 *> \param[in] N
86 *> \verbatim
87 *> N is INTEGER
88 *> The size of the matrices. If it is zero, SGET52 does
89 *> nothing. It must be at least zero.
90 *> \endverbatim
91 *>
92 *> \param[in] A
93 *> \verbatim
94 *> A is REAL array, dimension (LDA, N)
95 *> The matrix A.
96 *> \endverbatim
97 *>
98 *> \param[in] LDA
99 *> \verbatim
100 *> LDA is INTEGER
101 *> The leading dimension of A. It must be at least 1
102 *> and at least N.
103 *> \endverbatim
104 *>
105 *> \param[in] B
106 *> \verbatim
107 *> B is REAL array, dimension (LDB, N)
108 *> The matrix B.
109 *> \endverbatim
110 *>
111 *> \param[in] LDB
112 *> \verbatim
113 *> LDB is INTEGER
114 *> The leading dimension of B. It must be at least 1
115 *> and at least N.
116 *> \endverbatim
117 *>
118 *> \param[in] E
119 *> \verbatim
120 *> E is REAL array, dimension (LDE, N)
121 *> The matrix of eigenvectors. It must be O( 1 ). Complex
122 *> eigenvalues and eigenvectors always come in pairs, the
123 *> eigenvalue and its conjugate being stored in adjacent
124 *> elements of ALPHAR, ALPHAI, and BETA. Thus, if a(j)/b(j)
125 *> and a(j+1)/b(j+1) are a complex conjugate pair of
126 *> generalized eigenvalues, then E(,j) contains the real part
127 *> of the eigenvector and E(,j+1) contains the imaginary part.
128 *> Note that whether E(,j) is a real eigenvector or part of a
129 *> complex one is specified by whether ALPHAI(j) is zero or not.
130 *> \endverbatim
131 *>
132 *> \param[in] LDE
133 *> \verbatim
134 *> LDE is INTEGER
135 *> The leading dimension of E. It must be at least 1 and at
136 *> least N.
137 *> \endverbatim
138 *>
139 *> \param[in] ALPHAR
140 *> \verbatim
141 *> ALPHAR is REAL array, dimension (N)
142 *> The real parts of the values a(j) as described above, which,
143 *> along with b(j), define the generalized eigenvalues.
144 *> Complex eigenvalues always come in complex conjugate pairs
145 *> a(j)/b(j) and a(j+1)/b(j+1), which are stored in adjacent
146 *> elements in ALPHAR, ALPHAI, and BETA. Thus, if the j-th
147 *> and (j+1)-st eigenvalues form a pair, ALPHAR(j+1)/BETA(j+1)
148 *> is assumed to be equal to ALPHAR(j)/BETA(j).
149 *> \endverbatim
150 *>
151 *> \param[in] ALPHAI
152 *> \verbatim
153 *> ALPHAI is REAL array, dimension (N)
154 *> The imaginary parts of the values a(j) as described above,
155 *> which, along with b(j), define the generalized eigenvalues.
156 *> If ALPHAI(j)=0, then the eigenvalue is real, otherwise it
157 *> is part of a complex conjugate pair. Complex eigenvalues
158 *> always come in complex conjugate pairs a(j)/b(j) and
159 *> a(j+1)/b(j+1), which are stored in adjacent elements in
160 *> ALPHAR, ALPHAI, and BETA. Thus, if the j-th and (j+1)-st
161 *> eigenvalues form a pair, ALPHAI(j+1)/BETA(j+1) is assumed to
162 *> be equal to -ALPHAI(j)/BETA(j). Also, nonzero values in
163 *> ALPHAI are assumed to always come in adjacent pairs.
164 *> \endverbatim
165 *>
166 *> \param[in] BETA
167 *> \verbatim
168 *> BETA is REAL array, dimension (N)
169 *> The values b(j) as described above, which, along with a(j),
170 *> define the generalized eigenvalues.
171 *> \endverbatim
172 *>
173 *> \param[out] WORK
174 *> \verbatim
175 *> WORK is REAL array, dimension (N**2+N)
176 *> \endverbatim
177 *>
178 *> \param[out] RESULT
179 *> \verbatim
180 *> RESULT is REAL array, dimension (2)
181 *> The values computed by the test described above. If A E or
182 *> B E is likely to overflow, then RESULT(1:2) is set to
183 *> 10 / ulp.
184 *> \endverbatim
185 *
186 * Authors:
187 * ========
188 *
189 *> \author Univ. of Tennessee
190 *> \author Univ. of California Berkeley
191 *> \author Univ. of Colorado Denver
192 *> \author NAG Ltd.
193 *
194 *> \ingroup single_eig
195 *
196 * =====================================================================
197  SUBROUTINE sget52( LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR,
198  $ ALPHAI, BETA, WORK, RESULT )
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 *
380  END
subroutine sgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
SGEMV
Definition: sgemv.f:156
subroutine sget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
SGET52
Definition: sget52.f:199