 LAPACK 3.11.0 LAPACK: Linear Algebra PACKage
Searching...
No Matches

## ◆ dget22()

 subroutine dget22 ( character TRANSA, character TRANSE, character TRANSW, integer N, double precision, dimension( lda, * ) A, integer LDA, double precision, dimension( lde, * ) E, integer LDE, double precision, dimension( * ) WR, double precision, dimension( * ) WI, double precision, dimension( * ) WORK, double precision, dimension( 2 ) RESULT )

DGET22

Purpose:
``` DGET22 does an eigenvector check.

The basic test is:

RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )

using the 1-norm.  It also tests the normalization of E:

RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
j

where E(j) is the j-th eigenvector, and m-norm is the max-norm of a
vector.  If an eigenvector is complex, as determined from WI(j)
nonzero, then the max-norm of the vector ( er + i*ei ) is the maximum
of
|er(1)| + |ei(1)|, ... , |er(n)| + |ei(n)|

W is a block diagonal matrix, with a 1 by 1 block for each real
eigenvalue and a 2 by 2 block for each complex conjugate pair.
If eigenvalues j and j+1 are a complex conjugate pair, so that
WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the 2 by 2
block corresponding to the pair will be:

(  wr  wi  )
( -wi  wr  )

Such a block multiplying an n by 2 matrix ( ur ui ) on the right
will be the same as multiplying  ur + i*ui  by  wr + i*wi.

To handle various schemes for storage of left eigenvectors, there are
Parameters
 [in] TRANSA ``` TRANSA is CHARACTER*1 Specifies whether or not A is transposed. = 'N': No transpose = 'T': Transpose = 'C': Conjugate transpose (= Transpose)``` [in] TRANSE ``` TRANSE is CHARACTER*1 Specifies whether or not E is transposed. = 'N': No transpose, eigenvectors are in columns of E = 'T': Transpose, eigenvectors are in rows of E = 'C': Conjugate transpose (= Transpose)``` [in] TRANSW ``` TRANSW is CHARACTER*1 Specifies whether or not W is transposed. = 'N': No transpose = 'T': Transpose, use -WI(j) instead of WI(j) = 'C': Conjugate transpose, use -WI(j) instead of WI(j)``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in] A ``` A is DOUBLE PRECISION array, dimension (LDA,N) The matrix whose eigenvectors are in E.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [in] E ``` E is DOUBLE PRECISION array, dimension (LDE,N) The matrix of eigenvectors. If TRANSE = 'N', the eigenvectors are stored in the columns of E, if TRANSE = 'T' or 'C', the eigenvectors are stored in the rows of E.``` [in] LDE ``` LDE is INTEGER The leading dimension of the array E. LDE >= max(1,N).``` [in] WR ` WR is DOUBLE PRECISION array, dimension (N)` [in] WI ``` WI is DOUBLE PRECISION array, dimension (N) The real and imaginary parts of the eigenvalues of A. Purely real eigenvalues are indicated by WI(j) = 0. Complex conjugate pairs are indicated by WR(j)=WR(j+1) and WI(j) = - WI(j+1) non-zero; the real part is assumed to be stored in the j-th row/column and the imaginary part in the (j+1)-th row/column.``` [out] WORK ` WORK is DOUBLE PRECISION array, dimension (N*(N+1))` [out] RESULT ``` RESULT is DOUBLE PRECISION array, dimension (2) RESULT(1) = | A E - E W | / ( |A| |E| ulp ) RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp ) j```

Definition at line 166 of file dget22.f.

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*
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2970
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
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
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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:187
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: