LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> 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 double_eig
164*
165* =====================================================================
166 SUBROUTINE dget22( 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 DOUBLE PRECISION A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
179 $ work( * ), wr( * )
180* ..
181*
182* =====================================================================
183*
184* .. Parameters ..
185 DOUBLE PRECISION ZERO, ONE
186 parameter( zero = 0.0d0, one = 1.0d0 )
187* ..
188* .. Local Scalars ..
189 CHARACTER NORMA, NORME
190 INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
191 $ jvec
192 DOUBLE PRECISION ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
193 $ ulp, unfl
194* ..
195* .. Local Arrays ..
196 DOUBLE PRECISION WMAT( 2, 2 )
197* ..
198* .. External Functions ..
199 LOGICAL LSAME
200 DOUBLE PRECISION DLAMCH, DLANGE
201 EXTERNAL lsame, dlamch, dlange
202* ..
203* .. External Subroutines ..
204 EXTERNAL daxpy, dgemm, dlaset
205* ..
206* .. Intrinsic Functions ..
207 INTRINSIC abs, dble, max, min
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 = dlamch( 'Safe minimum' )
219 ulp = dlamch( '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( dlange( norma, n, n, a, lda, work ), unfl )
311*
312* Norm of E:
313*
314 enorm = max( dlange( norme, n, n, e, lde, work ), ulp )
315*
316* Norm of error:
317*
318* Error = AE - EW
319*
320 CALL dlaset( '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 dgemm( 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 daxpy( 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 dgemm( transa, transe, n, n, n, one, a, lda, e, lde, -one,
357 $ work, n )
358*
359 errnrm = dlange( '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 $ ( dble( n )*ulp )
377*
378 RETURN
379*
380* End of DGET22
381*
382 END
subroutine dget22(transa, transe, transw, n, a, lda, e, lde, wr, wi, work, result)
DGET22
Definition dget22.f:168
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
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:110