LAPACK 3.11.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, dlabad,
261 $ xerbla
262* ..
263* .. External Functions ..
264 LOGICAL LSAME
265 DOUBLE PRECISION DLAMCH, DLANGE
266 EXTERNAL lsame, dlamch, dlange
267* ..
268* .. Intrinsic Functions ..
269 INTRINSIC abs, max, sqrt
270* ..
271* .. Executable Statements ..
272*
273* Decode the input arguments
274*
275 IF( lsame( jobvl, 'N' ) ) THEN
276 ijobvl = 1
277 ilvl = .false.
278 ELSE IF( lsame( jobvl, 'V' ) ) THEN
279 ijobvl = 2
280 ilvl = .true.
281 ELSE
282 ijobvl = -1
283 ilvl = .false.
284 END IF
285*
286 IF( lsame( jobvr, 'N' ) ) THEN
287 ijobvr = 1
288 ilvr = .false.
289 ELSE IF( lsame( jobvr, 'V' ) ) THEN
290 ijobvr = 2
291 ilvr = .true.
292 ELSE
293 ijobvr = -1
294 ilvr = .false.
295 END IF
296 ilv = ilvl .OR. ilvr
297*
298* Test the input arguments
299*
300 info = 0
301 lquery = ( lwork.EQ.-1 )
302 IF( ijobvl.LE.0 ) THEN
303 info = -1
304 ELSE IF( ijobvr.LE.0 ) THEN
305 info = -2
306 ELSE IF( n.LT.0 ) THEN
307 info = -3
308 ELSE IF( lda.LT.max( 1, n ) ) THEN
309 info = -5
310 ELSE IF( ldb.LT.max( 1, n ) ) THEN
311 info = -7
312 ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
313 info = -12
314 ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
315 info = -14
316 ELSE IF( lwork.LT.max( 1, 8*n ) .AND. .NOT.lquery ) THEN
317 info = -16
318 END IF
319*
320* Compute workspace
321*
322 IF( info.EQ.0 ) THEN
323 CALL dgeqrf( n, n, b, ldb, work, work, -1, ierr )
324 lwkopt = max(1, 8*n, 3*n+int( work( 1 ) ) )
325 CALL dormqr( 'L', 'T', n, n, n, b, ldb, work, a, lda, work, -1,
326 $ ierr )
327 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
328 IF( ilvl ) THEN
329 CALL dorgqr( n, n, n, vl, ldvl, work, work, -1, ierr )
330 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
331 END IF
332 IF( ilv ) THEN
333 CALL dgghd3( jobvl, jobvr, n, 1, n, a, lda, b, ldb, vl,
334 $ ldvl, vr, ldvr, work, -1, ierr )
335 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
336 CALL dlaqz0( 'S', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
337 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
338 $ work, -1, 0, ierr )
339 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
340 ELSE
341 CALL dgghd3( 'N', 'N', n, 1, n, a, lda, b, ldb, vl, ldvl,
342 $ vr, ldvr, work, -1, ierr )
343 lwkopt = max( lwkopt, 3*n+int( work( 1 ) ) )
344 CALL dlaqz0( 'E', jobvl, jobvr, n, 1, n, a, lda, b, ldb,
345 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
346 $ work, -1, 0, ierr )
347 lwkopt = max( lwkopt, 2*n+int( work( 1 ) ) )
348 END IF
349
350 work( 1 ) = lwkopt
351 END IF
352*
353 IF( info.NE.0 ) THEN
354 CALL xerbla( 'DGGEV3 ', -info )
355 RETURN
356 ELSE IF( lquery ) THEN
357 RETURN
358 END IF
359*
360* Quick return if possible
361*
362 IF( n.EQ.0 )
363 $ RETURN
364*
365* Get machine constants
366*
367 eps = dlamch( 'P' )
368 smlnum = dlamch( 'S' )
369 bignum = one / smlnum
370 CALL dlabad( smlnum, bignum )
371 smlnum = sqrt( smlnum ) / eps
372 bignum = one / smlnum
373*
374* Scale A if max element outside range [SMLNUM,BIGNUM]
375*
376 anrm = dlange( 'M', n, n, a, lda, work )
377 ilascl = .false.
378 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
379 anrmto = smlnum
380 ilascl = .true.
381 ELSE IF( anrm.GT.bignum ) THEN
382 anrmto = bignum
383 ilascl = .true.
384 END IF
385 IF( ilascl )
386 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
387*
388* Scale B if max element outside range [SMLNUM,BIGNUM]
389*
390 bnrm = dlange( 'M', n, n, b, ldb, work )
391 ilbscl = .false.
392 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
393 bnrmto = smlnum
394 ilbscl = .true.
395 ELSE IF( bnrm.GT.bignum ) THEN
396 bnrmto = bignum
397 ilbscl = .true.
398 END IF
399 IF( ilbscl )
400 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
401*
402* Permute the matrices A, B to isolate eigenvalues if possible
403*
404 ileft = 1
405 iright = n + 1
406 iwrk = iright + n
407 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
408 $ work( iright ), work( iwrk ), ierr )
409*
410* Reduce B to triangular form (QR decomposition of B)
411*
412 irows = ihi + 1 - ilo
413 IF( ilv ) THEN
414 icols = n + 1 - ilo
415 ELSE
416 icols = irows
417 END IF
418 itau = iwrk
419 iwrk = itau + irows
420 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
421 $ work( iwrk ), lwork+1-iwrk, ierr )
422*
423* Apply the orthogonal transformation to matrix A
424*
425 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
426 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
427 $ lwork+1-iwrk, ierr )
428*
429* Initialize VL
430*
431 IF( ilvl ) THEN
432 CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
433 IF( irows.GT.1 ) THEN
434 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
435 $ vl( ilo+1, ilo ), ldvl )
436 END IF
437 CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
438 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
439 END IF
440*
441* Initialize VR
442*
443 IF( ilvr )
444 $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
445*
446* Reduce to generalized Hessenberg form
447*
448 IF( ilv ) THEN
449*
450* Eigenvectors requested -- work on whole matrix.
451*
452 CALL dgghd3( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
453 $ ldvl, vr, ldvr, work( iwrk ), lwork+1-iwrk, ierr )
454 ELSE
455 CALL dgghd3( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
456 $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr,
457 $ work( iwrk ), lwork+1-iwrk, ierr )
458 END IF
459*
460* Perform QZ algorithm (Compute eigenvalues, and optionally, the
461* Schur forms and Schur vectors)
462*
463 iwrk = itau
464 IF( ilv ) THEN
465 chtemp = 'S'
466 ELSE
467 chtemp = 'E'
468 END IF
469 CALL dlaqz0( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
470 $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
471 $ work( iwrk ), lwork+1-iwrk, 0, ierr )
472 IF( ierr.NE.0 ) THEN
473 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
474 info = ierr
475 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
476 info = ierr - n
477 ELSE
478 info = n + 1
479 END IF
480 GO TO 110
481 END IF
482*
483* Compute Eigenvectors
484*
485 IF( ilv ) THEN
486 IF( ilvl ) THEN
487 IF( ilvr ) THEN
488 chtemp = 'B'
489 ELSE
490 chtemp = 'L'
491 END IF
492 ELSE
493 chtemp = 'R'
494 END IF
495 CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
496 $ vr, ldvr, n, in, work( iwrk ), ierr )
497 IF( ierr.NE.0 ) THEN
498 info = n + 2
499 GO TO 110
500 END IF
501*
502* Undo balancing on VL and VR and normalization
503*
504 IF( ilvl ) THEN
505 CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
506 $ work( iright ), n, vl, ldvl, ierr )
507 DO 50 jc = 1, n
508 IF( alphai( jc ).LT.zero )
509 $ GO TO 50
510 temp = zero
511 IF( alphai( jc ).EQ.zero ) THEN
512 DO 10 jr = 1, n
513 temp = max( temp, abs( vl( jr, jc ) ) )
514 10 CONTINUE
515 ELSE
516 DO 20 jr = 1, n
517 temp = max( temp, abs( vl( jr, jc ) )+
518 $ abs( vl( jr, jc+1 ) ) )
519 20 CONTINUE
520 END IF
521 IF( temp.LT.smlnum )
522 $ GO TO 50
523 temp = one / temp
524 IF( alphai( jc ).EQ.zero ) THEN
525 DO 30 jr = 1, n
526 vl( jr, jc ) = vl( jr, jc )*temp
527 30 CONTINUE
528 ELSE
529 DO 40 jr = 1, n
530 vl( jr, jc ) = vl( jr, jc )*temp
531 vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
532 40 CONTINUE
533 END IF
534 50 CONTINUE
535 END IF
536 IF( ilvr ) THEN
537 CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
538 $ work( iright ), n, vr, ldvr, ierr )
539 DO 100 jc = 1, n
540 IF( alphai( jc ).LT.zero )
541 $ GO TO 100
542 temp = zero
543 IF( alphai( jc ).EQ.zero ) THEN
544 DO 60 jr = 1, n
545 temp = max( temp, abs( vr( jr, jc ) ) )
546 60 CONTINUE
547 ELSE
548 DO 70 jr = 1, n
549 temp = max( temp, abs( vr( jr, jc ) )+
550 $ abs( vr( jr, jc+1 ) ) )
551 70 CONTINUE
552 END IF
553 IF( temp.LT.smlnum )
554 $ GO TO 100
555 temp = one / temp
556 IF( alphai( jc ).EQ.zero ) THEN
557 DO 80 jr = 1, n
558 vr( jr, jc ) = vr( jr, jc )*temp
559 80 CONTINUE
560 ELSE
561 DO 90 jr = 1, n
562 vr( jr, jc ) = vr( jr, jc )*temp
563 vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
564 90 CONTINUE
565 END IF
566 100 CONTINUE
567 END IF
568*
569* End of eigenvector calculation
570*
571 END IF
572*
573* Undo scaling if necessary
574*
575 110 CONTINUE
576*
577 IF( ilascl ) THEN
578 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
579 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
580 END IF
581*
582 IF( ilbscl ) THEN
583 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
584 END IF
585*
586 work( 1 ) = lwkopt
587 RETURN
588*
589* End of DGGEV3
590*
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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 dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
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
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
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:146
subroutine dtgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTGEVC
Definition: dtgevc.f:295
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 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
subroutine dgghd3(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DGGHD3
Definition: dgghd3.f:230
Here is the call graph for this function:
Here is the caller graph for this function: