LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
sget22.f
Go to the documentation of this file.
1 *> \brief \b SGET22
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 SGET22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
12 * WI, WORK, RESULT )
13 *
14 * .. Scalar Arguments ..
15 * CHARACTER TRANSA, TRANSE, TRANSW
16 * INTEGER LDA, LDE, N
17 * ..
18 * .. Array Arguments ..
19 * REAL A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
20 * $ WORK( * ), WR( * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> SGET22 does an eigenvector check.
30 *>
31 *> The basic test is:
32 *>
33 *> RESULT(1) = | A E - E W | / ( |A| |E| ulp )
34 *>
35 *> using the 1-norm. It also tests the normalization of E:
36 *>
37 *> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
38 *> j
39 *>
40 *> where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
41 *> vector. If an eigenvector is complex, as determined from WI(j)
42 *> nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum
43 *> of
44 *> |er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)|
45 *>
46 *> W is a block diagonal matrix, with a 1 by 1 block for each real
47 *> eigenvalue and a 2 by 2 block for each complex conjugate pair.
48 *> If eigenvalues j and j+1 are a complex conjugate pair, so that
49 *> WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2
50 *> block corresponding to the pair will be:
51 *>
52 *> ( wr wi )
53 *> ( -wi wr )
54 *>
55 *> Such a block multiplying an n by 2 matrix ( ur ui ) on the right
56 *> will be the same as multiplying ur + i*ui by wr + i*wi.
57 *>
58 *> To handle various schemes for storage of left eigenvectors, there are
59 *> options to use A-transpose instead of A, E-transpose instead of E,
60 *> and/or W-transpose instead of W.
61 *> \endverbatim
62 *
63 * Arguments:
64 * ==========
65 *
66 *> \param[in] TRANSA
67 *> \verbatim
68 *> TRANSA is CHARACTER*1
69 *> Specifies whether or not A is transposed.
70 *> = 'N': No transpose
71 *> = 'T': Transpose
72 *> = 'C': Conjugate transpose (= Transpose)
73 *> \endverbatim
74 *>
75 *> \param[in] TRANSE
76 *> \verbatim
77 *> TRANSE is CHARACTER*1
78 *> Specifies whether or not E is transposed.
79 *> = 'N': No transpose, eigenvectors are in columns of E
80 *> = 'T': Transpose, eigenvectors are in rows of E
81 *> = 'C': Conjugate transpose (= Transpose)
82 *> \endverbatim
83 *>
84 *> \param[in] TRANSW
85 *> \verbatim
86 *> TRANSW is CHARACTER*1
87 *> Specifies whether or not W is transposed.
88 *> = 'N': No transpose
89 *> = 'T': Transpose, use -WI(j) instead of WI(j)
90 *> = 'C': Conjugate transpose, use -WI(j) instead of WI(j)
91 *> \endverbatim
92 *>
93 *> \param[in] N
94 *> \verbatim
95 *> N is INTEGER
96 *> The order of the matrix A. N >= 0.
97 *> \endverbatim
98 *>
99 *> \param[in] A
100 *> \verbatim
101 *> A is REAL array, dimension (LDA,N)
102 *> The matrix whose eigenvectors are in E.
103 *> \endverbatim
104 *>
105 *> \param[in] LDA
106 *> \verbatim
107 *> LDA is INTEGER
108 *> The leading dimension of the array A. LDA >= max(1,N).
109 *> \endverbatim
110 *>
111 *> \param[in] E
112 *> \verbatim
113 *> E is REAL array, dimension (LDE,N)
114 *> The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors
115 *> are stored in the columns of E, if TRANSE = 'T' or 'C', the
116 *> eigenvectors are stored in the rows of E.
117 *> \endverbatim
118 *>
119 *> \param[in] LDE
120 *> \verbatim
121 *> LDE is INTEGER
122 *> The leading dimension of the array E. LDE >= max(1,N).
123 *> \endverbatim
124 *>
125 *> \param[in] WR
126 *> \verbatim
127 *> WR is REAL array, dimension (N)
128 *> \endverbatim
129 *>
130 *> \param[in] WI
131 *> \verbatim
132 *> WI is REAL array, dimension (N)
133 *>
134 *> The real and imaginary parts of the eigenvalues of A.
135 *> Purely real eigenvalues are indicated by WI(j) = 0.
136 *> Complex conjugate pairs are indicated by WR(j)=WR(j+1) and
137 *> WI(j) = - WI(j+1) non-zero; the real part is assumed to be
138 *> stored in the j-th row/column and the imaginary part in
139 *> the (j+1)-th row/column.
140 *> \endverbatim
141 *>
142 *> \param[out] WORK
143 *> \verbatim
144 *> WORK is REAL array, dimension (N*(N+1))
145 *> \endverbatim
146 *>
147 *> \param[out] RESULT
148 *> \verbatim
149 *> RESULT is REAL array, dimension (2)
150 *> RESULT(1) = | A E - E W | / ( |A| |E| ulp )
151 *> RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
152 *> j
153 *> \endverbatim
154 *
155 * Authors:
156 * ========
157 *
158 *> \author Univ. of Tennessee
159 *> \author Univ. of California Berkeley
160 *> \author Univ. of Colorado Denver
161 *> \author NAG Ltd.
162 *
163 *> \ingroup single_eig
164 *
165 * =====================================================================
166  SUBROUTINE sget22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
167  $ WI, WORK, RESULT )
168 *
169 * -- LAPACK test routine --
170 * -- LAPACK is a software package provided by Univ. of Tennessee, --
171 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
172 *
173 * .. Scalar Arguments ..
174  CHARACTER TRANSA, TRANSE, TRANSW
175  INTEGER LDA, LDE, N
176 * ..
177 * .. Array Arguments ..
178  REAL A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
179  $ work( * ), wr( * )
180 * ..
181 *
182 * =====================================================================
183 *
184 * .. Parameters ..
185  REAL ZERO, ONE
186  parameter( zero = 0.0, one = 1.0 )
187 * ..
188 * .. Local Scalars ..
189  CHARACTER NORMA, NORME
190  INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
191  $ jvec
192  REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
193  $ ulp, unfl
194 * ..
195 * .. Local Arrays ..
196  REAL WMAT( 2, 2 )
197 * ..
198 * .. External Functions ..
199  LOGICAL LSAME
200  REAL SLAMCH, SLANGE
201  EXTERNAL lsame, slamch, slange
202 * ..
203 * .. External Subroutines ..
204  EXTERNAL saxpy, sgemm, slaset
205 * ..
206 * .. Intrinsic Functions ..
207  INTRINSIC abs, max, min, real
208 * ..
209 * .. Executable Statements ..
210 *
211 * Initialize RESULT (in case N=0)
212 *
213  result( 1 ) = zero
214  result( 2 ) = zero
215  IF( n.LE.0 )
216  $ RETURN
217 *
218  unfl = slamch( 'Safe minimum' )
219  ulp = slamch( 'Precision' )
220 *
221  itrnse = 0
222  ince = 1
223  norma = 'O'
224  norme = 'O'
225 *
226  IF( lsame( transa, 'T' ) .OR. lsame( transa, 'C' ) ) THEN
227  norma = 'I'
228  END IF
229  IF( lsame( transe, 'T' ) .OR. lsame( transe, 'C' ) ) THEN
230  norme = 'I'
231  itrnse = 1
232  ince = lde
233  END IF
234 *
235 * Check normalization of E
236 *
237  enrmin = one / ulp
238  enrmax = zero
239  IF( itrnse.EQ.0 ) THEN
240 *
241 * Eigenvectors are column vectors.
242 *
243  ipair = 0
244  DO 30 jvec = 1, n
245  temp1 = zero
246  IF( ipair.EQ.0 .AND. jvec.LT.n .AND. wi( jvec ).NE.zero )
247  $ ipair = 1
248  IF( ipair.EQ.1 ) THEN
249 *
250 * Complex eigenvector
251 *
252  DO 10 j = 1, n
253  temp1 = max( temp1, abs( e( j, jvec ) )+
254  $ abs( e( j, jvec+1 ) ) )
255  10 CONTINUE
256  enrmin = min( enrmin, temp1 )
257  enrmax = max( enrmax, temp1 )
258  ipair = 2
259  ELSE IF( ipair.EQ.2 ) THEN
260  ipair = 0
261  ELSE
262 *
263 * Real eigenvector
264 *
265  DO 20 j = 1, n
266  temp1 = max( temp1, abs( e( j, jvec ) ) )
267  20 CONTINUE
268  enrmin = min( enrmin, temp1 )
269  enrmax = max( enrmax, temp1 )
270  ipair = 0
271  END IF
272  30 CONTINUE
273 *
274  ELSE
275 *
276 * Eigenvectors are row vectors.
277 *
278  DO 40 jvec = 1, n
279  work( jvec ) = zero
280  40 CONTINUE
281 *
282  DO 60 j = 1, n
283  ipair = 0
284  DO 50 jvec = 1, n
285  IF( ipair.EQ.0 .AND. jvec.LT.n .AND. wi( jvec ).NE.zero )
286  $ ipair = 1
287  IF( ipair.EQ.1 ) THEN
288  work( jvec ) = max( work( jvec ),
289  $ abs( e( j, jvec ) )+abs( e( j,
290  $ jvec+1 ) ) )
291  work( jvec+1 ) = work( jvec )
292  ELSE IF( ipair.EQ.2 ) THEN
293  ipair = 0
294  ELSE
295  work( jvec ) = max( work( jvec ),
296  $ abs( e( j, jvec ) ) )
297  ipair = 0
298  END IF
299  50 CONTINUE
300  60 CONTINUE
301 *
302  DO 70 jvec = 1, n
303  enrmin = min( enrmin, work( jvec ) )
304  enrmax = max( enrmax, work( jvec ) )
305  70 CONTINUE
306  END IF
307 *
308 * Norm of A:
309 *
310  anorm = max( slange( norma, n, n, a, lda, work ), unfl )
311 *
312 * Norm of E:
313 *
314  enorm = max( slange( norme, n, n, e, lde, work ), ulp )
315 *
316 * Norm of error:
317 *
318 * Error = AE - EW
319 *
320  CALL slaset( 'Full', n, n, zero, zero, work, n )
321 *
322  ipair = 0
323  ierow = 1
324  iecol = 1
325 *
326  DO 80 jcol = 1, n
327  IF( itrnse.EQ.1 ) THEN
328  ierow = jcol
329  ELSE
330  iecol = jcol
331  END IF
332 *
333  IF( ipair.EQ.0 .AND. wi( jcol ).NE.zero )
334  $ ipair = 1
335 *
336  IF( ipair.EQ.1 ) THEN
337  wmat( 1, 1 ) = wr( jcol )
338  wmat( 2, 1 ) = -wi( jcol )
339  wmat( 1, 2 ) = wi( jcol )
340  wmat( 2, 2 ) = wr( jcol )
341  CALL sgemm( transe, transw, n, 2, 2, one, e( ierow, iecol ),
342  $ lde, wmat, 2, zero, work( n*( jcol-1 )+1 ), n )
343  ipair = 2
344  ELSE IF( ipair.EQ.2 ) THEN
345  ipair = 0
346 *
347  ELSE
348 *
349  CALL saxpy( n, wr( jcol ), e( ierow, iecol ), ince,
350  $ work( n*( jcol-1 )+1 ), 1 )
351  ipair = 0
352  END IF
353 *
354  80 CONTINUE
355 *
356  CALL sgemm( transa, transe, n, n, n, one, a, lda, e, lde, -one,
357  $ work, n )
358 *
359  errnrm = slange( 'One', n, n, work, n, work( n*n+1 ) ) / enorm
360 *
361 * Compute RESULT(1) (avoiding under/overflow)
362 *
363  IF( anorm.GT.errnrm ) THEN
364  result( 1 ) = ( errnrm / anorm ) / ulp
365  ELSE
366  IF( anorm.LT.one ) THEN
367  result( 1 ) = one / ulp
368  ELSE
369  result( 1 ) = min( errnrm / anorm, one ) / ulp
370  END IF
371  END IF
372 *
373 * Compute RESULT(2) : the normalization error in E.
374 *
375  result( 2 ) = max( abs( enrmax-one ), abs( enrmin-one ) ) /
376  $ ( real( n )*ulp )
377 *
378  RETURN
379 *
380 * End of SGET22
381 *
382  END
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: slaset.f:110
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
subroutine sget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
SGET22
Definition: sget22.f:168