LAPACK 3.11.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dgeev()

subroutine dgeev ( character  JOBVL,
character  JOBVR,
integer  N,
double precision, dimension( lda, * )  A,
integer  LDA,
double precision, dimension( * )  WR,
double precision, dimension( * )  WI,
double precision, dimension( ldvl, * )  VL,
integer  LDVL,
double precision, dimension( ldvr, * )  VR,
integer  LDVR,
double precision, dimension( * )  WORK,
integer  LWORK,
integer  INFO 
)

DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices

Download DGEEV + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 DGEEV computes for an N-by-N real nonsymmetric matrix A, the
 eigenvalues and, optionally, the left and/or right eigenvectors.

 The right eigenvector v(j) of A satisfies
                  A * v(j) = lambda(j) * v(j)
 where lambda(j) is its eigenvalue.
 The left eigenvector u(j) of A satisfies
               u(j)**H * A = lambda(j) * u(j)**H
 where u(j)**H denotes the conjugate-transpose of u(j).

 The computed eigenvectors are normalized to have Euclidean norm
 equal to 1 and largest component real.
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N': left eigenvectors of A are not computed;
          = 'V': left eigenvectors of A are computed.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N': right eigenvectors of A are not computed;
          = 'V': right eigenvectors of A are computed.
[in]N
          N is INTEGER
          The order of the matrix A. N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA,N)
          On entry, the N-by-N matrix A.
          On exit, A has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A.  LDA >= max(1,N).
[out]WR
          WR is DOUBLE PRECISION array, dimension (N)
[out]WI
          WI is DOUBLE PRECISION array, dimension (N)
          WR and WI contain the real and imaginary parts,
          respectively, of the computed eigenvalues.  Complex
          conjugate pairs of eigenvalues appear consecutively
          with the eigenvalue having the positive imaginary part
          first.
[out]VL
          VL is DOUBLE PRECISION array, dimension (LDVL,N)
          If JOBVL = 'V', the left eigenvectors u(j) are stored one
          after another in the columns of VL, in the same order
          as their eigenvalues.
          If JOBVL = 'N', VL is not referenced.
          If the j-th eigenvalue is real, then u(j) = VL(:,j),
          the j-th column of VL.
          If the j-th and (j+1)-st eigenvalues form a complex
          conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
          u(j+1) = VL(:,j) - i*VL(:,j+1).
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the array VL.  LDVL >= 1; if
          JOBVL = 'V', LDVL >= N.
[out]VR
          VR is DOUBLE PRECISION array, dimension (LDVR,N)
          If JOBVR = 'V', the right eigenvectors v(j) are stored one
          after another in the columns of VR, in the same order
          as their eigenvalues.
          If JOBVR = 'N', VR is not referenced.
          If the j-th eigenvalue is real, then v(j) = VR(:,j),
          the j-th column of VR.
          If the j-th and (j+1)-st eigenvalues form a complex
          conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
          v(j+1) = VR(:,j) - i*VR(:,j+1).
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the array VR.  LDVR >= 1; if
          JOBVR = 'V', LDVR >= N.
[out]WORK
          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
[in]LWORK
          LWORK is INTEGER
          The dimension of the array WORK.  LWORK >= max(1,3*N), and
          if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N.  For good
          performance, LWORK must generally be larger.

          If LWORK = -1, then a workspace query is assumed; the routine
          only calculates the optimal size of the WORK array, returns
          this value as the first entry of the WORK array, and no error
          message related to LWORK is issued by XERBLA.
[out]INFO
          INFO is INTEGER
          = 0:  successful exit
          < 0:  if INFO = -i, the i-th argument had an illegal value.
          > 0:  if INFO = i, the QR algorithm failed to compute all the
                eigenvalues, and no eigenvectors have been computed;
                elements i+1:N of WR and WI contain eigenvalues which
                have converged.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 190 of file dgeev.f.

192 implicit none
193*
194* -- LAPACK driver routine --
195* -- LAPACK is a software package provided by Univ. of Tennessee, --
196* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
197*
198* .. Scalar Arguments ..
199 CHARACTER JOBVL, JOBVR
200 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N
201* ..
202* .. Array Arguments ..
203 DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
204 $ WI( * ), WORK( * ), WR( * )
205* ..
206*
207* =====================================================================
208*
209* .. Parameters ..
210 DOUBLE PRECISION ZERO, ONE
211 parameter( zero = 0.0d0, one = 1.0d0 )
212* ..
213* .. Local Scalars ..
214 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR
215 CHARACTER SIDE
216 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K,
217 $ LWORK_TREVC, MAXWRK, MINWRK, NOUT
218 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM,
219 $ SN
220* ..
221* .. Local Arrays ..
222 LOGICAL SELECT( 1 )
223 DOUBLE PRECISION DUM( 1 )
224* ..
225* .. External Subroutines ..
226 EXTERNAL dgebak, dgebal, dgehrd, dhseqr, dlabad, dlacpy,
228 $ xerbla
229* ..
230* .. External Functions ..
231 LOGICAL LSAME
232 INTEGER IDAMAX, ILAENV
233 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2
234 EXTERNAL lsame, idamax, ilaenv, dlamch, dlange, dlapy2,
235 $ dnrm2
236* ..
237* .. Intrinsic Functions ..
238 INTRINSIC max, sqrt
239* ..
240* .. Executable Statements ..
241*
242* Test the input arguments
243*
244 info = 0
245 lquery = ( lwork.EQ.-1 )
246 wantvl = lsame( jobvl, 'V' )
247 wantvr = lsame( jobvr, 'V' )
248 IF( ( .NOT.wantvl ) .AND. ( .NOT.lsame( jobvl, 'N' ) ) ) THEN
249 info = -1
250 ELSE IF( ( .NOT.wantvr ) .AND. ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
251 info = -2
252 ELSE IF( n.LT.0 ) THEN
253 info = -3
254 ELSE IF( lda.LT.max( 1, n ) ) THEN
255 info = -5
256 ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
257 info = -9
258 ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
259 info = -11
260 END IF
261*
262* Compute workspace
263* (Note: Comments in the code beginning "Workspace:" describe the
264* minimal amount of workspace needed at that point in the code,
265* as well as the preferred amount for good performance.
266* NB refers to the optimal block size for the immediately
267* following subroutine, as returned by ILAENV.
268* HSWORK refers to the workspace preferred by DHSEQR, as
269* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
270* the worst case.)
271*
272 IF( info.EQ.0 ) THEN
273 IF( n.EQ.0 ) THEN
274 minwrk = 1
275 maxwrk = 1
276 ELSE
277 maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
278 IF( wantvl ) THEN
279 minwrk = 4*n
280 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
281 $ 'DORGHR', ' ', n, 1, n, -1 ) )
282 CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl, ldvl,
283 $ work, -1, info )
284 hswork = int( work(1) )
285 maxwrk = max( maxwrk, n + 1, n + hswork )
286 CALL dtrevc3( 'L', 'B', SELECT, n, a, lda,
287 $ vl, ldvl, vr, ldvr, n, nout,
288 $ work, -1, ierr )
289 lwork_trevc = int( work(1) )
290 maxwrk = max( maxwrk, n + lwork_trevc )
291 maxwrk = max( maxwrk, 4*n )
292 ELSE IF( wantvr ) THEN
293 minwrk = 4*n
294 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
295 $ 'DORGHR', ' ', n, 1, n, -1 ) )
296 CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
297 $ work, -1, info )
298 hswork = int( work(1) )
299 maxwrk = max( maxwrk, n + 1, n + hswork )
300 CALL dtrevc3( 'R', 'B', SELECT, n, a, lda,
301 $ vl, ldvl, vr, ldvr, n, nout,
302 $ work, -1, ierr )
303 lwork_trevc = int( work(1) )
304 maxwrk = max( maxwrk, n + lwork_trevc )
305 maxwrk = max( maxwrk, 4*n )
306 ELSE
307 minwrk = 3*n
308 CALL dhseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
309 $ work, -1, info )
310 hswork = int( work(1) )
311 maxwrk = max( maxwrk, n + 1, n + hswork )
312 END IF
313 maxwrk = max( maxwrk, minwrk )
314 END IF
315 work( 1 ) = maxwrk
316*
317 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
318 info = -13
319 END IF
320 END IF
321*
322 IF( info.NE.0 ) THEN
323 CALL xerbla( 'DGEEV ', -info )
324 RETURN
325 ELSE IF( lquery ) THEN
326 RETURN
327 END IF
328*
329* Quick return if possible
330*
331 IF( n.EQ.0 )
332 $ RETURN
333*
334* Get machine constants
335*
336 eps = dlamch( 'P' )
337 smlnum = dlamch( 'S' )
338 bignum = one / smlnum
339 CALL dlabad( smlnum, bignum )
340 smlnum = sqrt( smlnum ) / eps
341 bignum = one / smlnum
342*
343* Scale A if max element outside range [SMLNUM,BIGNUM]
344*
345 anrm = dlange( 'M', n, n, a, lda, dum )
346 scalea = .false.
347 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
348 scalea = .true.
349 cscale = smlnum
350 ELSE IF( anrm.GT.bignum ) THEN
351 scalea = .true.
352 cscale = bignum
353 END IF
354 IF( scalea )
355 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
356*
357* Balance the matrix
358* (Workspace: need N)
359*
360 ibal = 1
361 CALL dgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
362*
363* Reduce to upper Hessenberg form
364* (Workspace: need 3*N, prefer 2*N+N*NB)
365*
366 itau = ibal + n
367 iwrk = itau + n
368 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
369 $ lwork-iwrk+1, ierr )
370*
371 IF( wantvl ) THEN
372*
373* Want left eigenvectors
374* Copy Householder vectors to VL
375*
376 side = 'L'
377 CALL dlacpy( 'L', n, n, a, lda, vl, ldvl )
378*
379* Generate orthogonal matrix in VL
380* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
381*
382 CALL dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
383 $ lwork-iwrk+1, ierr )
384*
385* Perform QR iteration, accumulating Schur vectors in VL
386* (Workspace: need N+1, prefer N+HSWORK (see comments) )
387*
388 iwrk = itau
389 CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
390 $ work( iwrk ), lwork-iwrk+1, info )
391*
392 IF( wantvr ) THEN
393*
394* Want left and right eigenvectors
395* Copy Schur vectors to VR
396*
397 side = 'B'
398 CALL dlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
399 END IF
400*
401 ELSE IF( wantvr ) THEN
402*
403* Want right eigenvectors
404* Copy Householder vectors to VR
405*
406 side = 'R'
407 CALL dlacpy( 'L', n, n, a, lda, vr, ldvr )
408*
409* Generate orthogonal matrix in VR
410* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
411*
412 CALL dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
413 $ lwork-iwrk+1, ierr )
414*
415* Perform QR iteration, accumulating Schur vectors in VR
416* (Workspace: need N+1, prefer N+HSWORK (see comments) )
417*
418 iwrk = itau
419 CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
420 $ work( iwrk ), lwork-iwrk+1, info )
421*
422 ELSE
423*
424* Compute eigenvalues only
425* (Workspace: need N+1, prefer N+HSWORK (see comments) )
426*
427 iwrk = itau
428 CALL dhseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
429 $ work( iwrk ), lwork-iwrk+1, info )
430 END IF
431*
432* If INFO .NE. 0 from DHSEQR, then quit
433*
434 IF( info.NE.0 )
435 $ GO TO 50
436*
437 IF( wantvl .OR. wantvr ) THEN
438*
439* Compute left and/or right eigenvectors
440* (Workspace: need 4*N, prefer N + N + 2*N*NB)
441*
442 CALL dtrevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
443 $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
444 END IF
445*
446 IF( wantvl ) THEN
447*
448* Undo balancing of left eigenvectors
449* (Workspace: need N)
450*
451 CALL dgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,
452 $ ierr )
453*
454* Normalize left eigenvectors and make largest component real
455*
456 DO 20 i = 1, n
457 IF( wi( i ).EQ.zero ) THEN
458 scl = one / dnrm2( n, vl( 1, i ), 1 )
459 CALL dscal( n, scl, vl( 1, i ), 1 )
460 ELSE IF( wi( i ).GT.zero ) THEN
461 scl = one / dlapy2( dnrm2( n, vl( 1, i ), 1 ),
462 $ dnrm2( n, vl( 1, i+1 ), 1 ) )
463 CALL dscal( n, scl, vl( 1, i ), 1 )
464 CALL dscal( n, scl, vl( 1, i+1 ), 1 )
465 DO 10 k = 1, n
466 work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
467 10 CONTINUE
468 k = idamax( n, work( iwrk ), 1 )
469 CALL dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
470 CALL drot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
471 vl( k, i+1 ) = zero
472 END IF
473 20 CONTINUE
474 END IF
475*
476 IF( wantvr ) THEN
477*
478* Undo balancing of right eigenvectors
479* (Workspace: need N)
480*
481 CALL dgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,
482 $ ierr )
483*
484* Normalize right eigenvectors and make largest component real
485*
486 DO 40 i = 1, n
487 IF( wi( i ).EQ.zero ) THEN
488 scl = one / dnrm2( n, vr( 1, i ), 1 )
489 CALL dscal( n, scl, vr( 1, i ), 1 )
490 ELSE IF( wi( i ).GT.zero ) THEN
491 scl = one / dlapy2( dnrm2( n, vr( 1, i ), 1 ),
492 $ dnrm2( n, vr( 1, i+1 ), 1 ) )
493 CALL dscal( n, scl, vr( 1, i ), 1 )
494 CALL dscal( n, scl, vr( 1, i+1 ), 1 )
495 DO 30 k = 1, n
496 work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
497 30 CONTINUE
498 k = idamax( n, work( iwrk ), 1 )
499 CALL dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
500 CALL drot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
501 vr( k, i+1 ) = zero
502 END IF
503 40 CONTINUE
504 END IF
505*
506* Undo scaling if necessary
507*
508 50 CONTINUE
509 IF( scalea ) THEN
510 CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
511 $ max( n-info, 1 ), ierr )
512 CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
513 $ max( n-info, 1 ), ierr )
514 IF( info.GT.0 ) THEN
515 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
516 $ ierr )
517 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
518 $ ierr )
519 END IF
520 END IF
521*
522 work( 1 ) = maxwrk
523 RETURN
524*
525* End of DGEEV
526*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:111
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
double precision function dlapy2(X, Y)
DLAPY2 returns sqrt(x2+y2).
Definition: dlapy2.f:63
integer function idamax(N, DX, INCX)
IDAMAX
Definition: idamax.f:71
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
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
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:167
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
Definition: dgebal.f:160
subroutine dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
Definition: dgebak.f:130
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:316
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:126
subroutine dtrevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
DTREVC3
Definition: dtrevc3.f:237
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition: dnrm2.f90:89
Here is the call graph for this function:
Here is the caller graph for this function: