 LAPACK  3.10.0 LAPACK: Linear Algebra PACKage

◆ 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:2942
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: