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