LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ sget22()

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 )
                       j
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 166 of file sget22.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  REAL A( LDA, * ), E( LDE, * ), RESULT( 2 ), WI( * ),
179  $ WORK( * ), WR( * )
180 * ..
181 *
182 * =====================================================================
183 *
184 * .. Parameters ..
185  REAL ZERO, ONE
186  parameter( zero = 0.0, one = 1.0 )
187 * ..
188 * .. Local Scalars ..
189  CHARACTER NORMA, NORME
190  INTEGER IECOL, IEROW, INCE, IPAIR, ITRNSE, J, JCOL,
191  $ JVEC
192  REAL ANORM, ENORM, ENRMAX, ENRMIN, ERRNRM, TEMP1,
193  $ ULP, UNFL
194 * ..
195 * .. Local Arrays ..
196  REAL WMAT( 2, 2 )
197 * ..
198 * .. External Functions ..
199  LOGICAL LSAME
200  REAL SLAMCH, SLANGE
201  EXTERNAL lsame, slamch, slange
202 * ..
203 * .. External Subroutines ..
204  EXTERNAL saxpy, sgemm, slaset
205 * ..
206 * .. Intrinsic Functions ..
207  INTRINSIC abs, max, min, real
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 = slamch( 'Safe minimum' )
219  ulp = slamch( '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( slange( norma, n, n, a, lda, work ), unfl )
311 *
312 * Norm of E:
313 *
314  enorm = max( slange( norme, n, n, e, lde, work ), ulp )
315 *
316 * Norm of error:
317 *
318 * Error = AE - EW
319 *
320  CALL slaset( '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 sgemm( 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 saxpy( 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 sgemm( transa, transe, n, n, n, one, a, lda, e, lde, -one,
357  $ work, n )
358 *
359  errnrm = slange( '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  $ ( real( n )*ulp )
377 *
378  RETURN
379 *
380 * End of SGET22
381 *
logical function lde(RI, RJ, LR)
Definition: dblat2.f:2942
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:110
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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:114
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
Definition: saxpy.f:89
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:187
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: