LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \ingroup doubleGEeigen
222 *
223 * =====================================================================
224  SUBROUTINE dgegs( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
225  $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
226  $ LWORK, INFO )
227 *
228 * -- LAPACK driver routine --
229 * -- LAPACK is a software package provided by Univ. of Tennessee, --
230 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
231 *
232 * .. Scalar Arguments ..
233  CHARACTER JOBVSL, JOBVSR
234  INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
235 * ..
236 * .. Array Arguments ..
237  DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
238  $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
239  $ vsr( ldvsr, * ), work( * )
240 * ..
241 *
242 * =====================================================================
243 *
244 * .. Parameters ..
245  DOUBLE PRECISION ZERO, ONE
246  PARAMETER ( ZERO = 0.0d0, one = 1.0d0 )
247 * ..
248 * .. Local Scalars ..
249  LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
250  INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
251  $ iright, irows, itau, iwork, lopt, lwkmin,
252  $ lwkopt, nb, nb1, nb2, nb3
253  DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
254  $ SAFMIN, SMLNUM
255 * ..
256 * .. External Subroutines ..
257  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlacpy,
259 * ..
260 * .. External Functions ..
261  LOGICAL LSAME
262  INTEGER ILAENV
263  DOUBLE PRECISION DLAMCH, DLANGE
264  EXTERNAL lsame, ilaenv, dlamch, dlange
265 * ..
266 * .. Intrinsic Functions ..
267  INTRINSIC int, max
268 * ..
269 * .. Executable Statements ..
270 *
271 * Decode the input arguments
272 *
273  IF( lsame( jobvsl, 'N' ) ) THEN
274  ijobvl = 1
275  ilvsl = .false.
276  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
277  ijobvl = 2
278  ilvsl = .true.
279  ELSE
280  ijobvl = -1
281  ilvsl = .false.
282  END IF
283 *
284  IF( lsame( jobvsr, 'N' ) ) THEN
285  ijobvr = 1
286  ilvsr = .false.
287  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
288  ijobvr = 2
289  ilvsr = .true.
290  ELSE
291  ijobvr = -1
292  ilvsr = .false.
293  END IF
294 *
295 * Test the input arguments
296 *
297  lwkmin = max( 4*n, 1 )
298  lwkopt = lwkmin
299  work( 1 ) = lwkopt
300  lquery = ( lwork.EQ.-1 )
301  info = 0
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( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
313  info = -12
314  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
315  info = -14
316  ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
317  info = -16
318  END IF
319 *
320  IF( info.EQ.0 ) THEN
321  nb1 = ilaenv( 1, 'DGEQRF', ' ', n, n, -1, -1 )
322  nb2 = ilaenv( 1, 'DORMQR', ' ', n, n, n, -1 )
323  nb3 = ilaenv( 1, 'DORGQR', ' ', n, n, n, -1 )
324  nb = max( nb1, nb2, nb3 )
325  lopt = 2*n + n*( nb+1 )
326  work( 1 ) = lopt
327  END IF
328 *
329  IF( info.NE.0 ) THEN
330  CALL xerbla( 'DGEGS ', -info )
331  RETURN
332  ELSE IF( lquery ) THEN
333  RETURN
334  END IF
335 *
336 * Quick return if possible
337 *
338  IF( n.EQ.0 )
339  $ RETURN
340 *
341 * Get machine constants
342 *
343  eps = dlamch( 'E' )*dlamch( 'B' )
344  safmin = dlamch( 'S' )
345  smlnum = n*safmin / eps
346  bignum = one / smlnum
347 *
348 * Scale A if max element outside range [SMLNUM,BIGNUM]
349 *
350  anrm = dlange( 'M', n, n, a, lda, work )
351  ilascl = .false.
352  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
353  anrmto = smlnum
354  ilascl = .true.
355  ELSE IF( anrm.GT.bignum ) THEN
356  anrmto = bignum
357  ilascl = .true.
358  END IF
359 *
360  IF( ilascl ) THEN
361  CALL dlascl( 'G', -1, -1, anrm, anrmto, n, n, a, lda, iinfo )
362  IF( iinfo.NE.0 ) THEN
363  info = n + 9
364  RETURN
365  END IF
366  END IF
367 *
368 * Scale B if max element outside range [SMLNUM,BIGNUM]
369 *
370  bnrm = dlange( 'M', n, n, b, ldb, work )
371  ilbscl = .false.
372  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
373  bnrmto = smlnum
374  ilbscl = .true.
375  ELSE IF( bnrm.GT.bignum ) THEN
376  bnrmto = bignum
377  ilbscl = .true.
378  END IF
379 *
380  IF( ilbscl ) THEN
381  CALL dlascl( 'G', -1, -1, bnrm, bnrmto, n, n, b, ldb, iinfo )
382  IF( iinfo.NE.0 ) THEN
383  info = n + 9
384  RETURN
385  END IF
386  END IF
387 *
388 * Permute the matrix to make it more nearly triangular
389 * Workspace layout: (2*N words -- "work..." not actually used)
390 * left_permutation, right_permutation, work...
391 *
392  ileft = 1
393  iright = n + 1
394  iwork = iright + n
395  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
396  $ work( iright ), work( iwork ), iinfo )
397  IF( iinfo.NE.0 ) THEN
398  info = n + 1
399  GO TO 10
400  END IF
401 *
402 * Reduce B to triangular form, and initialize VSL and/or VSR
403 * Workspace layout: ("work..." must have at least N words)
404 * left_permutation, right_permutation, tau, work...
405 *
406  irows = ihi + 1 - ilo
407  icols = n + 1 - ilo
408  itau = iwork
409  iwork = itau + irows
410  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
411  $ work( iwork ), lwork+1-iwork, iinfo )
412  IF( iinfo.GE.0 )
413  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
414  IF( iinfo.NE.0 ) THEN
415  info = n + 2
416  GO TO 10
417  END IF
418 *
419  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
420  $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
421  $ lwork+1-iwork, iinfo )
422  IF( iinfo.GE.0 )
423  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
424  IF( iinfo.NE.0 ) THEN
425  info = n + 3
426  GO TO 10
427  END IF
428 *
429  IF( ilvsl ) THEN
430  CALL dlaset( 'Full', n, n, zero, one, vsl, ldvsl )
431  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
432  $ vsl( ilo+1, ilo ), ldvsl )
433  CALL dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
434  $ work( itau ), work( iwork ), lwork+1-iwork,
435  $ iinfo )
436  IF( iinfo.GE.0 )
437  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
438  IF( iinfo.NE.0 ) THEN
439  info = n + 4
440  GO TO 10
441  END IF
442  END IF
443 *
444  IF( ilvsr )
445  $ CALL dlaset( 'Full', n, n, zero, one, vsr, ldvsr )
446 *
447 * Reduce to generalized Hessenberg form
448 *
449  CALL dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
450  $ ldvsl, vsr, ldvsr, iinfo )
451  IF( iinfo.NE.0 ) THEN
452  info = n + 5
453  GO TO 10
454  END IF
455 *
456 * Perform QZ algorithm, computing Schur vectors if desired
457 * Workspace layout: ("work..." must have at least 1 word)
458 * left_permutation, right_permutation, work...
459 *
460  iwork = itau
461  CALL dhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
462  $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
463  $ work( iwork ), lwork+1-iwork, iinfo )
464  IF( iinfo.GE.0 )
465  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
466  IF( iinfo.NE.0 ) THEN
467  IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
468  info = iinfo
469  ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
470  info = iinfo - n
471  ELSE
472  info = n + 6
473  END IF
474  GO TO 10
475  END IF
476 *
477 * Apply permutation to VSL and VSR
478 *
479  IF( ilvsl ) THEN
480  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
481  $ work( iright ), n, vsl, ldvsl, iinfo )
482  IF( iinfo.NE.0 ) THEN
483  info = n + 7
484  GO TO 10
485  END IF
486  END IF
487  IF( ilvsr ) THEN
488  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
489  $ work( iright ), n, vsr, ldvsr, iinfo )
490  IF( iinfo.NE.0 ) THEN
491  info = n + 8
492  GO TO 10
493  END IF
494  END IF
495 *
496 * Undo scaling
497 *
498  IF( ilascl ) THEN
499  CALL dlascl( 'H', -1, -1, anrmto, anrm, n, n, a, lda, iinfo )
500  IF( iinfo.NE.0 ) THEN
501  info = n + 9
502  RETURN
503  END IF
504  CALL dlascl( 'G', -1, -1, anrmto, anrm, n, 1, alphar, n,
505  $ iinfo )
506  IF( iinfo.NE.0 ) THEN
507  info = n + 9
508  RETURN
509  END IF
510  CALL dlascl( 'G', -1, -1, anrmto, anrm, n, 1, alphai, n,
511  $ iinfo )
512  IF( iinfo.NE.0 ) THEN
513  info = n + 9
514  RETURN
515  END IF
516  END IF
517 *
518  IF( ilbscl ) THEN
519  CALL dlascl( 'U', -1, -1, bnrmto, bnrm, n, n, b, ldb, iinfo )
520  IF( iinfo.NE.0 ) THEN
521  info = n + 9
522  RETURN
523  END IF
524  CALL dlascl( 'G', -1, -1, bnrmto, bnrm, n, 1, beta, n, iinfo )
525  IF( iinfo.NE.0 ) THEN
526  info = n + 9
527  RETURN
528  END IF
529  END IF
530 *
531  10 CONTINUE
532  work( 1 ) = lwkopt
533 *
534  RETURN
535 *
536 * End of DGEGS
537 *
538  END
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
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
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
Definition: dhgeqz.f:304
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:145
subroutine dgegs(JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK, LWORK, INFO)
DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: dgegs.f:227
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
Definition: dorgqr.f:128
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD
Definition: dgghrd.f:207
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167