LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
cggev.f
Go to the documentation of this file.
1 *> \brief <b> CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CGGEV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cggev.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cggev.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cggev.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
22 * VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBVL, JOBVR
26 * INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * REAL RWORK( * )
30 * COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
31 * $ BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
32 * $ WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> CGGEV computes for a pair of N-by-N complex nonsymmetric matrices
42 *> (A,B), the generalized eigenvalues, and optionally, the left and/or
43 *> right generalized eigenvectors.
44 *>
45 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
46 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
47 *> singular. It is usually represented as the pair (alpha,beta), as
48 *> there is a reasonable interpretation for beta=0, and even for both
49 *> being zero.
50 *>
51 *> The right generalized eigenvector v(j) corresponding to the
52 *> generalized eigenvalue lambda(j) of (A,B) satisfies
53 *>
54 *> A * v(j) = lambda(j) * B * v(j).
55 *>
56 *> The left generalized eigenvector u(j) corresponding to the
57 *> generalized eigenvalues lambda(j) of (A,B) satisfies
58 *>
59 *> u(j)**H * A = lambda(j) * u(j)**H * B
60 *>
61 *> where u(j)**H is the conjugate-transpose of u(j).
62 *> \endverbatim
63 *
64 * Arguments:
65 * ==========
66 *
67 *> \param[in] JOBVL
68 *> \verbatim
69 *> JOBVL is CHARACTER*1
70 *> = 'N': do not compute the left generalized eigenvectors;
71 *> = 'V': compute the left generalized eigenvectors.
72 *> \endverbatim
73 *>
74 *> \param[in] JOBVR
75 *> \verbatim
76 *> JOBVR is CHARACTER*1
77 *> = 'N': do not compute the right generalized eigenvectors;
78 *> = 'V': compute the right generalized eigenvectors.
79 *> \endverbatim
80 *>
81 *> \param[in] N
82 *> \verbatim
83 *> N is INTEGER
84 *> The order of the matrices A, B, VL, and VR. N >= 0.
85 *> \endverbatim
86 *>
87 *> \param[in,out] A
88 *> \verbatim
89 *> A is COMPLEX array, dimension (LDA, N)
90 *> On entry, the matrix A in the pair (A,B).
91 *> On exit, A has been overwritten.
92 *> \endverbatim
93 *>
94 *> \param[in] LDA
95 *> \verbatim
96 *> LDA is INTEGER
97 *> The leading dimension of A. LDA >= max(1,N).
98 *> \endverbatim
99 *>
100 *> \param[in,out] B
101 *> \verbatim
102 *> B is COMPLEX array, dimension (LDB, N)
103 *> On entry, the matrix B in the pair (A,B).
104 *> On exit, B has been overwritten.
105 *> \endverbatim
106 *>
107 *> \param[in] LDB
108 *> \verbatim
109 *> LDB is INTEGER
110 *> The leading dimension of B. LDB >= max(1,N).
111 *> \endverbatim
112 *>
113 *> \param[out] ALPHA
114 *> \verbatim
115 *> ALPHA is COMPLEX array, dimension (N)
116 *> \endverbatim
117 *>
118 *> \param[out] BETA
119 *> \verbatim
120 *> BETA is COMPLEX array, dimension (N)
121 *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
122 *> generalized eigenvalues.
123 *>
124 *> Note: the quotients ALPHA(j)/BETA(j) may easily over- or
125 *> underflow, and BETA(j) may even be zero. Thus, the user
126 *> should avoid naively computing the ratio alpha/beta.
127 *> However, ALPHA will be always less than and usually
128 *> comparable with norm(A) in magnitude, and BETA always less
129 *> than and usually comparable with norm(B).
130 *> \endverbatim
131 *>
132 *> \param[out] VL
133 *> \verbatim
134 *> VL is COMPLEX array, dimension (LDVL,N)
135 *> If JOBVL = 'V', the left generalized eigenvectors u(j) are
136 *> stored one after another in the columns of VL, in the same
137 *> order as their eigenvalues.
138 *> Each eigenvector is scaled so the largest component has
139 *> abs(real part) + abs(imag. part) = 1.
140 *> Not referenced if JOBVL = 'N'.
141 *> \endverbatim
142 *>
143 *> \param[in] LDVL
144 *> \verbatim
145 *> LDVL is INTEGER
146 *> The leading dimension of the matrix VL. LDVL >= 1, and
147 *> if JOBVL = 'V', LDVL >= N.
148 *> \endverbatim
149 *>
150 *> \param[out] VR
151 *> \verbatim
152 *> VR is COMPLEX array, dimension (LDVR,N)
153 *> If JOBVR = 'V', the right generalized eigenvectors v(j) are
154 *> stored one after another in the columns of VR, in the same
155 *> order as their eigenvalues.
156 *> Each eigenvector is scaled so the largest component has
157 *> abs(real part) + abs(imag. part) = 1.
158 *> Not referenced if JOBVR = 'N'.
159 *> \endverbatim
160 *>
161 *> \param[in] LDVR
162 *> \verbatim
163 *> LDVR is INTEGER
164 *> The leading dimension of the matrix VR. LDVR >= 1, and
165 *> if JOBVR = 'V', LDVR >= N.
166 *> \endverbatim
167 *>
168 *> \param[out] WORK
169 *> \verbatim
170 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
171 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
172 *> \endverbatim
173 *>
174 *> \param[in] LWORK
175 *> \verbatim
176 *> LWORK is INTEGER
177 *> The dimension of the array WORK. LWORK >= max(1,2*N).
178 *> For good performance, LWORK must generally be larger.
179 *>
180 *> If LWORK = -1, then a workspace query is assumed; the routine
181 *> only calculates the optimal size of the WORK array, returns
182 *> this value as the first entry of the WORK array, and no error
183 *> message related to LWORK is issued by XERBLA.
184 *> \endverbatim
185 *>
186 *> \param[out] RWORK
187 *> \verbatim
188 *> RWORK is REAL array, dimension (8*N)
189 *> \endverbatim
190 *>
191 *> \param[out] INFO
192 *> \verbatim
193 *> INFO is INTEGER
194 *> = 0: successful exit
195 *> < 0: if INFO = -i, the i-th argument had an illegal value.
196 *> =1,...,N:
197 *> The QZ iteration failed. No eigenvectors have been
198 *> calculated, but ALPHA(j) and BETA(j) should be
199 *> correct for j=INFO+1,...,N.
200 *> > N: =N+1: other then QZ iteration failed in CHGEQZ,
201 *> =N+2: error return from CTGEVC.
202 *> \endverbatim
203 *
204 * Authors:
205 * ========
206 *
207 *> \author Univ. of Tennessee
208 *> \author Univ. of California Berkeley
209 *> \author Univ. of Colorado Denver
210 *> \author NAG Ltd.
211 *
212 *> \ingroup complexGEeigen
213 *
214 * =====================================================================
215  SUBROUTINE cggev( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
216  $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
217 *
218 * -- LAPACK driver routine --
219 * -- LAPACK is a software package provided by Univ. of Tennessee, --
220 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
221 *
222 * .. Scalar Arguments ..
223  CHARACTER JOBVL, JOBVR
224  INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
225 * ..
226 * .. Array Arguments ..
227  REAL RWORK( * )
228  COMPLEX A( LDA, * ), ALPHA( * ), B( LDB, * ),
229  $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
230  $ work( * )
231 * ..
232 *
233 * =====================================================================
234 *
235 * .. Parameters ..
236  REAL ZERO, ONE
237  parameter( zero = 0.0e0, one = 1.0e0 )
238  COMPLEX CZERO, CONE
239  parameter( czero = ( 0.0e0, 0.0e0 ),
240  $ cone = ( 1.0e0, 0.0e0 ) )
241 * ..
242 * .. Local Scalars ..
243  LOGICAL ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
244  CHARACTER CHTEMP
245  INTEGER ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
246  $ in, iright, irows, irwrk, itau, iwrk, jc, jr,
247  $ lwkmin, lwkopt
248  REAL ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
249  $ smlnum, temp
250  COMPLEX X
251 * ..
252 * .. Local Arrays ..
253  LOGICAL LDUMMA( 1 )
254 * ..
255 * .. External Subroutines ..
256  EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
258  $ xerbla
259 * ..
260 * .. External Functions ..
261  LOGICAL LSAME
262  INTEGER ILAENV
263  REAL CLANGE, SLAMCH
264  EXTERNAL lsame, ilaenv, clange, slamch
265 * ..
266 * .. Intrinsic Functions ..
267  INTRINSIC abs, aimag, max, real, sqrt
268 * ..
269 * .. Statement Functions ..
270  REAL ABS1
271 * ..
272 * .. Statement Function definitions ..
273  abs1( x ) = abs( real( x ) ) + abs( aimag( x ) )
274 * ..
275 * .. Executable Statements ..
276 *
277 * Decode the input arguments
278 *
279  IF( lsame( jobvl, 'N' ) ) THEN
280  ijobvl = 1
281  ilvl = .false.
282  ELSE IF( lsame( jobvl, 'V' ) ) THEN
283  ijobvl = 2
284  ilvl = .true.
285  ELSE
286  ijobvl = -1
287  ilvl = .false.
288  END IF
289 *
290  IF( lsame( jobvr, 'N' ) ) THEN
291  ijobvr = 1
292  ilvr = .false.
293  ELSE IF( lsame( jobvr, 'V' ) ) THEN
294  ijobvr = 2
295  ilvr = .true.
296  ELSE
297  ijobvr = -1
298  ilvr = .false.
299  END IF
300  ilv = ilvl .OR. ilvr
301 *
302 * Test the input arguments
303 *
304  info = 0
305  lquery = ( lwork.EQ.-1 )
306  IF( ijobvl.LE.0 ) THEN
307  info = -1
308  ELSE IF( ijobvr.LE.0 ) THEN
309  info = -2
310  ELSE IF( n.LT.0 ) THEN
311  info = -3
312  ELSE IF( lda.LT.max( 1, n ) ) THEN
313  info = -5
314  ELSE IF( ldb.LT.max( 1, n ) ) THEN
315  info = -7
316  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
317  info = -11
318  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
319  info = -13
320  END IF
321 *
322 * Compute workspace
323 * (Note: Comments in the code beginning "Workspace:" describe the
324 * minimal amount of workspace needed at that point in the code,
325 * as well as the preferred amount for good performance.
326 * NB refers to the optimal block size for the immediately
327 * following subroutine, as returned by ILAENV. The workspace is
328 * computed assuming ILO = 1 and IHI = N, the worst case.)
329 *
330  IF( info.EQ.0 ) THEN
331  lwkmin = max( 1, 2*n )
332  lwkopt = max( 1, n + n*ilaenv( 1, 'CGEQRF', ' ', n, 1, n, 0 ) )
333  lwkopt = max( lwkopt, n +
334  $ n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, 0 ) )
335  IF( ilvl ) THEN
336  lwkopt = max( lwkopt, n +
337  $ n*ilaenv( 1, 'CUNGQR', ' ', n, 1, n, -1 ) )
338  END IF
339  work( 1 ) = lwkopt
340 *
341  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
342  $ info = -15
343  END IF
344 *
345  IF( info.NE.0 ) THEN
346  CALL xerbla( 'CGGEV ', -info )
347  RETURN
348  ELSE IF( lquery ) THEN
349  RETURN
350  END IF
351 *
352 * Quick return if possible
353 *
354  IF( n.EQ.0 )
355  $ RETURN
356 *
357 * Get machine constants
358 *
359  eps = slamch( 'E' )*slamch( 'B' )
360  smlnum = slamch( 'S' )
361  bignum = one / smlnum
362  CALL slabad( smlnum, bignum )
363  smlnum = sqrt( smlnum ) / eps
364  bignum = one / smlnum
365 *
366 * Scale A if max element outside range [SMLNUM,BIGNUM]
367 *
368  anrm = clange( 'M', n, n, a, lda, rwork )
369  ilascl = .false.
370  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
371  anrmto = smlnum
372  ilascl = .true.
373  ELSE IF( anrm.GT.bignum ) THEN
374  anrmto = bignum
375  ilascl = .true.
376  END IF
377  IF( ilascl )
378  $ CALL clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
379 *
380 * Scale B if max element outside range [SMLNUM,BIGNUM]
381 *
382  bnrm = clange( 'M', n, n, b, ldb, rwork )
383  ilbscl = .false.
384  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
385  bnrmto = smlnum
386  ilbscl = .true.
387  ELSE IF( bnrm.GT.bignum ) THEN
388  bnrmto = bignum
389  ilbscl = .true.
390  END IF
391  IF( ilbscl )
392  $ CALL clascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
393 *
394 * Permute the matrices A, B to isolate eigenvalues if possible
395 * (Real Workspace: need 6*N)
396 *
397  ileft = 1
398  iright = n + 1
399  irwrk = iright + n
400  CALL cggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
401  $ rwork( iright ), rwork( irwrk ), ierr )
402 *
403 * Reduce B to triangular form (QR decomposition of B)
404 * (Complex Workspace: need N, prefer N*NB)
405 *
406  irows = ihi + 1 - ilo
407  IF( ilv ) THEN
408  icols = n + 1 - ilo
409  ELSE
410  icols = irows
411  END IF
412  itau = 1
413  iwrk = itau + irows
414  CALL cgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
415  $ work( iwrk ), lwork+1-iwrk, ierr )
416 *
417 * Apply the orthogonal transformation to matrix A
418 * (Complex Workspace: need N, prefer N*NB)
419 *
420  CALL cunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
421  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
422  $ lwork+1-iwrk, ierr )
423 *
424 * Initialize VL
425 * (Complex Workspace: need N, prefer N*NB)
426 *
427  IF( ilvl ) THEN
428  CALL claset( 'Full', n, n, czero, cone, vl, ldvl )
429  IF( irows.GT.1 ) THEN
430  CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
431  $ vl( ilo+1, ilo ), ldvl )
432  END IF
433  CALL cungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
434  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
435  END IF
436 *
437 * Initialize VR
438 *
439  IF( ilvr )
440  $ CALL claset( 'Full', n, n, czero, cone, vr, ldvr )
441 *
442 * Reduce to generalized Hessenberg form
443 *
444  IF( ilv ) THEN
445 *
446 * Eigenvectors requested -- work on whole matrix.
447 *
448  CALL cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
449  $ ldvl, vr, ldvr, ierr )
450  ELSE
451  CALL cgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
452  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
453  END IF
454 *
455 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
456 * Schur form and Schur vectors)
457 * (Complex Workspace: need N)
458 * (Real Workspace: need N)
459 *
460  iwrk = itau
461  IF( ilv ) THEN
462  chtemp = 'S'
463  ELSE
464  chtemp = 'E'
465  END IF
466  CALL chgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
467  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
468  $ lwork+1-iwrk, rwork( irwrk ), ierr )
469  IF( ierr.NE.0 ) THEN
470  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
471  info = ierr
472  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
473  info = ierr - n
474  ELSE
475  info = n + 1
476  END IF
477  GO TO 70
478  END IF
479 *
480 * Compute Eigenvectors
481 * (Real Workspace: need 2*N)
482 * (Complex Workspace: need 2*N)
483 *
484  IF( ilv ) THEN
485  IF( ilvl ) THEN
486  IF( ilvr ) THEN
487  chtemp = 'B'
488  ELSE
489  chtemp = 'L'
490  END IF
491  ELSE
492  chtemp = 'R'
493  END IF
494 *
495  CALL ctgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
496  $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
497  $ ierr )
498  IF( ierr.NE.0 ) THEN
499  info = n + 2
500  GO TO 70
501  END IF
502 *
503 * Undo balancing on VL and VR and normalization
504 * (Workspace: none needed)
505 *
506  IF( ilvl ) THEN
507  CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
508  $ rwork( iright ), n, vl, ldvl, ierr )
509  DO 30 jc = 1, n
510  temp = zero
511  DO 10 jr = 1, n
512  temp = max( temp, abs1( vl( jr, jc ) ) )
513  10 CONTINUE
514  IF( temp.LT.smlnum )
515  $ GO TO 30
516  temp = one / temp
517  DO 20 jr = 1, n
518  vl( jr, jc ) = vl( jr, jc )*temp
519  20 CONTINUE
520  30 CONTINUE
521  END IF
522  IF( ilvr ) THEN
523  CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
524  $ rwork( iright ), n, vr, ldvr, ierr )
525  DO 60 jc = 1, n
526  temp = zero
527  DO 40 jr = 1, n
528  temp = max( temp, abs1( vr( jr, jc ) ) )
529  40 CONTINUE
530  IF( temp.LT.smlnum )
531  $ GO TO 60
532  temp = one / temp
533  DO 50 jr = 1, n
534  vr( jr, jc ) = vr( jr, jc )*temp
535  50 CONTINUE
536  60 CONTINUE
537  END IF
538  END IF
539 *
540 * Undo scaling if necessary
541 *
542  70 CONTINUE
543 *
544  IF( ilascl )
545  $ CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
546 *
547  IF( ilbscl )
548  $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
549 *
550  work( 1 ) = lwkopt
551  RETURN
552 *
553 * End of CGGEV
554 *
555  END
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
CGGBAL
Definition: cggbal.f:177
subroutine cggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
CGGBAK
Definition: cggbak.f:148
subroutine chgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
CHGEQZ
Definition: chgeqz.f:284
subroutine cgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
CGEQRF
Definition: cgeqrf.f:145
subroutine ctgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
CTGEVC
Definition: ctgevc.f:219
subroutine cggev(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO)
CGGEV computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: cggev.f:217
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: claset.f:106
subroutine clascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
CLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: clascl.f:143
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
subroutine cunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
CUNMQR
Definition: cunmqr.f:168
subroutine cgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
CGGHRD
Definition: cgghrd.f:204
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
Definition: cungqr.f:128