LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dgegs.f
Go to the documentation of this file.
1 *> \brief <b> DGEEVX 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 DGEGS + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgegs.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgegs.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgegs.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
22 * ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
23 * LWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBVSL, JOBVSR
27 * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
28 * ..
29 * .. Array Arguments ..
30 * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
31 * $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
32 * $ VSR( LDVSR, * ), WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> This routine is deprecated and has been replaced by routine DGGES.
42 *>
43 *> DGEGS computes the eigenvalues, real Schur form, and, optionally,
44 *> left and or/right Schur vectors of a real matrix pair (A,B).
45 *> Given two square matrices A and B, the generalized real Schur
46 *> factorization has the form
47 *>
48 *> A = Q*S*Z**T, B = Q*T*Z**T
49 *>
50 *> where Q and Z are orthogonal matrices, T is upper triangular, and S
51 *> is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
52 *> blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
53 *> of eigenvalues of (A,B). The columns of Q are the left Schur vectors
54 *> and the columns of Z are the right Schur vectors.
55 *>
56 *> If only the eigenvalues of (A,B) are needed, the driver routine
57 *> DGEGV should be used instead. See DGEGV for a description of the
58 *> eigenvalues of the generalized nonsymmetric eigenvalue problem
59 *> (GNEP).
60 *> \endverbatim
61 *
62 * Arguments:
63 * ==========
64 *
65 *> \param[in] JOBVSL
66 *> \verbatim
67 *> JOBVSL is CHARACTER*1
68 *> = 'N': do not compute the left Schur vectors;
69 *> = 'V': compute the left Schur vectors (returned in VSL).
70 *> \endverbatim
71 *>
72 *> \param[in] JOBVSR
73 *> \verbatim
74 *> JOBVSR is CHARACTER*1
75 *> = 'N': do not compute the right Schur vectors;
76 *> = 'V': compute the right Schur vectors (returned in VSR).
77 *> \endverbatim
78 *>
79 *> \param[in] N
80 *> \verbatim
81 *> N is INTEGER
82 *> The order of the matrices A, B, VSL, and VSR. N >= 0.
83 *> \endverbatim
84 *>
85 *> \param[in,out] A
86 *> \verbatim
87 *> A is DOUBLE PRECISION array, dimension (LDA, N)
88 *> On entry, the matrix A.
89 *> On exit, the upper quasi-triangular matrix S from the
90 *> generalized real Schur factorization.
91 *> \endverbatim
92 *>
93 *> \param[in] LDA
94 *> \verbatim
95 *> LDA is INTEGER
96 *> The leading dimension of A. LDA >= max(1,N).
97 *> \endverbatim
98 *>
99 *> \param[in,out] B
100 *> \verbatim
101 *> B is DOUBLE PRECISION array, dimension (LDB, N)
102 *> On entry, the matrix B.
103 *> On exit, the upper triangular matrix T from the generalized
104 *> real Schur factorization.
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] ALPHAR
114 *> \verbatim
115 *> ALPHAR is DOUBLE PRECISION array, dimension (N)
116 *> The real parts of each scalar alpha defining an eigenvalue
117 *> of GNEP.
118 *> \endverbatim
119 *>
120 *> \param[out] ALPHAI
121 *> \verbatim
122 *> ALPHAI is DOUBLE PRECISION array, dimension (N)
123 *> The imaginary parts of each scalar alpha defining an
124 *> eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
125 *> eigenvalue is real; if positive, then the j-th and (j+1)-st
126 *> eigenvalues are a complex conjugate pair, with
127 *> ALPHAI(j+1) = -ALPHAI(j).
128 *> \endverbatim
129 *>
130 *> \param[out] BETA
131 *> \verbatim
132 *> BETA is DOUBLE PRECISION array, dimension (N)
133 *> The scalars beta that define the eigenvalues of GNEP.
134 *> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
135 *> beta = BETA(j) represent the j-th eigenvalue of the matrix
136 *> pair (A,B), in one of the forms lambda = alpha/beta or
137 *> mu = beta/alpha. Since either lambda or mu may overflow,
138 *> they should not, in general, be computed.
139 *> \endverbatim
140 *>
141 *> \param[out] VSL
142 *> \verbatim
143 *> VSL is DOUBLE PRECISION array, dimension (LDVSL,N)
144 *> If JOBVSL = 'V', the matrix of left Schur vectors Q.
145 *> Not referenced if JOBVSL = 'N'.
146 *> \endverbatim
147 *>
148 *> \param[in] LDVSL
149 *> \verbatim
150 *> LDVSL is INTEGER
151 *> The leading dimension of the matrix VSL. LDVSL >=1, and
152 *> if JOBVSL = 'V', LDVSL >= N.
153 *> \endverbatim
154 *>
155 *> \param[out] VSR
156 *> \verbatim
157 *> VSR is DOUBLE PRECISION array, dimension (LDVSR,N)
158 *> If JOBVSR = 'V', the matrix of right Schur vectors Z.
159 *> Not referenced if JOBVSR = 'N'.
160 *> \endverbatim
161 *>
162 *> \param[in] LDVSR
163 *> \verbatim
164 *> LDVSR is INTEGER
165 *> The leading dimension of the matrix VSR. LDVSR >= 1, and
166 *> if JOBVSR = 'V', LDVSR >= N.
167 *> \endverbatim
168 *>
169 *> \param[out] WORK
170 *> \verbatim
171 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
172 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
173 *> \endverbatim
174 *>
175 *> \param[in] LWORK
176 *> \verbatim
177 *> LWORK is INTEGER
178 *> The dimension of the array WORK. LWORK >= max(1,4*N).
179 *> For good performance, LWORK must generally be larger.
180 *> To compute the optimal value of LWORK, call ILAENV to get
181 *> blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
182 *> NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
183 *> The optimal LWORK is 2*N + N*(NB+1).
184 *>
185 *> If LWORK = -1, then a workspace query is assumed; the routine
186 *> only calculates the optimal size of the WORK array, returns
187 *> this value as the first entry of the WORK array, and no error
188 *> message related to LWORK is issued by XERBLA.
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. (A,B) are not in Schur
198 *> form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
199 *> be correct for j=INFO+1,...,N.
200 *> > N: errors that usually indicate LAPACK problems:
201 *> =N+1: error return from DGGBAL
202 *> =N+2: error return from DGEQRF
203 *> =N+3: error return from DORMQR
204 *> =N+4: error return from DORGQR
205 *> =N+5: error return from DGGHRD
206 *> =N+6: error return from DHGEQZ (other than failed
207 *> iteration)
208 *> =N+7: error return from DGGBAK (computing VSL)
209 *> =N+8: error return from DGGBAK (computing VSR)
210 *> =N+9: error return from DLASCL (various places)
211 *> \endverbatim
212 *
213 * Authors:
214 * ========
215 *
216 *> \author Univ. of Tennessee
217 *> \author Univ. of California Berkeley
218 *> \author Univ. of Colorado Denver
219 *> \author NAG Ltd.
220 *
221 *> \date November 2011
222 *
223 *> \ingroup doubleGEeigen
224 *
225 * =====================================================================
226  SUBROUTINE dgegs( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
227  $ alphai, beta, vsl, ldvsl, vsr, ldvsr, work,
228  $ lwork, info )
229 *
230 * -- LAPACK driver routine (version 3.4.0) --
231 * -- LAPACK is a software package provided by Univ. of Tennessee, --
232 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
233 * November 2011
234 *
235 * .. Scalar Arguments ..
236  CHARACTER jobvsl, jobvsr
237  INTEGER info, lda, ldb, ldvsl, ldvsr, lwork, n
238 * ..
239 * .. Array Arguments ..
240  DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
241  $ b( ldb, * ), beta( * ), vsl( ldvsl, * ),
242  $ vsr( ldvsr, * ), work( * )
243 * ..
244 *
245 * =====================================================================
246 *
247 * .. Parameters ..
248  DOUBLE PRECISION zero, one
249  parameter( zero = 0.0d0, one = 1.0d0 )
250 * ..
251 * .. Local Scalars ..
252  LOGICAL ilascl, ilbscl, ilvsl, ilvsr, lquery
253  INTEGER icols, ihi, iinfo, ijobvl, ijobvr, ileft, ilo,
254  $ iright, irows, itau, iwork, lopt, lwkmin,
255  $ lwkopt, nb, nb1, nb2, nb3
256  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps,
257  $ safmin, smlnum
258 * ..
259 * .. External Subroutines ..
260  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlacpy,
262 * ..
263 * .. External Functions ..
264  LOGICAL lsame
265  INTEGER ilaenv
266  DOUBLE PRECISION dlamch, dlange
267  EXTERNAL lsame, ilaenv, dlamch, dlange
268 * ..
269 * .. Intrinsic Functions ..
270  INTRINSIC int, max
271 * ..
272 * .. Executable Statements ..
273 *
274 * Decode the input arguments
275 *
276  IF( lsame( jobvsl, 'N' ) ) THEN
277  ijobvl = 1
278  ilvsl = .false.
279  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
280  ijobvl = 2
281  ilvsl = .true.
282  ELSE
283  ijobvl = -1
284  ilvsl = .false.
285  END IF
286 *
287  IF( lsame( jobvsr, 'N' ) ) THEN
288  ijobvr = 1
289  ilvsr = .false.
290  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
291  ijobvr = 2
292  ilvsr = .true.
293  ELSE
294  ijobvr = -1
295  ilvsr = .false.
296  END IF
297 *
298 * Test the input arguments
299 *
300  lwkmin = max( 4*n, 1 )
301  lwkopt = lwkmin
302  work( 1 ) = lwkopt
303  lquery = ( lwork.EQ.-1 )
304  info = 0
305  IF( ijobvl.LE.0 ) THEN
306  info = -1
307  ELSE IF( ijobvr.LE.0 ) THEN
308  info = -2
309  ELSE IF( n.LT.0 ) THEN
310  info = -3
311  ELSE IF( lda.LT.max( 1, n ) ) THEN
312  info = -5
313  ELSE IF( ldb.LT.max( 1, n ) ) THEN
314  info = -7
315  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
316  info = -12
317  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
318  info = -14
319  ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
320  info = -16
321  END IF
322 *
323  IF( info.EQ.0 ) THEN
324  nb1 = ilaenv( 1, 'DGEQRF', ' ', n, n, -1, -1 )
325  nb2 = ilaenv( 1, 'DORMQR', ' ', n, n, n, -1 )
326  nb3 = ilaenv( 1, 'DORGQR', ' ', n, n, n, -1 )
327  nb = max( nb1, nb2, nb3 )
328  lopt = 2*n + n*( nb+1 )
329  work( 1 ) = lopt
330  END IF
331 *
332  IF( info.NE.0 ) THEN
333  CALL xerbla( 'DGEGS ', -info )
334  RETURN
335  ELSE IF( lquery ) THEN
336  RETURN
337  END IF
338 *
339 * Quick return if possible
340 *
341  IF( n.EQ.0 )
342  $ RETURN
343 *
344 * Get machine constants
345 *
346  eps = dlamch( 'E' )*dlamch( 'B' )
347  safmin = dlamch( 'S' )
348  smlnum = n*safmin / eps
349  bignum = one / smlnum
350 *
351 * Scale A if max element outside range [SMLNUM,BIGNUM]
352 *
353  anrm = dlange( 'M', n, n, a, lda, work )
354  ilascl = .false.
355  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
356  anrmto = smlnum
357  ilascl = .true.
358  ELSE IF( anrm.GT.bignum ) THEN
359  anrmto = bignum
360  ilascl = .true.
361  END IF
362 *
363  IF( ilascl ) THEN
364  CALL dlascl( 'G', -1, -1, anrm, anrmto, n, n, a, lda, iinfo )
365  IF( iinfo.NE.0 ) THEN
366  info = n + 9
367  RETURN
368  END IF
369  END IF
370 *
371 * Scale B if max element outside range [SMLNUM,BIGNUM]
372 *
373  bnrm = dlange( 'M', n, n, b, ldb, work )
374  ilbscl = .false.
375  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
376  bnrmto = smlnum
377  ilbscl = .true.
378  ELSE IF( bnrm.GT.bignum ) THEN
379  bnrmto = bignum
380  ilbscl = .true.
381  END IF
382 *
383  IF( ilbscl ) THEN
384  CALL dlascl( 'G', -1, -1, bnrm, bnrmto, n, n, b, ldb, iinfo )
385  IF( iinfo.NE.0 ) THEN
386  info = n + 9
387  RETURN
388  END IF
389  END IF
390 *
391 * Permute the matrix to make it more nearly triangular
392 * Workspace layout: (2*N words -- "work..." not actually used)
393 * left_permutation, right_permutation, work...
394 *
395  ileft = 1
396  iright = n + 1
397  iwork = iright + n
398  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
399  $ work( iright ), work( iwork ), iinfo )
400  IF( iinfo.NE.0 ) THEN
401  info = n + 1
402  go to 10
403  END IF
404 *
405 * Reduce B to triangular form, and initialize VSL and/or VSR
406 * Workspace layout: ("work..." must have at least N words)
407 * left_permutation, right_permutation, tau, work...
408 *
409  irows = ihi + 1 - ilo
410  icols = n + 1 - ilo
411  itau = iwork
412  iwork = itau + irows
413  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
414  $ work( iwork ), lwork+1-iwork, iinfo )
415  IF( iinfo.GE.0 )
416  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
417  IF( iinfo.NE.0 ) THEN
418  info = n + 2
419  go to 10
420  END IF
421 *
422  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
423  $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
424  $ lwork+1-iwork, iinfo )
425  IF( iinfo.GE.0 )
426  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
427  IF( iinfo.NE.0 ) THEN
428  info = n + 3
429  go to 10
430  END IF
431 *
432  IF( ilvsl ) THEN
433  CALL dlaset( 'Full', n, n, zero, one, vsl, ldvsl )
434  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
435  $ vsl( ilo+1, ilo ), ldvsl )
436  CALL dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
437  $ work( itau ), work( iwork ), lwork+1-iwork,
438  $ iinfo )
439  IF( iinfo.GE.0 )
440  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
441  IF( iinfo.NE.0 ) THEN
442  info = n + 4
443  go to 10
444  END IF
445  END IF
446 *
447  IF( ilvsr )
448  $ CALL dlaset( 'Full', n, n, zero, one, vsr, ldvsr )
449 *
450 * Reduce to generalized Hessenberg form
451 *
452  CALL dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
453  $ ldvsl, vsr, ldvsr, iinfo )
454  IF( iinfo.NE.0 ) THEN
455  info = n + 5
456  go to 10
457  END IF
458 *
459 * Perform QZ algorithm, computing Schur vectors if desired
460 * Workspace layout: ("work..." must have at least 1 word)
461 * left_permutation, right_permutation, work...
462 *
463  iwork = itau
464  CALL dhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
465  $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
466  $ work( iwork ), lwork+1-iwork, iinfo )
467  IF( iinfo.GE.0 )
468  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
469  IF( iinfo.NE.0 ) THEN
470  IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
471  info = iinfo
472  ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
473  info = iinfo - n
474  ELSE
475  info = n + 6
476  END IF
477  go to 10
478  END IF
479 *
480 * Apply permutation to VSL and VSR
481 *
482  IF( ilvsl ) THEN
483  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
484  $ work( iright ), n, vsl, ldvsl, iinfo )
485  IF( iinfo.NE.0 ) THEN
486  info = n + 7
487  go to 10
488  END IF
489  END IF
490  IF( ilvsr ) THEN
491  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
492  $ work( iright ), n, vsr, ldvsr, iinfo )
493  IF( iinfo.NE.0 ) THEN
494  info = n + 8
495  go to 10
496  END IF
497  END IF
498 *
499 * Undo scaling
500 *
501  IF( ilascl ) THEN
502  CALL dlascl( 'H', -1, -1, anrmto, anrm, n, n, a, lda, iinfo )
503  IF( iinfo.NE.0 ) THEN
504  info = n + 9
505  RETURN
506  END IF
507  CALL dlascl( 'G', -1, -1, anrmto, anrm, n, 1, alphar, n,
508  $ iinfo )
509  IF( iinfo.NE.0 ) THEN
510  info = n + 9
511  RETURN
512  END IF
513  CALL dlascl( 'G', -1, -1, anrmto, anrm, n, 1, alphai, n,
514  $ iinfo )
515  IF( iinfo.NE.0 ) THEN
516  info = n + 9
517  RETURN
518  END IF
519  END IF
520 *
521  IF( ilbscl ) THEN
522  CALL dlascl( 'U', -1, -1, bnrmto, bnrm, n, n, b, ldb, iinfo )
523  IF( iinfo.NE.0 ) THEN
524  info = n + 9
525  RETURN
526  END IF
527  CALL dlascl( 'G', -1, -1, bnrmto, bnrm, n, 1, beta, n, iinfo )
528  IF( iinfo.NE.0 ) THEN
529  info = n + 9
530  RETURN
531  END IF
532  END IF
533 *
534  10 CONTINUE
535  work( 1 ) = lwkopt
536 *
537  RETURN
538 *
539 * End of DGEGS
540 *
541  END