LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zgges.f
Go to the documentation of this file.
1 *> \brief <b> ZGGES computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors 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 ZGGES + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgges.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgges.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgges.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGGES( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
22 * SDIM, ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
23 * LWORK, RWORK, BWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBVSL, JOBVSR, SORT
27 * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
28 * ..
29 * .. Array Arguments ..
30 * LOGICAL BWORK( * )
31 * DOUBLE PRECISION RWORK( * )
32 * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
33 * $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
34 * $ WORK( * )
35 * ..
36 * .. Function Arguments ..
37 * LOGICAL SELCTG
38 * EXTERNAL SELCTG
39 * ..
40 *
41 *
42 *> \par Purpose:
43 * =============
44 *>
45 *> \verbatim
46 *>
47 *> ZGGES computes for a pair of N-by-N complex nonsymmetric matrices
48 *> (A,B), the generalized eigenvalues, the generalized complex Schur
49 *> form (S, T), and optionally left and/or right Schur vectors (VSL
50 *> and VSR). This gives the generalized Schur factorization
51 *>
52 *> (A,B) = ( (VSL)*S*(VSR)**H, (VSL)*T*(VSR)**H )
53 *>
54 *> where (VSR)**H is the conjugate-transpose of VSR.
55 *>
56 *> Optionally, it also orders the eigenvalues so that a selected cluster
57 *> of eigenvalues appears in the leading diagonal blocks of the upper
58 *> triangular matrix S and the upper triangular matrix T. The leading
59 *> columns of VSL and VSR then form an unitary basis for the
60 *> corresponding left and right eigenspaces (deflating subspaces).
61 *>
62 *> (If only the generalized eigenvalues are needed, use the driver
63 *> ZGGEV instead, which is faster.)
64 *>
65 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
66 *> or a ratio alpha/beta = w, such that A - w*B is singular. It is
67 *> usually represented as the pair (alpha,beta), as there is a
68 *> reasonable interpretation for beta=0, and even for both being zero.
69 *>
70 *> A pair of matrices (S,T) is in generalized complex Schur form if S
71 *> and T are upper triangular and, in addition, the diagonal elements
72 *> of T are non-negative real numbers.
73 *> \endverbatim
74 *
75 * Arguments:
76 * ==========
77 *
78 *> \param[in] JOBVSL
79 *> \verbatim
80 *> JOBVSL is CHARACTER*1
81 *> = 'N': do not compute the left Schur vectors;
82 *> = 'V': compute the left Schur vectors.
83 *> \endverbatim
84 *>
85 *> \param[in] JOBVSR
86 *> \verbatim
87 *> JOBVSR is CHARACTER*1
88 *> = 'N': do not compute the right Schur vectors;
89 *> = 'V': compute the right Schur vectors.
90 *> \endverbatim
91 *>
92 *> \param[in] SORT
93 *> \verbatim
94 *> SORT is CHARACTER*1
95 *> Specifies whether or not to order the eigenvalues on the
96 *> diagonal of the generalized Schur form.
97 *> = 'N': Eigenvalues are not ordered;
98 *> = 'S': Eigenvalues are ordered (see SELCTG).
99 *> \endverbatim
100 *>
101 *> \param[in] SELCTG
102 *> \verbatim
103 *> SELCTG is a LOGICAL FUNCTION of two COMPLEX*16 arguments
104 *> SELCTG must be declared EXTERNAL in the calling subroutine.
105 *> If SORT = 'N', SELCTG is not referenced.
106 *> If SORT = 'S', SELCTG is used to select eigenvalues to sort
107 *> to the top left of the Schur form.
108 *> An eigenvalue ALPHA(j)/BETA(j) is selected if
109 *> SELCTG(ALPHA(j),BETA(j)) is true.
110 *>
111 *> Note that a selected complex eigenvalue may no longer satisfy
112 *> SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
113 *> ordering may change the value of complex eigenvalues
114 *> (especially if the eigenvalue is ill-conditioned), in this
115 *> case INFO is set to N+2 (See INFO below).
116 *> \endverbatim
117 *>
118 *> \param[in] N
119 *> \verbatim
120 *> N is INTEGER
121 *> The order of the matrices A, B, VSL, and VSR. N >= 0.
122 *> \endverbatim
123 *>
124 *> \param[in,out] A
125 *> \verbatim
126 *> A is COMPLEX*16 array, dimension (LDA, N)
127 *> On entry, the first of the pair of matrices.
128 *> On exit, A has been overwritten by its generalized Schur
129 *> form S.
130 *> \endverbatim
131 *>
132 *> \param[in] LDA
133 *> \verbatim
134 *> LDA is INTEGER
135 *> The leading dimension of A. LDA >= max(1,N).
136 *> \endverbatim
137 *>
138 *> \param[in,out] B
139 *> \verbatim
140 *> B is COMPLEX*16 array, dimension (LDB, N)
141 *> On entry, the second of the pair of matrices.
142 *> On exit, B has been overwritten by its generalized Schur
143 *> form T.
144 *> \endverbatim
145 *>
146 *> \param[in] LDB
147 *> \verbatim
148 *> LDB is INTEGER
149 *> The leading dimension of B. LDB >= max(1,N).
150 *> \endverbatim
151 *>
152 *> \param[out] SDIM
153 *> \verbatim
154 *> SDIM is INTEGER
155 *> If SORT = 'N', SDIM = 0.
156 *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
157 *> for which SELCTG is true.
158 *> \endverbatim
159 *>
160 *> \param[out] ALPHA
161 *> \verbatim
162 *> ALPHA is COMPLEX*16 array, dimension (N)
163 *> \endverbatim
164 *>
165 *> \param[out] BETA
166 *> \verbatim
167 *> BETA is COMPLEX*16 array, dimension (N)
168 *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
169 *> generalized eigenvalues. ALPHA(j), j=1,...,N and BETA(j),
170 *> j=1,...,N are the diagonals of the complex Schur form (A,B)
171 *> output by ZGGES. The BETA(j) will be non-negative real.
172 *>
173 *> Note: the quotients ALPHA(j)/BETA(j) may easily over- or
174 *> underflow, and BETA(j) may even be zero. Thus, the user
175 *> should avoid naively computing the ratio alpha/beta.
176 *> However, ALPHA will be always less than and usually
177 *> comparable with norm(A) in magnitude, and BETA always less
178 *> than and usually comparable with norm(B).
179 *> \endverbatim
180 *>
181 *> \param[out] VSL
182 *> \verbatim
183 *> VSL is COMPLEX*16 array, dimension (LDVSL,N)
184 *> If JOBVSL = 'V', VSL will contain the left Schur vectors.
185 *> Not referenced if JOBVSL = 'N'.
186 *> \endverbatim
187 *>
188 *> \param[in] LDVSL
189 *> \verbatim
190 *> LDVSL is INTEGER
191 *> The leading dimension of the matrix VSL. LDVSL >= 1, and
192 *> if JOBVSL = 'V', LDVSL >= N.
193 *> \endverbatim
194 *>
195 *> \param[out] VSR
196 *> \verbatim
197 *> VSR is COMPLEX*16 array, dimension (LDVSR,N)
198 *> If JOBVSR = 'V', VSR will contain the right Schur vectors.
199 *> Not referenced if JOBVSR = 'N'.
200 *> \endverbatim
201 *>
202 *> \param[in] LDVSR
203 *> \verbatim
204 *> LDVSR is INTEGER
205 *> The leading dimension of the matrix VSR. LDVSR >= 1, and
206 *> if JOBVSR = 'V', LDVSR >= N.
207 *> \endverbatim
208 *>
209 *> \param[out] WORK
210 *> \verbatim
211 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
212 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
213 *> \endverbatim
214 *>
215 *> \param[in] LWORK
216 *> \verbatim
217 *> LWORK is INTEGER
218 *> The dimension of the array WORK. LWORK >= max(1,2*N).
219 *> For good performance, LWORK must generally be larger.
220 *>
221 *> If LWORK = -1, then a workspace query is assumed; the routine
222 *> only calculates the optimal size of the WORK array, returns
223 *> this value as the first entry of the WORK array, and no error
224 *> message related to LWORK is issued by XERBLA.
225 *> \endverbatim
226 *>
227 *> \param[out] RWORK
228 *> \verbatim
229 *> RWORK is DOUBLE PRECISION array, dimension (8*N)
230 *> \endverbatim
231 *>
232 *> \param[out] BWORK
233 *> \verbatim
234 *> BWORK is LOGICAL array, dimension (N)
235 *> Not referenced if SORT = 'N'.
236 *> \endverbatim
237 *>
238 *> \param[out] INFO
239 *> \verbatim
240 *> INFO is INTEGER
241 *> = 0: successful exit
242 *> < 0: if INFO = -i, the i-th argument had an illegal value.
243 *> =1,...,N:
244 *> The QZ iteration failed. (A,B) are not in Schur
245 *> form, but ALPHA(j) and BETA(j) should be correct for
246 *> j=INFO+1,...,N.
247 *> > N: =N+1: other than QZ iteration failed in ZHGEQZ
248 *> =N+2: after reordering, roundoff changed values of
249 *> some complex eigenvalues so that leading
250 *> eigenvalues in the Generalized Schur form no
251 *> longer satisfy SELCTG=.TRUE. This could also
252 *> be caused due to scaling.
253 *> =N+3: reordering falied in ZTGSEN.
254 *> \endverbatim
255 *
256 * Authors:
257 * ========
258 *
259 *> \author Univ. of Tennessee
260 *> \author Univ. of California Berkeley
261 *> \author Univ. of Colorado Denver
262 *> \author NAG Ltd.
263 *
264 *> \date November 2011
265 *
266 *> \ingroup complex16GEeigen
267 *
268 * =====================================================================
269  SUBROUTINE zgges( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B, LDB,
270  $ sdim, alpha, beta, vsl, ldvsl, vsr, ldvsr, work,
271  $ lwork, rwork, bwork, info )
272 *
273 * -- LAPACK driver routine (version 3.4.0) --
274 * -- LAPACK is a software package provided by Univ. of Tennessee, --
275 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
276 * November 2011
277 *
278 * .. Scalar Arguments ..
279  CHARACTER jobvsl, jobvsr, sort
280  INTEGER info, lda, ldb, ldvsl, ldvsr, lwork, n, sdim
281 * ..
282 * .. Array Arguments ..
283  LOGICAL bwork( * )
284  DOUBLE PRECISION rwork( * )
285  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
286  $ beta( * ), vsl( ldvsl, * ), vsr( ldvsr, * ),
287  $ work( * )
288 * ..
289 * .. Function Arguments ..
290  LOGICAL selctg
291  EXTERNAL selctg
292 * ..
293 *
294 * =====================================================================
295 *
296 * .. Parameters ..
297  DOUBLE PRECISION zero, one
298  parameter( zero = 0.0d0, one = 1.0d0 )
299  COMPLEX*16 czero, cone
300  parameter( czero = ( 0.0d0, 0.0d0 ),
301  $ cone = ( 1.0d0, 0.0d0 ) )
302 * ..
303 * .. Local Scalars ..
304  LOGICAL cursl, ilascl, ilbscl, ilvsl, ilvsr, lastsl,
305  $ lquery, wantst
306  INTEGER i, icols, ierr, ihi, ijobvl, ijobvr, ileft,
307  $ ilo, iright, irows, irwrk, itau, iwrk, lwkmin,
308  $ lwkopt
309  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps, pvsl,
310  $ pvsr, smlnum
311 * ..
312 * .. Local Arrays ..
313  INTEGER idum( 1 )
314  DOUBLE PRECISION dif( 2 )
315 * ..
316 * .. External Subroutines ..
317  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghrd,
319  $ zunmqr
320 * ..
321 * .. External Functions ..
322  LOGICAL lsame
323  INTEGER ilaenv
324  DOUBLE PRECISION dlamch, zlange
325  EXTERNAL lsame, ilaenv, dlamch, zlange
326 * ..
327 * .. Intrinsic Functions ..
328  INTRINSIC max, sqrt
329 * ..
330 * .. Executable Statements ..
331 *
332 * Decode the input arguments
333 *
334  IF( lsame( jobvsl, 'N' ) ) THEN
335  ijobvl = 1
336  ilvsl = .false.
337  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
338  ijobvl = 2
339  ilvsl = .true.
340  ELSE
341  ijobvl = -1
342  ilvsl = .false.
343  END IF
344 *
345  IF( lsame( jobvsr, 'N' ) ) THEN
346  ijobvr = 1
347  ilvsr = .false.
348  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
349  ijobvr = 2
350  ilvsr = .true.
351  ELSE
352  ijobvr = -1
353  ilvsr = .false.
354  END IF
355 *
356  wantst = lsame( sort, 'S' )
357 *
358 * Test the input arguments
359 *
360  info = 0
361  lquery = ( lwork.EQ.-1 )
362  IF( ijobvl.LE.0 ) THEN
363  info = -1
364  ELSE IF( ijobvr.LE.0 ) THEN
365  info = -2
366  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
367  info = -3
368  ELSE IF( n.LT.0 ) THEN
369  info = -5
370  ELSE IF( lda.LT.max( 1, n ) ) THEN
371  info = -7
372  ELSE IF( ldb.LT.max( 1, n ) ) THEN
373  info = -9
374  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
375  info = -14
376  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
377  info = -16
378  END IF
379 *
380 * Compute workspace
381 * (Note: Comments in the code beginning "Workspace:" describe the
382 * minimal amount of workspace needed at that point in the code,
383 * as well as the preferred amount for good performance.
384 * NB refers to the optimal block size for the immediately
385 * following subroutine, as returned by ILAENV.)
386 *
387  IF( info.EQ.0 ) THEN
388  lwkmin = max( 1, 2*n )
389  lwkopt = max( 1, n + n*ilaenv( 1, 'ZGEQRF', ' ', n, 1, n, 0 ) )
390  lwkopt = max( lwkopt, n +
391  $ n*ilaenv( 1, 'ZUNMQR', ' ', n, 1, n, -1 ) )
392  IF( ilvsl ) THEN
393  lwkopt = max( lwkopt, n +
394  $ n*ilaenv( 1, 'ZUNGQR', ' ', n, 1, n, -1 ) )
395  END IF
396  work( 1 ) = lwkopt
397 *
398  IF( lwork.LT.lwkmin .AND. .NOT.lquery )
399  $ info = -18
400  END IF
401 *
402  IF( info.NE.0 ) THEN
403  CALL xerbla( 'ZGGES ', -info )
404  return
405  ELSE IF( lquery ) THEN
406  return
407  END IF
408 *
409 * Quick return if possible
410 *
411  IF( n.EQ.0 ) THEN
412  sdim = 0
413  return
414  END IF
415 *
416 * Get machine constants
417 *
418  eps = dlamch( 'P' )
419  smlnum = dlamch( 'S' )
420  bignum = one / smlnum
421  CALL dlabad( smlnum, bignum )
422  smlnum = sqrt( smlnum ) / eps
423  bignum = one / smlnum
424 *
425 * Scale A if max element outside range [SMLNUM,BIGNUM]
426 *
427  anrm = zlange( 'M', n, n, a, lda, rwork )
428  ilascl = .false.
429  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
430  anrmto = smlnum
431  ilascl = .true.
432  ELSE IF( anrm.GT.bignum ) THEN
433  anrmto = bignum
434  ilascl = .true.
435  END IF
436 *
437  IF( ilascl )
438  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
439 *
440 * Scale B if max element outside range [SMLNUM,BIGNUM]
441 *
442  bnrm = zlange( 'M', n, n, b, ldb, rwork )
443  ilbscl = .false.
444  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
445  bnrmto = smlnum
446  ilbscl = .true.
447  ELSE IF( bnrm.GT.bignum ) THEN
448  bnrmto = bignum
449  ilbscl = .true.
450  END IF
451 *
452  IF( ilbscl )
453  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
454 *
455 * Permute the matrix to make it more nearly triangular
456 * (Real Workspace: need 6*N)
457 *
458  ileft = 1
459  iright = n + 1
460  irwrk = iright + n
461  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
462  $ rwork( iright ), rwork( irwrk ), ierr )
463 *
464 * Reduce B to triangular form (QR decomposition of B)
465 * (Complex Workspace: need N, prefer N*NB)
466 *
467  irows = ihi + 1 - ilo
468  icols = n + 1 - ilo
469  itau = 1
470  iwrk = itau + irows
471  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
472  $ work( iwrk ), lwork+1-iwrk, ierr )
473 *
474 * Apply the orthogonal transformation to matrix A
475 * (Complex Workspace: need N, prefer N*NB)
476 *
477  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
478  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
479  $ lwork+1-iwrk, ierr )
480 *
481 * Initialize VSL
482 * (Complex Workspace: need N, prefer N*NB)
483 *
484  IF( ilvsl ) THEN
485  CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
486  IF( irows.GT.1 ) THEN
487  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
488  $ vsl( ilo+1, ilo ), ldvsl )
489  END IF
490  CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
491  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
492  END IF
493 *
494 * Initialize VSR
495 *
496  IF( ilvsr )
497  $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
498 *
499 * Reduce to generalized Hessenberg form
500 * (Workspace: none needed)
501 *
502  CALL zgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
503  $ ldvsl, vsr, ldvsr, ierr )
504 *
505  sdim = 0
506 *
507 * Perform QZ algorithm, computing Schur vectors if desired
508 * (Complex Workspace: need N)
509 * (Real Workspace: need N)
510 *
511  iwrk = itau
512  CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
513  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
514  $ lwork+1-iwrk, rwork( irwrk ), ierr )
515  IF( ierr.NE.0 ) THEN
516  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
517  info = ierr
518  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
519  info = ierr - n
520  ELSE
521  info = n + 1
522  END IF
523  go to 30
524  END IF
525 *
526 * Sort eigenvalues ALPHA/BETA if desired
527 * (Workspace: none needed)
528 *
529  IF( wantst ) THEN
530 *
531 * Undo scaling on eigenvalues before selecting
532 *
533  IF( ilascl )
534  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, 1, alpha, n, ierr )
535  IF( ilbscl )
536  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, 1, beta, n, ierr )
537 *
538 * Select eigenvalues
539 *
540  DO 10 i = 1, n
541  bwork( i ) = selctg( alpha( i ), beta( i ) )
542  10 continue
543 *
544  CALL ztgsen( 0, ilvsl, ilvsr, bwork, n, a, lda, b, ldb, alpha,
545  $ beta, vsl, ldvsl, vsr, ldvsr, sdim, pvsl, pvsr,
546  $ dif, work( iwrk ), lwork-iwrk+1, idum, 1, ierr )
547  IF( ierr.EQ.1 )
548  $ info = n + 3
549 *
550  END IF
551 *
552 * Apply back-permutation to VSL and VSR
553 * (Workspace: none needed)
554 *
555  IF( ilvsl )
556  $ CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
557  $ rwork( iright ), n, vsl, ldvsl, ierr )
558  IF( ilvsr )
559  $ CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
560  $ rwork( iright ), n, vsr, ldvsr, ierr )
561 *
562 * Undo scaling
563 *
564  IF( ilascl ) THEN
565  CALL zlascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
566  CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
567  END IF
568 *
569  IF( ilbscl ) THEN
570  CALL zlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
571  CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
572  END IF
573 *
574  IF( wantst ) THEN
575 *
576 * Check if reordering is correct
577 *
578  lastsl = .true.
579  sdim = 0
580  DO 20 i = 1, n
581  cursl = selctg( alpha( i ), beta( i ) )
582  IF( cursl )
583  $ sdim = sdim + 1
584  IF( cursl .AND. .NOT.lastsl )
585  $ info = n + 2
586  lastsl = cursl
587  20 continue
588 *
589  END IF
590 *
591  30 continue
592 *
593  work( 1 ) = lwkopt
594 *
595  return
596 *
597 * End of ZGGES
598 *
599  END