LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
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
 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 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
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

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*
lde
logical function lde(ri, rj, lr)
Definition dblat2.f:2970
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
dlamch
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
dlange
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
dlaset
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
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
Here is the call graph for this function:
Here is the caller graph for this function: