LAPACK  3.8.0
LAPACK: Linear Algebra PACKage

◆ 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.
Date
June 2016

Definition at line 193 of file dgeev.f.

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