LAPACK 3.12.1
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 188 of file dgeev.f.

191 implicit none
192*
193* -- LAPACK driver routine --
194* -- LAPACK is a software package provided by Univ. of Tennessee, --
195* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
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 $ LWORK_TREVC, 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, dlacpy,
226 $ dlartg,
227 $ dlascl, dorghr, drot, dscal, dtrevc3, 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,
234 $ 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.
251 $ ( .NOT.lsame( jobvr, 'N' ) ) ) THEN
252 info = -2
253 ELSE IF( n.LT.0 ) THEN
254 info = -3
255 ELSE IF( lda.LT.max( 1, n ) ) THEN
256 info = -5
257 ELSE IF( ldvl.LT.1 .OR. ( wantvl .AND. ldvl.LT.n ) ) THEN
258 info = -9
259 ELSE IF( ldvr.LT.1 .OR. ( wantvr .AND. ldvr.LT.n ) ) THEN
260 info = -11
261 END IF
262*
263* Compute workspace
264* (Note: Comments in the code beginning "Workspace:" describe the
265* minimal amount of workspace needed at that point in the code,
266* as well as the preferred amount for good performance.
267* NB refers to the optimal block size for the immediately
268* following subroutine, as returned by ILAENV.
269* HSWORK refers to the workspace preferred by DHSEQR, as
270* calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
271* the worst case.)
272*
273 IF( info.EQ.0 ) THEN
274 IF( n.EQ.0 ) THEN
275 minwrk = 1
276 maxwrk = 1
277 ELSE
278 maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
279 IF( wantvl ) THEN
280 minwrk = 4*n
281 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
282 $ 'DORGHR', ' ', n, 1, n, -1 ) )
283 CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vl,
284 $ 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,
299 $ ldvr,
300 $ work, -1, info )
301 hswork = int( work(1) )
302 maxwrk = max( maxwrk, n + 1, n + hswork )
303 CALL dtrevc3( 'R', 'B', SELECT, n, a, lda,
304 $ vl, ldvl, vr, ldvr, n, nout,
305 $ work, -1, ierr )
306 lwork_trevc = int( work(1) )
307 maxwrk = max( maxwrk, n + lwork_trevc )
308 maxwrk = max( maxwrk, 4*n )
309 ELSE
310 minwrk = 3*n
311 CALL dhseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr,
312 $ ldvr,
313 $ work, -1, info )
314 hswork = int( work(1) )
315 maxwrk = max( maxwrk, n + 1, n + hswork )
316 END IF
317 maxwrk = max( maxwrk, minwrk )
318 END IF
319 work( 1 ) = maxwrk
320*
321 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
322 info = -13
323 END IF
324 END IF
325*
326 IF( info.NE.0 ) THEN
327 CALL xerbla( 'DGEEV ', -info )
328 RETURN
329 ELSE IF( lquery ) THEN
330 RETURN
331 END IF
332*
333* Quick return if possible
334*
335 IF( n.EQ.0 )
336 $ RETURN
337*
338* Get machine constants
339*
340 eps = dlamch( 'P' )
341 smlnum = dlamch( 'S' )
342 bignum = one / smlnum
343 smlnum = sqrt( smlnum ) / eps
344 bignum = one / smlnum
345*
346* Scale A if max element outside range [SMLNUM,BIGNUM]
347*
348 anrm = dlange( 'M', n, n, a, lda, dum )
349 scalea = .false.
350 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
351 scalea = .true.
352 cscale = smlnum
353 ELSE IF( anrm.GT.bignum ) THEN
354 scalea = .true.
355 cscale = bignum
356 END IF
357 IF( scalea )
358 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
359*
360* Balance the matrix
361* (Workspace: need N)
362*
363 ibal = 1
364 CALL dgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
365*
366* Reduce to upper Hessenberg form
367* (Workspace: need 3*N, prefer 2*N+N*NB)
368*
369 itau = ibal + n
370 iwrk = itau + n
371 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
372 $ lwork-iwrk+1, ierr )
373*
374 IF( wantvl ) THEN
375*
376* Want left eigenvectors
377* Copy Householder vectors to VL
378*
379 side = 'L'
380 CALL dlacpy( 'L', n, n, a, lda, vl, ldvl )
381*
382* Generate orthogonal matrix in VL
383* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
384*
385 CALL dorghr( n, ilo, ihi, vl, ldvl, work( itau ),
386 $ work( iwrk ),
387 $ lwork-iwrk+1, ierr )
388*
389* Perform QR iteration, accumulating Schur vectors in VL
390* (Workspace: need N+1, prefer N+HSWORK (see comments) )
391*
392 iwrk = itau
393 CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl,
394 $ ldvl,
395 $ work( iwrk ), lwork-iwrk+1, info )
396*
397 IF( wantvr ) THEN
398*
399* Want left and right eigenvectors
400* Copy Schur vectors to VR
401*
402 side = 'B'
403 CALL dlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
404 END IF
405*
406 ELSE IF( wantvr ) THEN
407*
408* Want right eigenvectors
409* Copy Householder vectors to VR
410*
411 side = 'R'
412 CALL dlacpy( 'L', n, n, a, lda, vr, ldvr )
413*
414* Generate orthogonal matrix in VR
415* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
416*
417 CALL dorghr( n, ilo, ihi, vr, ldvr, work( itau ),
418 $ work( iwrk ),
419 $ lwork-iwrk+1, ierr )
420*
421* Perform QR iteration, accumulating Schur vectors in VR
422* (Workspace: need N+1, prefer N+HSWORK (see comments) )
423*
424 iwrk = itau
425 CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr,
426 $ ldvr,
427 $ work( iwrk ), lwork-iwrk+1, info )
428*
429 ELSE
430*
431* Compute eigenvalues only
432* (Workspace: need N+1, prefer N+HSWORK (see comments) )
433*
434 iwrk = itau
435 CALL dhseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr,
436 $ ldvr,
437 $ work( iwrk ), lwork-iwrk+1, info )
438 END IF
439*
440* If INFO .NE. 0 from DHSEQR, then quit
441*
442 IF( info.NE.0 )
443 $ GO TO 50
444*
445 IF( wantvl .OR. wantvr ) THEN
446*
447* Compute left and/or right eigenvectors
448* (Workspace: need 4*N, prefer N + N + 2*N*NB)
449*
450 CALL dtrevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr,
451 $ ldvr,
452 $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
453 END IF
454*
455 IF( wantvl ) THEN
456*
457* Undo balancing of left eigenvectors
458* (Workspace: need N)
459*
460 CALL dgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl,
461 $ ldvl,
462 $ ierr )
463*
464* Normalize left eigenvectors and make largest component real
465*
466 DO 20 i = 1, n
467 IF( wi( i ).EQ.zero ) THEN
468 scl = one / dnrm2( n, vl( 1, i ), 1 )
469 CALL dscal( n, scl, vl( 1, i ), 1 )
470 ELSE IF( wi( i ).GT.zero ) THEN
471 scl = one / dlapy2( dnrm2( n, vl( 1, i ), 1 ),
472 $ dnrm2( n, vl( 1, i+1 ), 1 ) )
473 CALL dscal( n, scl, vl( 1, i ), 1 )
474 CALL dscal( n, scl, vl( 1, i+1 ), 1 )
475 DO 10 k = 1, n
476 work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
477 10 CONTINUE
478 k = idamax( n, work( iwrk ), 1 )
479 CALL dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
480 CALL drot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
481 vl( k, i+1 ) = zero
482 END IF
483 20 CONTINUE
484 END IF
485*
486 IF( wantvr ) THEN
487*
488* Undo balancing of right eigenvectors
489* (Workspace: need N)
490*
491 CALL dgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr,
492 $ ldvr,
493 $ ierr )
494*
495* Normalize right eigenvectors and make largest component real
496*
497 DO 40 i = 1, n
498 IF( wi( i ).EQ.zero ) THEN
499 scl = one / dnrm2( n, vr( 1, i ), 1 )
500 CALL dscal( n, scl, vr( 1, i ), 1 )
501 ELSE IF( wi( i ).GT.zero ) THEN
502 scl = one / dlapy2( dnrm2( n, vr( 1, i ), 1 ),
503 $ dnrm2( n, vr( 1, i+1 ), 1 ) )
504 CALL dscal( n, scl, vr( 1, i ), 1 )
505 CALL dscal( n, scl, vr( 1, i+1 ), 1 )
506 DO 30 k = 1, n
507 work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
508 30 CONTINUE
509 k = idamax( n, work( iwrk ), 1 )
510 CALL dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
511 CALL drot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
512 vr( k, i+1 ) = zero
513 END IF
514 40 CONTINUE
515 END IF
516*
517* Undo scaling if necessary
518*
519 50 CONTINUE
520 IF( scalea ) THEN
521 CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1,
522 $ wr( info+1 ),
523 $ max( n-info, 1 ), ierr )
524 CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1,
525 $ wi( info+1 ),
526 $ max( n-info, 1 ), ierr )
527 IF( info.GT.0 ) THEN
528 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
529 $ ierr )
530 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
531 $ ierr )
532 END IF
533 END IF
534*
535 work( 1 ) = maxwrk
536 RETURN
537*
538* End of DGEEV
539*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dgebak
subroutine dgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
DGEBAK
Definition dgebak.f:128
dgebal
subroutine dgebal(job, n, a, lda, ilo, ihi, scale, info)
DGEBAL
Definition dgebal.f:161
dgehrd
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:166
dhseqr
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
Definition dhseqr.f:314
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
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
dlapy2
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:61
dlartg
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
dlascl
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:142
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
dnrm2
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
drot
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92
dscal
subroutine dscal(n, da, dx, incx)
DSCAL
Definition dscal.f:79
dtrevc3
subroutine dtrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
DTREVC3
Definition dtrevc3.f:235
dorghr
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR
Definition dorghr.f:125
Here is the call graph for this function:
Here is the caller graph for this function: