LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
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 SHGEQZ,
201 *> =N+2: error return from STGEVC.
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 *> \date April 2012
213 *
214 *> \ingroup complexGEeigen
215 *
216 * =====================================================================
217  SUBROUTINE cggev( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
218  $ vl, ldvl, vr, ldvr, work, lwork, rwork, info )
219 *
220 * -- LAPACK driver routine (version 3.4.1) --
221 * -- LAPACK is a software package provided by Univ. of Tennessee, --
222 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
223 * April 2012
224 *
225 * .. Scalar Arguments ..
226  CHARACTER jobvl, jobvr
227  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
228 * ..
229 * .. Array Arguments ..
230  REAL rwork( * )
231  COMPLEX a( lda, * ), alpha( * ), b( ldb, * ),
232  $ beta( * ), vl( ldvl, * ), vr( ldvr, * ),
233  $ work( * )
234 * ..
235 *
236 * =====================================================================
237 *
238 * .. Parameters ..
239  REAL zero, one
240  parameter( zero = 0.0e0, one = 1.0e0 )
241  COMPLEX czero, cone
242  parameter( czero = ( 0.0e0, 0.0e0 ),
243  $ cone = ( 1.0e0, 0.0e0 ) )
244 * ..
245 * .. Local Scalars ..
246  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery
247  CHARACTER chtemp
248  INTEGER icols, ierr, ihi, ijobvl, ijobvr, ileft, ilo,
249  $ in, iright, irows, irwrk, itau, iwrk, jc, jr,
250  $ lwkmin, lwkopt
251  REAL anrm, anrmto, bignum, bnrm, bnrmto, eps,
252  $ smlnum, temp
253  COMPLEX x
254 * ..
255 * .. Local Arrays ..
256  LOGICAL ldumma( 1 )
257 * ..
258 * .. External Subroutines ..
259  EXTERNAL cgeqrf, cggbak, cggbal, cgghrd, chgeqz, clacpy,
261  $ xerbla
262 * ..
263 * .. External Functions ..
264  LOGICAL lsame
265  INTEGER ilaenv
266  REAL clange, slamch
267  EXTERNAL lsame, ilaenv, clange, slamch
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC abs, aimag, max, REAL, sqrt
271 * ..
272 * .. Statement Functions ..
273  REAL abs1
274 * ..
275 * .. Statement Function definitions ..
276  abs1( x ) = abs( REAL( X ) ) + abs( aimag( x ) )
277 * ..
278 * .. Executable Statements ..
279 *
280 * Decode the input arguments
281 *
282  IF( lsame( jobvl, 'N' ) ) THEN
283  ijobvl = 1
284  ilvl = .false.
285  ELSE IF( lsame( jobvl, 'V' ) ) THEN
286  ijobvl = 2
287  ilvl = .true.
288  ELSE
289  ijobvl = -1
290  ilvl = .false.
291  END IF
292 *
293  IF( lsame( jobvr, 'N' ) ) THEN
294  ijobvr = 1
295  ilvr = .false.
296  ELSE IF( lsame( jobvr, 'V' ) ) THEN
297  ijobvr = 2
298  ilvr = .true.
299  ELSE
300  ijobvr = -1
301  ilvr = .false.
302  END IF
303  ilv = ilvl .OR. ilvr
304 *
305 * Test the input arguments
306 *
307  info = 0
308  lquery = ( lwork.EQ.-1 )
309  IF( ijobvl.LE.0 ) THEN
310  info = -1
311  ELSE IF( ijobvr.LE.0 ) THEN
312  info = -2
313  ELSE IF( n.LT.0 ) THEN
314  info = -3
315  ELSE IF( lda.LT.max( 1, n ) ) THEN
316  info = -5
317  ELSE IF( ldb.LT.max( 1, n ) ) THEN
318  info = -7
319  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
320  info = -11
321  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
322  info = -13
323  END IF
324 *
325 * Compute workspace
326 * (Note: Comments in the code beginning "Workspace:" describe the
327 * minimal amount of workspace needed at that point in the code,
328 * as well as the preferred amount for good performance.
329 * NB refers to the optimal block size for the immediately
330 * following subroutine, as returned by ILAENV. The workspace is
331 * computed assuming ILO = 1 and IHI = N, the worst case.)
332 *
333  IF( info.EQ.0 ) THEN
334  lwkmin = max( 1, 2*n )
335  lwkopt = max( 1, n + n*ilaenv( 1, 'CGEQRF', ' ', n, 1, n, 0 ) )
336  lwkopt = max( lwkopt, n +
337  $ n*ilaenv( 1, 'CUNMQR', ' ', n, 1, n, 0 ) )
338  IF( ilvl ) THEN
339  lwkopt = max( lwkopt, n +
340  $ n*ilaenv( 1, 'CUNGQR', ' ', n, 1, n, -1 ) )
341  END IF
342  work( 1 ) = lwkopt
343 *
344  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
345  $ info = -15
346  END IF
347 *
348  IF( info.NE.0 ) THEN
349  CALL xerbla( 'CGGEV ', -info )
350  return
351  ELSE IF( lquery ) THEN
352  return
353  END IF
354 *
355 * Quick return if possible
356 *
357  IF( n.EQ.0 )
358  $ return
359 *
360 * Get machine constants
361 *
362  eps = slamch( 'E' )*slamch( 'B' )
363  smlnum = slamch( 'S' )
364  bignum = one / smlnum
365  CALL slabad( smlnum, bignum )
366  smlnum = sqrt( smlnum ) / eps
367  bignum = one / smlnum
368 *
369 * Scale A if max element outside range [SMLNUM,BIGNUM]
370 *
371  anrm = clange( 'M', n, n, a, lda, rwork )
372  ilascl = .false.
373  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
374  anrmto = smlnum
375  ilascl = .true.
376  ELSE IF( anrm.GT.bignum ) THEN
377  anrmto = bignum
378  ilascl = .true.
379  END IF
380  IF( ilascl )
381  $ CALL clascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
382 *
383 * Scale B if max element outside range [SMLNUM,BIGNUM]
384 *
385  bnrm = clange( 'M', n, n, b, ldb, rwork )
386  ilbscl = .false.
387  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
388  bnrmto = smlnum
389  ilbscl = .true.
390  ELSE IF( bnrm.GT.bignum ) THEN
391  bnrmto = bignum
392  ilbscl = .true.
393  END IF
394  IF( ilbscl )
395  $ CALL clascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
396 *
397 * Permute the matrices A, B to isolate eigenvalues if possible
398 * (Real Workspace: need 6*N)
399 *
400  ileft = 1
401  iright = n + 1
402  irwrk = iright + n
403  CALL cggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
404  $ rwork( iright ), rwork( irwrk ), ierr )
405 *
406 * Reduce B to triangular form (QR decomposition of B)
407 * (Complex Workspace: need N, prefer N*NB)
408 *
409  irows = ihi + 1 - ilo
410  IF( ilv ) THEN
411  icols = n + 1 - ilo
412  ELSE
413  icols = irows
414  END IF
415  itau = 1
416  iwrk = itau + irows
417  CALL cgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
418  $ work( iwrk ), lwork+1-iwrk, ierr )
419 *
420 * Apply the orthogonal transformation to matrix A
421 * (Complex Workspace: need N, prefer N*NB)
422 *
423  CALL cunmqr( 'L', 'C', 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 * (Complex Workspace: need N, prefer N*NB)
429 *
430  IF( ilvl ) THEN
431  CALL claset( 'Full', n, n, czero, cone, vl, ldvl )
432  IF( irows.GT.1 ) THEN
433  CALL clacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
434  $ vl( ilo+1, ilo ), ldvl )
435  END IF
436  CALL cungqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
437  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
438  END IF
439 *
440 * Initialize VR
441 *
442  IF( ilvr )
443  $ CALL claset( 'Full', n, n, czero, cone, vr, ldvr )
444 *
445 * Reduce to generalized Hessenberg form
446 *
447  IF( ilv ) THEN
448 *
449 * Eigenvectors requested -- work on whole matrix.
450 *
451  CALL cgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
452  $ ldvl, vr, ldvr, ierr )
453  ELSE
454  CALL cgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
455  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
456  END IF
457 *
458 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
459 * Schur form and Schur vectors)
460 * (Complex Workspace: need N)
461 * (Real Workspace: need N)
462 *
463  iwrk = itau
464  IF( ilv ) THEN
465  chtemp = 'S'
466  ELSE
467  chtemp = 'E'
468  END IF
469  CALL chgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
470  $ alpha, beta, vl, ldvl, vr, ldvr, work( iwrk ),
471  $ lwork+1-iwrk, rwork( irwrk ), 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 70
481  END IF
482 *
483 * Compute Eigenvectors
484 * (Real Workspace: need 2*N)
485 * (Complex Workspace: need 2*N)
486 *
487  IF( ilv ) THEN
488  IF( ilvl ) THEN
489  IF( ilvr ) THEN
490  chtemp = 'B'
491  ELSE
492  chtemp = 'L'
493  END IF
494  ELSE
495  chtemp = 'R'
496  END IF
497 *
498  CALL ctgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
499  $ vr, ldvr, n, in, work( iwrk ), rwork( irwrk ),
500  $ ierr )
501  IF( ierr.NE.0 ) THEN
502  info = n + 2
503  go to 70
504  END IF
505 *
506 * Undo balancing on VL and VR and normalization
507 * (Workspace: none needed)
508 *
509  IF( ilvl ) THEN
510  CALL cggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
511  $ rwork( iright ), n, vl, ldvl, ierr )
512  DO 30 jc = 1, n
513  temp = zero
514  DO 10 jr = 1, n
515  temp = max( temp, abs1( vl( jr, jc ) ) )
516  10 continue
517  IF( temp.LT.smlnum )
518  $ go to 30
519  temp = one / temp
520  DO 20 jr = 1, n
521  vl( jr, jc ) = vl( jr, jc )*temp
522  20 continue
523  30 continue
524  END IF
525  IF( ilvr ) THEN
526  CALL cggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
527  $ rwork( iright ), n, vr, ldvr, ierr )
528  DO 60 jc = 1, n
529  temp = zero
530  DO 40 jr = 1, n
531  temp = max( temp, abs1( vr( jr, jc ) ) )
532  40 continue
533  IF( temp.LT.smlnum )
534  $ go to 60
535  temp = one / temp
536  DO 50 jr = 1, n
537  vr( jr, jc ) = vr( jr, jc )*temp
538  50 continue
539  60 continue
540  END IF
541  END IF
542 *
543 * Undo scaling if necessary
544 *
545  70 continue
546 *
547  IF( ilascl )
548  $ CALL clascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
549 *
550  IF( ilbscl )
551  $ CALL clascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
552 *
553  work( 1 ) = lwkopt
554  return
555 *
556 * End of CGGEV
557 *
558  END