LAPACK 3.12.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, dlacpy, dlartg,
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 = int( work(1) )
284 maxwrk = max( maxwrk, n + 1, n + hswork )
285 CALL dtrevc3( 'L', 'B', SELECT, n, a, lda,
286 $ vl, ldvl, vr, ldvr, n, nout,
287 $ work, -1, ierr )
288 lwork_trevc = int( work(1) )
289 maxwrk = max( maxwrk, n + lwork_trevc )
290 maxwrk = max( maxwrk, 4*n )
291 ELSE IF( wantvr ) THEN
292 minwrk = 4*n
293 maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
294 $ 'DORGHR', ' ', n, 1, n, -1 ) )
295 CALL dhseqr( 'S', 'V', n, 1, n, a, lda, wr, wi, vr, ldvr,
296 $ work, -1, info )
297 hswork = int( work(1) )
298 maxwrk = max( maxwrk, n + 1, n + hswork )
299 CALL dtrevc3( 'R', 'B', SELECT, n, a, lda,
300 $ vl, ldvl, vr, ldvr, n, nout,
301 $ work, -1, ierr )
302 lwork_trevc = int( work(1) )
303 maxwrk = max( maxwrk, n + lwork_trevc )
304 maxwrk = max( maxwrk, 4*n )
305 ELSE
306 minwrk = 3*n
307 CALL dhseqr( 'E', 'N', n, 1, n, a, lda, wr, wi, vr, ldvr,
308 $ work, -1, info )
309 hswork = int( work(1) )
310 maxwrk = max( maxwrk, n + 1, n + hswork )
311 END IF
312 maxwrk = max( maxwrk, minwrk )
313 END IF
314 work( 1 ) = maxwrk
315*
316 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
317 info = -13
318 END IF
319 END IF
320*
321 IF( info.NE.0 ) THEN
322 CALL xerbla( 'DGEEV ', -info )
323 RETURN
324 ELSE IF( lquery ) THEN
325 RETURN
326 END IF
327*
328* Quick return if possible
329*
330 IF( n.EQ.0 )
331 $ RETURN
332*
333* Get machine constants
334*
335 eps = dlamch( 'P' )
336 smlnum = dlamch( 'S' )
337 bignum = one / smlnum
338 smlnum = sqrt( smlnum ) / eps
339 bignum = one / smlnum
340*
341* Scale A if max element outside range [SMLNUM,BIGNUM]
342*
343 anrm = dlange( 'M', n, n, a, lda, dum )
344 scalea = .false.
345 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
346 scalea = .true.
347 cscale = smlnum
348 ELSE IF( anrm.GT.bignum ) THEN
349 scalea = .true.
350 cscale = bignum
351 END IF
352 IF( scalea )
353 $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
354*
355* Balance the matrix
356* (Workspace: need N)
357*
358 ibal = 1
359 CALL dgebal( 'B', n, a, lda, ilo, ihi, work( ibal ), ierr )
360*
361* Reduce to upper Hessenberg form
362* (Workspace: need 3*N, prefer 2*N+N*NB)
363*
364 itau = ibal + n
365 iwrk = itau + n
366 CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
367 $ lwork-iwrk+1, ierr )
368*
369 IF( wantvl ) THEN
370*
371* Want left eigenvectors
372* Copy Householder vectors to VL
373*
374 side = 'L'
375 CALL dlacpy( 'L', n, n, a, lda, vl, ldvl )
376*
377* Generate orthogonal matrix in VL
378* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
379*
380 CALL dorghr( n, ilo, ihi, vl, ldvl, work( itau ), work( iwrk ),
381 $ lwork-iwrk+1, ierr )
382*
383* Perform QR iteration, accumulating Schur vectors in VL
384* (Workspace: need N+1, prefer N+HSWORK (see comments) )
385*
386 iwrk = itau
387 CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vl, ldvl,
388 $ work( iwrk ), lwork-iwrk+1, info )
389*
390 IF( wantvr ) THEN
391*
392* Want left and right eigenvectors
393* Copy Schur vectors to VR
394*
395 side = 'B'
396 CALL dlacpy( 'F', n, n, vl, ldvl, vr, ldvr )
397 END IF
398*
399 ELSE IF( wantvr ) THEN
400*
401* Want right eigenvectors
402* Copy Householder vectors to VR
403*
404 side = 'R'
405 CALL dlacpy( 'L', n, n, a, lda, vr, ldvr )
406*
407* Generate orthogonal matrix in VR
408* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
409*
410 CALL dorghr( n, ilo, ihi, vr, ldvr, work( itau ), work( iwrk ),
411 $ lwork-iwrk+1, ierr )
412*
413* Perform QR iteration, accumulating Schur vectors in VR
414* (Workspace: need N+1, prefer N+HSWORK (see comments) )
415*
416 iwrk = itau
417 CALL dhseqr( 'S', 'V', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
418 $ work( iwrk ), lwork-iwrk+1, info )
419*
420 ELSE
421*
422* Compute eigenvalues only
423* (Workspace: need N+1, prefer N+HSWORK (see comments) )
424*
425 iwrk = itau
426 CALL dhseqr( 'E', 'N', n, ilo, ihi, a, lda, wr, wi, vr, ldvr,
427 $ work( iwrk ), lwork-iwrk+1, info )
428 END IF
429*
430* If INFO .NE. 0 from DHSEQR, then quit
431*
432 IF( info.NE.0 )
433 $ GO TO 50
434*
435 IF( wantvl .OR. wantvr ) THEN
436*
437* Compute left and/or right eigenvectors
438* (Workspace: need 4*N, prefer N + N + 2*N*NB)
439*
440 CALL dtrevc3( side, 'B', SELECT, n, a, lda, vl, ldvl, vr, ldvr,
441 $ n, nout, work( iwrk ), lwork-iwrk+1, ierr )
442 END IF
443*
444 IF( wantvl ) THEN
445*
446* Undo balancing of left eigenvectors
447* (Workspace: need N)
448*
449 CALL dgebak( 'B', 'L', n, ilo, ihi, work( ibal ), n, vl, ldvl,
450 $ ierr )
451*
452* Normalize left eigenvectors and make largest component real
453*
454 DO 20 i = 1, n
455 IF( wi( i ).EQ.zero ) THEN
456 scl = one / dnrm2( n, vl( 1, i ), 1 )
457 CALL dscal( n, scl, vl( 1, i ), 1 )
458 ELSE IF( wi( i ).GT.zero ) THEN
459 scl = one / dlapy2( dnrm2( n, vl( 1, i ), 1 ),
460 $ dnrm2( n, vl( 1, i+1 ), 1 ) )
461 CALL dscal( n, scl, vl( 1, i ), 1 )
462 CALL dscal( n, scl, vl( 1, i+1 ), 1 )
463 DO 10 k = 1, n
464 work( iwrk+k-1 ) = vl( k, i )**2 + vl( k, i+1 )**2
465 10 CONTINUE
466 k = idamax( n, work( iwrk ), 1 )
467 CALL dlartg( vl( k, i ), vl( k, i+1 ), cs, sn, r )
468 CALL drot( n, vl( 1, i ), 1, vl( 1, i+1 ), 1, cs, sn )
469 vl( k, i+1 ) = zero
470 END IF
471 20 CONTINUE
472 END IF
473*
474 IF( wantvr ) THEN
475*
476* Undo balancing of right eigenvectors
477* (Workspace: need N)
478*
479 CALL dgebak( 'B', 'R', n, ilo, ihi, work( ibal ), n, vr, ldvr,
480 $ ierr )
481*
482* Normalize right eigenvectors and make largest component real
483*
484 DO 40 i = 1, n
485 IF( wi( i ).EQ.zero ) THEN
486 scl = one / dnrm2( n, vr( 1, i ), 1 )
487 CALL dscal( n, scl, vr( 1, i ), 1 )
488 ELSE IF( wi( i ).GT.zero ) THEN
489 scl = one / dlapy2( dnrm2( n, vr( 1, i ), 1 ),
490 $ dnrm2( n, vr( 1, i+1 ), 1 ) )
491 CALL dscal( n, scl, vr( 1, i ), 1 )
492 CALL dscal( n, scl, vr( 1, i+1 ), 1 )
493 DO 30 k = 1, n
494 work( iwrk+k-1 ) = vr( k, i )**2 + vr( k, i+1 )**2
495 30 CONTINUE
496 k = idamax( n, work( iwrk ), 1 )
497 CALL dlartg( vr( k, i ), vr( k, i+1 ), cs, sn, r )
498 CALL drot( n, vr( 1, i ), 1, vr( 1, i+1 ), 1, cs, sn )
499 vr( k, i+1 ) = zero
500 END IF
501 40 CONTINUE
502 END IF
503*
504* Undo scaling if necessary
505*
506 50 CONTINUE
507 IF( scalea ) THEN
508 CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wr( info+1 ),
509 $ max( n-info, 1 ), ierr )
510 CALL dlascl( 'G', 0, 0, cscale, anrm, n-info, 1, wi( info+1 ),
511 $ max( n-info, 1 ), ierr )
512 IF( info.GT.0 ) THEN
513 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wr, n,
514 $ ierr )
515 CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
516 $ ierr )
517 END IF
518 END IF
519*
520 work( 1 ) = maxwrk
521 RETURN
522*
523* End of DGEEV
524*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgebak(job, side, n, ilo, ihi, scale, m, v, ldv, info)
DGEBAK
Definition dgebak.f:130
subroutine dgebal(job, n, a, lda, ilo, ihi, scale, info)
DGEBAL
Definition dgebal.f:163
subroutine dgehrd(n, ilo, ihi, a, lda, tau, work, lwork, info)
DGEHRD
Definition dgehrd.f:167
subroutine dhseqr(job, compz, n, ilo, ihi, h, ldh, wr, wi, z, ldz, work, lwork, info)
DHSEQR
Definition dhseqr.f:316
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:162
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
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
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
double precision function dlapy2(x, y)
DLAPY2 returns sqrt(x2+y2).
Definition dlapy2.f:63
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
real(wp) function dnrm2(n, x, incx)
DNRM2
Definition dnrm2.f90:89
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
subroutine dtrevc3(side, howmny, select, n, t, ldt, vl, ldvl, vr, ldvr, mm, m, work, lwork, info)
DTREVC3
Definition dtrevc3.f:237
subroutine dorghr(n, ilo, ihi, a, lda, tau, work, lwork, info)
DORGHR
Definition dorghr.f:126
Here is the call graph for this function:
Here is the caller graph for this function: