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

◆ dggev3()

subroutine dggev3 ( character  jobvl,
character  jobvr,
integer  n,
double precision, dimension( lda, * )  a,
integer  lda,
double precision, dimension( ldb, * )  b,
integer  ldb,
double precision, dimension( * )  alphar,
double precision, dimension( * )  alphai,
double precision, dimension( * )  beta,
double precision, dimension( ldvl, * )  vl,
integer  ldvl,
double precision, dimension( ldvr, * )  vr,
integer  ldvr,
double precision, dimension( * )  work,
integer  lwork,
integer  info 
)

DGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)

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

Purpose:
 DGGEV3 computes for a pair of N-by-N real nonsymmetric matrices (A,B)
 the generalized eigenvalues, and optionally, the left and/or right
 generalized eigenvectors.

 A generalized eigenvalue for a pair of matrices (A,B) is a scalar
 lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
 singular. It is usually represented as the pair (alpha,beta), as
 there is a reasonable interpretation for beta=0, and even for both
 being zero.

 The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
 of (A,B) satisfies

                  A * v(j) = lambda(j) * B * v(j).

 The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
 of (A,B) satisfies

                  u(j)**H * A  = lambda(j) * u(j)**H * B .

 where u(j)**H is the conjugate-transpose of u(j).
Parameters
[in]JOBVL
          JOBVL is CHARACTER*1
          = 'N':  do not compute the left generalized eigenvectors;
          = 'V':  compute the left generalized eigenvectors.
[in]JOBVR
          JOBVR is CHARACTER*1
          = 'N':  do not compute the right generalized eigenvectors;
          = 'V':  compute the right generalized eigenvectors.
[in]N
          N is INTEGER
          The order of the matrices A, B, VL, and VR.  N >= 0.
[in,out]A
          A is DOUBLE PRECISION array, dimension (LDA, N)
          On entry, the matrix A in the pair (A,B).
          On exit, A has been overwritten.
[in]LDA
          LDA is INTEGER
          The leading dimension of A.  LDA >= max(1,N).
[in,out]B
          B is DOUBLE PRECISION array, dimension (LDB, N)
          On entry, the matrix B in the pair (A,B).
          On exit, B has been overwritten.
[in]LDB
          LDB is INTEGER
          The leading dimension of B.  LDB >= max(1,N).
[out]ALPHAR
          ALPHAR is DOUBLE PRECISION array, dimension (N)
[out]ALPHAI
          ALPHAI is DOUBLE PRECISION array, dimension (N)
[out]BETA
          BETA is DOUBLE PRECISION array, dimension (N)
          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
          the j-th eigenvalue is real; if positive, then the j-th and
          (j+1)-st eigenvalues are a complex conjugate pair, with
          ALPHAI(j+1) negative.

          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
          may easily over- or underflow, and BETA(j) may even be zero.
          Thus, the user should avoid naively computing the ratio
          alpha/beta.  However, ALPHAR and ALPHAI will be always less
          than and usually comparable with norm(A) in magnitude, and
          BETA always less than and usually comparable with norm(B).
[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 the j-th eigenvalue is real, then
          u(j) = VL(:,j), the j-th column of VL. If the j-th and
          (j+1)-th 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).
          Each eigenvector is scaled so the largest component has
          abs(real part)+abs(imag. part)=1.
          Not referenced if JOBVL = 'N'.
[in]LDVL
          LDVL is INTEGER
          The leading dimension of the matrix VL. LDVL >= 1, and
          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 the j-th eigenvalue is real, then
          v(j) = VR(:,j), the j-th column of VR. If the j-th and
          (j+1)-th 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).
          Each eigenvector is scaled so the largest component has
          abs(real part)+abs(imag. part)=1.
          Not referenced if JOBVR = 'N'.
[in]LDVR
          LDVR is INTEGER
          The leading dimension of the matrix VR. LDVR >= 1, and
          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

          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.
          = 1,...,N:
                The QZ iteration failed.  No eigenvectors have been
                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
                should be correct for j=INFO+1,...,N.
          > N:  =N+1: other than QZ iteration failed in DLAQZ0.
                =N+2: error return from DTGEVC.
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 223 of file dggev3.f.

226*
227* -- LAPACK driver routine --
228* -- LAPACK is a software package provided by Univ. of Tennessee, --
229* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231* .. Scalar Arguments ..
232 CHARACTER JOBVL, JOBVR
233 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
234* ..
235* .. Array Arguments ..
236 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
237 $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
238 $ VR( LDVR, * ), WORK( * )
239* ..
240*
241* =====================================================================
242*
243* .. Parameters ..
244 DOUBLE PRECISION ZERO, ONE
245 parameter( zero = 0.0d+0, one = 1.0d+0 )
246* ..
247* .. Local Scalars ..
248 LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
249 CHARACTER CHTEMP
250 INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
251 $ IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, LWKOPT
252 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
253 $ SMLNUM, TEMP
254* ..
255* .. Local Arrays ..
256 LOGICAL LDUMMA( 1 )
257* ..
258* .. External Subroutines ..
259 EXTERNAL dgeqrf, dggbak, dggbal, dgghd3, dlaqz0, dlacpy,
261* ..
262* .. External Functions ..
263 LOGICAL LSAME
264 DOUBLE PRECISION DLAMCH, DLANGE
265 EXTERNAL lsame, dlamch, dlange
266* ..
267* .. Intrinsic Functions ..
268 INTRINSIC abs, max, sqrt
269* ..
270* .. Executable Statements ..
271*
272* Decode the input arguments
273*
274 IF( lsame( jobvl, 'N' ) ) THEN
275 ijobvl = 1
276 ilvl = .false.
277 ELSE IF( lsame( jobvl, 'V' ) ) THEN
278 ijobvl = 2
279 ilvl = .true.
280 ELSE
281 ijobvl = -1
282 ilvl = .false.
283 END IF
284*
285 IF( lsame( jobvr, 'N' ) ) THEN
286 ijobvr = 1
287 ilvr = .false.
288 ELSE IF( lsame( jobvr, 'V' ) ) THEN
289 ijobvr = 2
290 ilvr = .true.
291 ELSE
292 ijobvr = -1
293 ilvr = .false.
294 END IF
295 ilv = ilvl .OR. ilvr
296*
297* Test the input arguments
298*
299 info = 0
300 lquery = ( lwork.EQ.-1 )
301 IF( ijobvl.LE.0 ) THEN
302 info = -1
303 ELSE IF( ijobvr.LE.0 ) THEN
304 info = -2
305 ELSE IF( n.LT.0 ) THEN
306 info = -3
307 ELSE IF( lda.LT.max( 1, n ) ) THEN
308 info = -5
309 ELSE IF( ldb.LT.max( 1, n ) ) THEN
310 info = -7
311 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
312 info = -12
313 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
314 info = -14
315 ELSE IF( lwork.LT.max( 1, 8*n ) .AND. .NOT.lquery ) THEN
316 info = -16
317 END IF
318*
319* Compute workspace
320*
321 IF( info.EQ.0 ) THEN
322 CALL dgeqrf( n, n, b, ldb, work, work, -1, ierr )
323 lwkopt = max(1, 8*n, 3*n+int( work( 1 ) ) )
324 CALL dormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work, -1,
325 $ ierr )
326 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
327 IF( ilvl ) THEN
328 CALL dorgqr( n, n, n, vl, ldvl, work, work, -1, ierr )
329 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
330 END IF
331 IF( ilv ) THEN
332 CALL dgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
333 $ ldvl, vr, ldvr, work, -1, ierr )
334 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
335 CALL dlaqz0( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
336 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
337 $ work, -1, 0, ierr )
338 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
339 ELSE
340 CALL dgghd3( 'N', 'N', n, 1, n, a, lda, b, ldb, vl, ldvl,
341 $ vr, ldvr, work, -1, ierr )
342 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
343 CALL dlaqz0( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
344 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
345 $ work, -1, 0, ierr )
346 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
347 END IF
348
349 work( 1 ) = lwkopt
350 END IF
351*
352 IF( info.NE.0 ) THEN
353 CALL xerbla( 'DGGEV3 ', -info )
354 RETURN
355 ELSE IF( lquery ) THEN
356 RETURN
357 END IF
358*
359* Quick return if possible
360*
361 IF( n.EQ.0 )
362 $ RETURN
363*
364* Get machine constants
365*
366 eps = dlamch( 'P' )
367 smlnum = dlamch( 'S' )
368 bignum = one / smlnum
369 smlnum = sqrt( smlnum ) / eps
370 bignum = one / smlnum
371*
372* Scale A if max element outside range [SMLNUM,BIGNUM]
373*
374 anrm = dlange( 'M', n, n, a, lda, work )
375 ilascl = .false.
376 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
377 anrmto = smlnum
378 ilascl = .true.
379 ELSE IF( anrm.GT.bignum ) THEN
380 anrmto = bignum
381 ilascl = .true.
382 END IF
383 IF( ilascl )
384 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
385*
386* Scale B if max element outside range [SMLNUM,BIGNUM]
387*
388 bnrm = dlange( 'M', n, n, b, ldb, work )
389 ilbscl = .false.
390 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
391 bnrmto = smlnum
392 ilbscl = .true.
393 ELSE IF( bnrm.GT.bignum ) THEN
394 bnrmto = bignum
395 ilbscl = .true.
396 END IF
397 IF( ilbscl )
398 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
399*
400* Permute the matrices A, B to isolate eigenvalues if possible
401*
402 ileft = 1
403 iright = n + 1
404 iwrk = iright + n
405 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
406 $ work( iright ), work( iwrk ), ierr )
407*
408* Reduce B to triangular form (QR decomposition of B)
409*
410 irows = ihi + 1 - ilo
411 IF( ilv ) THEN
412 icols = n + 1 - ilo
413 ELSE
414 icols = irows
415 END IF
416 itau = iwrk
417 iwrk = itau + irows
418 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
419 $ work( iwrk ), lwork+1-iwrk, ierr )
420*
421* Apply the orthogonal transformation to matrix A
422*
423 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
424 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
425 $ lwork+1-iwrk, ierr )
426*
427* Initialize VL
428*
429 IF( ilvl ) THEN
430 CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
431 IF( irows.GT.1 ) THEN
432 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
433 $ vl( ilo+1, ilo ), ldvl )
434 END IF
435 CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
436 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
437 END IF
438*
439* Initialize VR
440*
441 IF( ilvr )
442 $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
443*
444* Reduce to generalized Hessenberg form
445*
446 IF( ilv ) THEN
447*
448* Eigenvectors requested -- work on whole matrix.
449*
450 CALL dgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
451 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
452 ELSE
453 CALL dgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
454 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
455 $ work( iwrk ), lwork+1-iwrk, ierr )
456 END IF
457*
458* Perform QZ algorithm (Compute eigenvalues, and optionally, the
459* Schur forms and Schur vectors)
460*
461 iwrk = itau
462 IF( ilv ) THEN
463 chtemp = 'S'
464 ELSE
465 chtemp = 'E'
466 END IF
467 CALL dlaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
468 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
469 $ work( iwrk ), lwork+1-iwrk, 0, ierr )
470 IF( ierr.NE.0 ) THEN
471 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
472 info = ierr
473 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
474 info = ierr - n
475 ELSE
476 info = n + 1
477 END IF
478 GO TO 110
479 END IF
480*
481* Compute Eigenvectors
482*
483 IF( ilv ) THEN
484 IF( ilvl ) THEN
485 IF( ilvr ) THEN
486 chtemp = 'B'
487 ELSE
488 chtemp = 'L'
489 END IF
490 ELSE
491 chtemp = 'R'
492 END IF
493 CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
494 $ vr, ldvr, n, in, work( iwrk ), ierr )
495 IF( ierr.NE.0 ) THEN
496 info = n + 2
497 GO TO 110
498 END IF
499*
500* Undo balancing on VL and VR and normalization
501*
502 IF( ilvl ) THEN
503 CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
504 $ work( iright ), n, vl, ldvl, ierr )
505 DO 50 jc = 1, n
506 IF( alphai( jc ).LT.zero )
507 $ GO TO 50
508 temp = zero
509 IF( alphai( jc ).EQ.zero ) THEN
510 DO 10 jr = 1, n
511 temp = max( temp, abs( vl( jr, jc ) ) )
512 10 CONTINUE
513 ELSE
514 DO 20 jr = 1, n
515 temp = max( temp, abs( vl( jr, jc ) )+
516 $ abs( vl( jr, jc+1 ) ) )
517 20 CONTINUE
518 END IF
519 IF( temp.LT.smlnum )
520 $ GO TO 50
521 temp = one / temp
522 IF( alphai( jc ).EQ.zero ) THEN
523 DO 30 jr = 1, n
524 vl( jr, jc ) = vl( jr, jc )*temp
525 30 CONTINUE
526 ELSE
527 DO 40 jr = 1, n
528 vl( jr, jc ) = vl( jr, jc )*temp
529 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
530 40 CONTINUE
531 END IF
532 50 CONTINUE
533 END IF
534 IF( ilvr ) THEN
535 CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
536 $ work( iright ), n, vr, ldvr, ierr )
537 DO 100 jc = 1, n
538 IF( alphai( jc ).LT.zero )
539 $ GO TO 100
540 temp = zero
541 IF( alphai( jc ).EQ.zero ) THEN
542 DO 60 jr = 1, n
543 temp = max( temp, abs( vr( jr, jc ) ) )
544 60 CONTINUE
545 ELSE
546 DO 70 jr = 1, n
547 temp = max( temp, abs( vr( jr, jc ) )+
548 $ abs( vr( jr, jc+1 ) ) )
549 70 CONTINUE
550 END IF
551 IF( temp.LT.smlnum )
552 $ GO TO 100
553 temp = one / temp
554 IF( alphai( jc ).EQ.zero ) THEN
555 DO 80 jr = 1, n
556 vr( jr, jc ) = vr( jr, jc )*temp
557 80 CONTINUE
558 ELSE
559 DO 90 jr = 1, n
560 vr( jr, jc ) = vr( jr, jc )*temp
561 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
562 90 CONTINUE
563 END IF
564 100 CONTINUE
565 END IF
566*
567* End of eigenvector calculation
568*
569 END IF
570*
571* Undo scaling if necessary
572*
573 110 CONTINUE
574*
575 IF( ilascl ) THEN
576 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
577 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
578 END IF
579*
580 IF( ilbscl ) THEN
581 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
582 END IF
583*
584 work( 1 ) = lwkopt
585 RETURN
586*
587* End of DGGEV3
588*
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dgeqrf(m, n, a, lda, tau, work, lwork, info)
DGEQRF
Definition dgeqrf.f:146
subroutine dggbak(job, side, n, ilo, ihi, lscale, rscale, m, v, ldv, info)
DGGBAK
Definition dggbak.f:147
subroutine dggbal(job, n, a, lda, b, ldb, ilo, ihi, lscale, rscale, work, info)
DGGBAL
Definition dggbal.f:177
subroutine dgghd3(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, work, lwork, info)
DGGHD3
Definition dgghd3.f:230
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
recursive subroutine dlaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, work, lwork, rec, info)
DLAQZ0
Definition dlaqz0.f:306
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 dlaset(uplo, m, n, alpha, beta, a, lda)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition dlaset.f:110
logical function lsame(ca, cb)
LSAME
Definition lsame.f:48
subroutine dtgevc(side, howmny, select, n, s, lds, p, ldp, vl, ldvl, vr, ldvr, mm, m, work, info)
DTGEVC
Definition dtgevc.f:295
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:128
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:167
Here is the call graph for this function:
Here is the caller graph for this function: