OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
cgegs.f
Go to the documentation of this file.
1*> \brief <b> CGEEVX 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 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.GE. IF( IINFO0 )
436 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
437.NE. IF( IINFO0 ) 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.NE. IF( IINFO0 ) 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.GE. IF( IINFO0 )
462 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
463.NE. IF( IINFO0 ) THEN
464.GT..AND..LE. IF( IINFO0 IINFON ) THEN
465 INFO = IINFO
466.GT..AND..LE. ELSE IF( IINFON IINFO2*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.NE. IF( IINFO0 ) 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.NE. IF( IINFO0 ) 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.NE. IF( IINFO0 ) THEN
498 INFO = N + 9
499 RETURN
500 END IF
501 CALL CLASCL( 'g', -1, -1, ANRMTO, ANRM, N, 1, ALPHA, N, IINFO )
502.NE. IF( IINFO0 ) 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.NE. IF( IINFO0 ) THEN
511 INFO = N + 9
512 RETURN
513 END IF
514 CALL CLASCL( 'g', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
515.NE. IF( IINFO0 ) 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)
XERBLA
Definition xerbla.f:60
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
cgeqrf
subroutine cgeqrf(m, n, a, lda, tau, work, lwork, info)
CGEQRF
Definition cgeqrf.f:146
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
cgegs
subroutine cgegs(jobvsl, jobvsr, n, a, lda, b, ldb, alpha, beta, vsl, ldvsl, vsr, ldvsr, work, lwork, rwork, info)
CGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition cgegs.f:225
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
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
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
cunmqr
subroutine cunmqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
CUNMQR
Definition cunmqr.f:168
cgghrd
subroutine cgghrd(compq, compz, n, ilo, ihi, a, lda, b, ldb, q, ldq, z, ldz, info)
CGGHRD
Definition cgghrd.f:204
cungqr
subroutine cungqr(m, n, k, a, lda, tau, work, lwork, info)
CUNGQR
Definition cungqr.f:128
max
#define max(a, b)
Definition macros.h:21