LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
dgeev.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dgeev (JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
  DGEEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices More...
 

Function/Subroutine Documentation

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
September 2012

Definition at line 191 of file dgeev.f.

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

Here is the call graph for this function:

Here is the caller graph for this function: