LAPACK  3.9.0
LAPACK: Linear Algebra PACKage
dget22.f
Go to the documentation of this file.
1 *> \brief \b DGET22
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 DGET22( 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 * DOUBLE PRECISION A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
20 * $ WORK( * ), WR( * )
21 * ..
22 *
23 *
24 *> \par Purpose:
25 * =============
26 *>
27 *> \verbatim
28 *>
29 *> DGET22 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N)
128 *> \endverbatim
129 *>
130 *> \param[in] WI
131 *> \verbatim
132 *> WI is DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N*(N+1))
145 *> \endverbatim
146 *>
147 *> \param[out] RESULT
148 *> \verbatim
149 *> RESULT is DOUBLE PRECISION 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 *> \endverbatim
153 *
154 * Authors:
155 * ========
156 *
157 *> \author Univ. of Tennessee
158 *> \author Univ. of California Berkeley
159 *> \author Univ. of Colorado Denver
160 *> \author NAG Ltd.
161 *
162 *> \date December 2016
163 *
164 *> \ingroup double_eig
165 *
166 * =====================================================================
167  SUBROUTINE dget22( TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR,
168  $ WI, WORK, RESULT )
169 *
170 * -- LAPACK test routine (version 3.7.0) --
171 * -- LAPACK is a software package provided by Univ. of Tennessee, --
172 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
173 * December 2016
174 *
175 * .. Scalar Arguments ..
176  CHARACTER TRANSA, TRANSE, TRANSW
177  INTEGER LDA, LDE, N
178 * ..
179 * .. Array Arguments ..
180  DOUBLE PRECISION A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
181  $ work( * ), wr( * )
182 * ..
183 *
184 * =====================================================================
185 *
186 * .. Parameters ..
187  DOUBLE PRECISION ZERO, ONE
188  parameter( zero = 0.0d0, one = 1.0d0 )
189 * ..
190 * .. Local Scalars ..
191  CHARACTER NORMA, NORME
192  INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
193  $ jvec
194  DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
195  $ ulp, unfl
196 * ..
197 * .. Local Arrays ..
198  DOUBLE PRECISION WMAT( 2, 2 )
199 * ..
200 * .. External Functions ..
201  LOGICAL LSAME
202  DOUBLE PRECISION DLAMCH, DLANGE
203  EXTERNAL lsame, dlamch, dlange
204 * ..
205 * .. External Subroutines ..
206  EXTERNAL daxpy, dgemm, dlaset
207 * ..
208 * .. Intrinsic Functions ..
209  INTRINSIC abs, dble, max, min
210 * ..
211 * .. Executable Statements ..
212 *
213 * Initialize RESULT (in case N=0)
214 *
215  result( 1 ) = zero
216  result( 2 ) = zero
217  IF( n.LE.0 )
218  $ RETURN
219 *
220  unfl = dlamch( 'Safe minimum' )
221  ulp = dlamch( 'Precision' )
222 *
223  itrnse = 0
224  ince = 1
225  norma = 'O'
226  norme = 'O'
227 *
228  IF( lsame( transa, 'T' ) .OR. lsame( transa, 'C' ) ) THEN
229  norma = 'I'
230  END IF
231  IF( lsame( transe, 'T' ) .OR. lsame( transe, 'C' ) ) THEN
232  norme = 'I'
233  itrnse = 1
234  ince = lde
235  END IF
236 *
237 * Check normalization of E
238 *
239  enrmin = one / ulp
240  enrmax = zero
241  IF( itrnse.EQ.0 ) THEN
242 *
243 * Eigenvectors are column vectors.
244 *
245  ipair = 0
246  DO 30 jvec = 1, n
247  temp1 = zero
248  IF( ipair.EQ.0 .AND. jvec.LT.n .AND. wi( jvec ).NE.zero )
249  $ ipair = 1
250  IF( ipair.EQ.1 ) THEN
251 *
252 * Complex eigenvector
253 *
254  DO 10 j = 1, n
255  temp1 = max( temp1, abs( e( j, jvec ) )+
256  $ abs( e( j, jvec+1 ) ) )
257  10 CONTINUE
258  enrmin = min( enrmin, temp1 )
259  enrmax = max( enrmax, temp1 )
260  ipair = 2
261  ELSE IF( ipair.EQ.2 ) THEN
262  ipair = 0
263  ELSE
264 *
265 * Real eigenvector
266 *
267  DO 20 j = 1, n
268  temp1 = max( temp1, abs( e( j, jvec ) ) )
269  20 CONTINUE
270  enrmin = min( enrmin, temp1 )
271  enrmax = max( enrmax, temp1 )
272  ipair = 0
273  END IF
274  30 CONTINUE
275 *
276  ELSE
277 *
278 * Eigenvectors are row vectors.
279 *
280  DO 40 jvec = 1, n
281  work( jvec ) = zero
282  40 CONTINUE
283 *
284  DO 60 j = 1, n
285  ipair = 0
286  DO 50 jvec = 1, n
287  IF( ipair.EQ.0 .AND. jvec.LT.n .AND. wi( jvec ).NE.zero )
288  $ ipair = 1
289  IF( ipair.EQ.1 ) THEN
290  work( jvec ) = max( work( jvec ),
291  $ abs( e( j, jvec ) )+abs( e( j,
292  $ jvec+1 ) ) )
293  work( jvec+1 ) = work( jvec )
294  ELSE IF( ipair.EQ.2 ) THEN
295  ipair = 0
296  ELSE
297  work( jvec ) = max( work( jvec ),
298  $ abs( e( j, jvec ) ) )
299  ipair = 0
300  END IF
301  50 CONTINUE
302  60 CONTINUE
303 *
304  DO 70 jvec = 1, n
305  enrmin = min( enrmin, work( jvec ) )
306  enrmax = max( enrmax, work( jvec ) )
307  70 CONTINUE
308  END IF
309 *
310 * Norm of A:
311 *
312  anorm = max( dlange( norma, n, n, a, lda, work ), unfl )
313 *
314 * Norm of E:
315 *
316  enorm = max( dlange( norme, n, n, e, lde, work ), ulp )
317 *
318 * Norm of error:
319 *
320 * Error = AE - EW
321 *
322  CALL dlaset( 'Full', n, n, zero, zero, work, n )
323 *
324  ipair = 0
325  ierow = 1
326  iecol = 1
327 *
328  DO 80 jcol = 1, n
329  IF( itrnse.EQ.1 ) THEN
330  ierow = jcol
331  ELSE
332  iecol = jcol
333  END IF
334 *
335  IF( ipair.EQ.0 .AND. wi( jcol ).NE.zero )
336  $ ipair = 1
337 *
338  IF( ipair.EQ.1 ) THEN
339  wmat( 1, 1 ) = wr( jcol )
340  wmat( 2, 1 ) = -wi( jcol )
341  wmat( 1, 2 ) = wi( jcol )
342  wmat( 2, 2 ) = wr( jcol )
343  CALL dgemm( transe, transw, n, 2, 2, one, e( ierow, iecol ),
344  $ lde, wmat, 2, zero, work( n*( jcol-1 )+1 ), n )
345  ipair = 2
346  ELSE IF( ipair.EQ.2 ) THEN
347  ipair = 0
348 *
349  ELSE
350 *
351  CALL daxpy( n, wr( jcol ), e( ierow, iecol ), ince,
352  $ work( n*( jcol-1 )+1 ), 1 )
353  ipair = 0
354  END IF
355 *
356  80 CONTINUE
357 *
358  CALL dgemm( transa, transe, n, n, n, one, a, lda, e, lde, -one,
359  $ work, n )
360 *
361  errnrm = dlange( 'One', n, n, work, n, work( n*n+1 ) ) / enorm
362 *
363 * Compute RESULT(1) (avoiding under/overflow)
364 *
365  IF( anorm.GT.errnrm ) THEN
366  result( 1 ) = ( errnrm / anorm ) / ulp
367  ELSE
368  IF( anorm.LT.one ) THEN
369  result( 1 ) = ( min( errnrm, anorm ) / anorm ) / ulp
370  ELSE
371  result( 1 ) = min( errnrm / anorm, one ) / ulp
372  END IF
373  END IF
374 *
375 * Compute RESULT(2) : the normalization error in E.
376 *
377  result( 2 ) = max( abs( enrmax-one ), abs( enrmin-one ) ) /
378  $ ( dble( n )*ulp )
379 *
380  RETURN
381 *
382 * End of DGET22
383 *
384  END
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:189
dget22
subroutine dget22(TRANSA, TRANSE, TRANSW, N, A, LDA, E, LDE, WR, WI, WORK, RESULT)
DGET22
Definition: dget22.f:169
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:112
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:91