LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine sget22 ( character  TRANSA,
character  TRANSE,
character  TRANSW,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( lde, * )  E,
integer  LDE,
real, dimension( * )  WR,
real, dimension( * )  WI,
real, dimension( * )  WORK,
real, dimension( 2 )  RESULT 
)

SGET22

Purpose: 
 SGET22 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
 options to use A-transpose instead of A, E-transpose instead of E,
 and/or W-transpose instead of W.
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 REAL 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 REAL 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 REAL array, dimension (N)
[in]WI
          WI is REAL 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 REAL array, dimension (N*(N+1))
[out]RESULT
          RESULT is REAL array, dimension (2)
          RESULT(1) = | A E  -  E W | / ( |A| |E| ulp )
          RESULT(2) = max | m-norm(E(j)) - 1 | / ( n ulp )
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Date
November 2011

Definition at line 169 of file sget22.f.

169 *
170 * -- LAPACK test routine (version 3.4.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 * November 2011
174 *
175 * .. Scalar Arguments ..
176  CHARACTER transa, transe, transw
177  INTEGER lda, lde, n
178 * ..
179 * .. Array Arguments ..
180  REAL a( lda, * ), e( lde, * ), result( 2 ), wi( * ),
181  $ work( * ), wr( * )
182 * ..
183 *
184 * =====================================================================
185 *
186 * .. Parameters ..
187  REAL zero, one
188  parameter ( zero = 0.0, one = 1.0 )
189 * ..
190 * .. Local Scalars ..
191  CHARACTER norma, norme
192  INTEGER iecol, ierow, ince, ipair, itrnse, j, jcol,
193  $ jvec
194  REAL anorm, enorm, enrmax, enrmin, errnrm, temp1,
195  $ ulp, unfl
196 * ..
197 * .. Local Arrays ..
198  REAL wmat( 2, 2 )
199 * ..
200 * .. External Functions ..
201  LOGICAL lsame
202  REAL slamch, slange
203  EXTERNAL lsame, slamch, slange
204 * ..
205 * .. External Subroutines ..
206  EXTERNAL saxpy, sgemm, slaset
207 * ..
208 * .. Intrinsic Functions ..
209  INTRINSIC abs, max, min, real
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 = slamch( 'Safe minimum' )
221  ulp = slamch( '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( slange( norma, n, n, a, lda, work ), unfl )
313 *
314 * Norm of E:
315 *
316  enorm = max( slange( norme, n, n, e, lde, work ), ulp )
317 *
318 * Norm of error:
319 *
320 * Error = AE - EW
321 *
322  CALL slaset( '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 sgemm( 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 saxpy( 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 sgemm( transa, transe, n, n, n, one, a, lda, e, lde, -one,
359  $ work, n )
360 *
361  errnrm = slange( '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  $ ( REAL( n )*ulp )
379 *
380  RETURN
381 *
382 * End of SGET22
383 *
lde
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2945
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
slaset
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:112
slange
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:116
saxpy
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:54
slamch
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:69
lsame
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:55

Here is the call graph for this function:

Here is the caller graph for this function: