LAPACK 3.12.1
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:112
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:108
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: