LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dggesx.f
Go to the documentation of this file.
1*> \brief <b> DGGESX 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*> Download DGGESX + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggesx.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggesx.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggesx.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
20* B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
21* VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
22* LIWORK, BWORK, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER JOBVSL, JOBVSR, SENSE, SORT
26* INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
27* $ SDIM
28* ..
29* .. Array Arguments ..
30* LOGICAL BWORK( * )
31* INTEGER IWORK( * )
32* DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
33* $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
34* $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
35* $ WORK( * )
36* ..
37* .. Function Arguments ..
38* LOGICAL SELCTG
39* EXTERNAL SELCTG
40* ..
41*
42*
43*> \par Purpose:
44* =============
45*>
46*> \verbatim
47*>
48*> DGGESX computes for a pair of N-by-N real nonsymmetric matrices
49*> (A,B), the generalized eigenvalues, the real Schur form (S,T), and,
50*> optionally, the left and/or right matrices of Schur vectors (VSL and
51*> VSR). This gives the generalized Schur factorization
52*>
53*> (A,B) = ( (VSL) S (VSR)**T, (VSL) T (VSR)**T )
54*>
55*> Optionally, it also orders the eigenvalues so that a selected cluster
56*> of eigenvalues appears in the leading diagonal blocks of the upper
57*> quasi-triangular matrix S and the upper triangular matrix T; computes
58*> a reciprocal condition number for the average of the selected
59*> eigenvalues (RCONDE); and computes a reciprocal condition number for
60*> the right and left deflating subspaces corresponding to the selected
61*> eigenvalues (RCONDV). The leading columns of VSL and VSR then form
62*> an orthonormal basis for the corresponding left and right eigenspaces
63*> (deflating subspaces).
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 or for both being zero.
69*>
70*> A pair of matrices (S,T) is in generalized real Schur form if T is
71*> upper triangular with non-negative diagonal and S is block upper
72*> triangular with 1-by-1 and 2-by-2 blocks. 1-by-1 blocks correspond
73*> to real generalized eigenvalues, while 2-by-2 blocks of S will be
74*> "standardized" by making the corresponding elements of T have the
75*> form:
76*> [ a 0 ]
77*> [ 0 b ]
78*>
79*> and the pair of corresponding 2-by-2 blocks in S and T will have a
80*> complex conjugate pair of generalized eigenvalues.
81*>
82*> \endverbatim
83*
84* Arguments:
85* ==========
86*
87*> \param[in] JOBVSL
88*> \verbatim
89*> JOBVSL is CHARACTER*1
90*> = 'N': do not compute the left Schur vectors;
91*> = 'V': compute the left Schur vectors.
92*> \endverbatim
93*>
94*> \param[in] JOBVSR
95*> \verbatim
96*> JOBVSR is CHARACTER*1
97*> = 'N': do not compute the right Schur vectors;
98*> = 'V': compute the right Schur vectors.
99*> \endverbatim
100*>
101*> \param[in] SORT
102*> \verbatim
103*> SORT is CHARACTER*1
104*> Specifies whether or not to order the eigenvalues on the
105*> diagonal of the generalized Schur form.
106*> = 'N': Eigenvalues are not ordered;
107*> = 'S': Eigenvalues are ordered (see SELCTG).
108*> \endverbatim
109*>
110*> \param[in] SELCTG
111*> \verbatim
112*> SELCTG is a LOGICAL FUNCTION of three DOUBLE PRECISION arguments
113*> SELCTG must be declared EXTERNAL in the calling subroutine.
114*> If SORT = 'N', SELCTG is not referenced.
115*> If SORT = 'S', SELCTG is used to select eigenvalues to sort
116*> to the top left of the Schur form.
117*> An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
118*> SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
119*> one of a complex conjugate pair of eigenvalues is selected,
120*> then both complex eigenvalues are selected.
121*> Note that a selected complex eigenvalue may no longer satisfy
122*> SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) = .TRUE. after ordering,
123*> since ordering may change the value of complex eigenvalues
124*> (especially if the eigenvalue is ill-conditioned), in this
125*> case INFO is set to N+3.
126*> \endverbatim
127*>
128*> \param[in] SENSE
129*> \verbatim
130*> SENSE is CHARACTER*1
131*> Determines which reciprocal condition numbers are computed.
132*> = 'N': None are computed;
133*> = 'E': Computed for average of selected eigenvalues only;
134*> = 'V': Computed for selected deflating subspaces only;
135*> = 'B': Computed for both.
136*> If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
137*> \endverbatim
138*>
139*> \param[in] N
140*> \verbatim
141*> N is INTEGER
142*> The order of the matrices A, B, VSL, and VSR. N >= 0.
143*> \endverbatim
144*>
145*> \param[in,out] A
146*> \verbatim
147*> A is DOUBLE PRECISION array, dimension (LDA, N)
148*> On entry, the first of the pair of matrices.
149*> On exit, A has been overwritten by its generalized Schur
150*> form S.
151*> \endverbatim
152*>
153*> \param[in] LDA
154*> \verbatim
155*> LDA is INTEGER
156*> The leading dimension of A. LDA >= max(1,N).
157*> \endverbatim
158*>
159*> \param[in,out] B
160*> \verbatim
161*> B is DOUBLE PRECISION array, dimension (LDB, N)
162*> On entry, the second of the pair of matrices.
163*> On exit, B has been overwritten by its generalized Schur
164*> form T.
165*> \endverbatim
166*>
167*> \param[in] LDB
168*> \verbatim
169*> LDB is INTEGER
170*> The leading dimension of B. LDB >= max(1,N).
171*> \endverbatim
172*>
173*> \param[out] SDIM
174*> \verbatim
175*> SDIM is INTEGER
176*> If SORT = 'N', SDIM = 0.
177*> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
178*> for which SELCTG is true. (Complex conjugate pairs for which
179*> SELCTG is true for either eigenvalue count as 2.)
180*> \endverbatim
181*>
182*> \param[out] ALPHAR
183*> \verbatim
184*> ALPHAR is DOUBLE PRECISION array, dimension (N)
185*> \endverbatim
186*>
187*> \param[out] ALPHAI
188*> \verbatim
189*> ALPHAI is DOUBLE PRECISION array, dimension (N)
190*> \endverbatim
191*>
192*> \param[out] BETA
193*> \verbatim
194*> BETA is DOUBLE PRECISION array, dimension (N)
195*> On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
196*> be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i
197*> and BETA(j),j=1,...,N are the diagonals of the complex Schur
198*> form (S,T) that would result if the 2-by-2 diagonal blocks of
199*> the real Schur form of (A,B) were further reduced to
200*> triangular form using 2-by-2 complex unitary transformations.
201*> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
202*> positive, then the j-th and (j+1)-st eigenvalues are a
203*> complex conjugate pair, with ALPHAI(j+1) negative.
204*>
205*> Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
206*> may easily over- or underflow, and BETA(j) may even be zero.
207*> Thus, the user should avoid naively computing the ratio.
208*> However, ALPHAR and ALPHAI will be always less than and
209*> usually comparable with norm(A) in magnitude, and BETA always
210*> less than and usually comparable with norm(B).
211*> \endverbatim
212*>
213*> \param[out] VSL
214*> \verbatim
215*> VSL is DOUBLE PRECISION array, dimension (LDVSL,N)
216*> If JOBVSL = 'V', VSL will contain the left Schur vectors.
217*> Not referenced if JOBVSL = 'N'.
218*> \endverbatim
219*>
220*> \param[in] LDVSL
221*> \verbatim
222*> LDVSL is INTEGER
223*> The leading dimension of the matrix VSL. LDVSL >=1, and
224*> if JOBVSL = 'V', LDVSL >= N.
225*> \endverbatim
226*>
227*> \param[out] VSR
228*> \verbatim
229*> VSR is DOUBLE PRECISION array, dimension (LDVSR,N)
230*> If JOBVSR = 'V', VSR will contain the right Schur vectors.
231*> Not referenced if JOBVSR = 'N'.
232*> \endverbatim
233*>
234*> \param[in] LDVSR
235*> \verbatim
236*> LDVSR is INTEGER
237*> The leading dimension of the matrix VSR. LDVSR >= 1, and
238*> if JOBVSR = 'V', LDVSR >= N.
239*> \endverbatim
240*>
241*> \param[out] RCONDE
242*> \verbatim
243*> RCONDE is DOUBLE PRECISION array, dimension ( 2 )
244*> If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
245*> reciprocal condition numbers for the average of the selected
246*> eigenvalues.
247*> Not referenced if SENSE = 'N' or 'V'.
248*> \endverbatim
249*>
250*> \param[out] RCONDV
251*> \verbatim
252*> RCONDV is DOUBLE PRECISION array, dimension ( 2 )
253*> If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
254*> reciprocal condition numbers for the selected deflating
255*> subspaces.
256*> Not referenced if SENSE = 'N' or 'E'.
257*> \endverbatim
258*>
259*> \param[out] WORK
260*> \verbatim
261*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
262*> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
263*> \endverbatim
264*>
265*> \param[in] LWORK
266*> \verbatim
267*> LWORK is INTEGER
268*> The dimension of the array WORK.
269*> If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
270*> LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
271*> LWORK >= max( 8*N, 6*N+16 ).
272*> Note that 2*SDIM*(N-SDIM) <= N*N/2.
273*> Note also that an error is only returned if
274*> LWORK < max( 8*N, 6*N+16), but if SENSE = 'E' or 'V' or 'B'
275*> this may not be large enough.
276*>
277*> If LWORK = -1, then a workspace query is assumed; the routine
278*> only calculates the bound on the optimal size of the WORK
279*> array and the minimum size of the IWORK array, returns these
280*> values as the first entries of the WORK and IWORK arrays, and
281*> no error message related to LWORK or LIWORK is issued by
282*> XERBLA.
283*> \endverbatim
284*>
285*> \param[out] IWORK
286*> \verbatim
287*> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
288*> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
289*> \endverbatim
290*>
291*> \param[in] LIWORK
292*> \verbatim
293*> LIWORK is INTEGER
294*> The dimension of the array IWORK.
295*> If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
296*> LIWORK >= N+6.
297*>
298*> If LIWORK = -1, then a workspace query is assumed; the
299*> routine only calculates the bound on the optimal size of the
300*> WORK array and the minimum size of the IWORK array, returns
301*> these values as the first entries of the WORK and IWORK
302*> arrays, and no error message related to LWORK or LIWORK is
303*> issued by XERBLA.
304*> \endverbatim
305*>
306*> \param[out] BWORK
307*> \verbatim
308*> BWORK is LOGICAL array, dimension (N)
309*> Not referenced if SORT = 'N'.
310*> \endverbatim
311*>
312*> \param[out] INFO
313*> \verbatim
314*> INFO is INTEGER
315*> = 0: successful exit
316*> < 0: if INFO = -i, the i-th argument had an illegal value.
317*> = 1,...,N:
318*> The QZ iteration failed. (A,B) are not in Schur
319*> form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
320*> be correct for j=INFO+1,...,N.
321*> > N: =N+1: other than QZ iteration failed in DHGEQZ
322*> =N+2: after reordering, roundoff changed values of
323*> some complex eigenvalues so that leading
324*> eigenvalues in the Generalized Schur form no
325*> longer satisfy SELCTG=.TRUE. This could also
326*> be caused due to scaling.
327*> =N+3: reordering failed in DTGSEN.
328*> \endverbatim
329*
330* Authors:
331* ========
332*
333*> \author Univ. of Tennessee
334*> \author Univ. of California Berkeley
335*> \author Univ. of Colorado Denver
336*> \author NAG Ltd.
337*
338*> \ingroup ggesx
339*
340*> \par Further Details:
341* =====================
342*>
343*> \verbatim
344*>
345*> An approximate (asymptotic) bound on the average absolute error of
346*> the selected eigenvalues is
347*>
348*> EPS * norm((A, B)) / RCONDE( 1 ).
349*>
350*> An approximate (asymptotic) bound on the maximum angular error in
351*> the computed deflating subspaces is
352*>
353*> EPS * norm((A, B)) / RCONDV( 2 ).
354*>
355*> See LAPACK User's Guide, section 4.11 for more information.
356*> \endverbatim
357*>
358* =====================================================================
359 SUBROUTINE dggesx( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A,
360 $ LDA,
361 $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
362 $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
363 $ LIWORK, BWORK, INFO )
364*
365* -- LAPACK driver routine --
366* -- LAPACK is a software package provided by Univ. of Tennessee, --
367* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
368*
369* .. Scalar Arguments ..
370 CHARACTER JOBVSL, JOBVSR, SENSE, SORT
371 INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
372 $ SDIM
373* ..
374* .. Array Arguments ..
375 LOGICAL BWORK( * )
376 INTEGER IWORK( * )
377 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
378 $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
379 $ rcondv( 2 ), vsl( ldvsl, * ), vsr( ldvsr, * ),
380 $ work( * )
381* ..
382* .. Function Arguments ..
383 LOGICAL SELCTG
384 EXTERNAL SELCTG
385* ..
386*
387* =====================================================================
388*
389* .. Parameters ..
390 DOUBLE PRECISION ZERO, ONE
391 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
392* ..
393* .. Local Scalars ..
394 LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
395 $ lquery, lst2sl, wantsb, wantse, wantsn, wantst,
396 $ wantsv
397 INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
398 $ ileft, ilo, ip, iright, irows, itau, iwrk,
399 $ liwmin, lwrk, maxwrk, minwrk
400 DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
401 $ pr, safmax, safmin, smlnum
402* ..
403* .. Local Arrays ..
404 DOUBLE PRECISION DIF( 2 )
405* ..
406* .. External Subroutines ..
407 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ,
408 $ dlacpy,
410* ..
411* .. External Functions ..
412 LOGICAL LSAME
413 INTEGER ILAENV
414 DOUBLE PRECISION DLAMCH, DLANGE
415 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
416* ..
417* .. Intrinsic Functions ..
418 INTRINSIC abs, max, sqrt
419* ..
420* .. Executable Statements ..
421*
422* Decode the input arguments
423*
424 IF( lsame( jobvsl, 'N' ) ) THEN
425 ijobvl = 1
426 ilvsl = .false.
427 ELSE IF( lsame( jobvsl, 'V' ) ) THEN
428 ijobvl = 2
429 ilvsl = .true.
430 ELSE
431 ijobvl = -1
432 ilvsl = .false.
433 END IF
434*
435 IF( lsame( jobvsr, 'N' ) ) THEN
436 ijobvr = 1
437 ilvsr = .false.
438 ELSE IF( lsame( jobvsr, 'V' ) ) THEN
439 ijobvr = 2
440 ilvsr = .true.
441 ELSE
442 ijobvr = -1
443 ilvsr = .false.
444 END IF
445*
446 wantst = lsame( sort, 'S' )
447 wantsn = lsame( sense, 'N' )
448 wantse = lsame( sense, 'E' )
449 wantsv = lsame( sense, 'V' )
450 wantsb = lsame( sense, 'B' )
451 lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
452 IF( wantsn ) THEN
453 ijob = 0
454 ELSE IF( wantse ) THEN
455 ijob = 1
456 ELSE IF( wantsv ) THEN
457 ijob = 2
458 ELSE IF( wantsb ) THEN
459 ijob = 4
460 END IF
461*
462* Test the input arguments
463*
464 info = 0
465 IF( ijobvl.LE.0 ) THEN
466 info = -1
467 ELSE IF( ijobvr.LE.0 ) THEN
468 info = -2
469 ELSE IF( ( .NOT.wantst ) .AND.
470 $ ( .NOT.lsame( sort, 'N' ) ) ) THEN
471 info = -3
472 ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
473 $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
474 info = -5
475 ELSE IF( n.LT.0 ) THEN
476 info = -6
477 ELSE IF( lda.LT.max( 1, n ) ) THEN
478 info = -8
479 ELSE IF( ldb.LT.max( 1, n ) ) THEN
480 info = -10
481 ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
482 info = -16
483 ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
484 info = -18
485 END IF
486*
487* Compute workspace
488* (Note: Comments in the code beginning "Workspace:" describe the
489* minimal amount of workspace needed at that point in the code,
490* as well as the preferred amount for good performance.
491* NB refers to the optimal block size for the immediately
492* following subroutine, as returned by ILAENV.)
493*
494 IF( info.EQ.0 ) THEN
495 IF( n.GT.0) THEN
496 minwrk = max( 8*n, 6*n + 16 )
497 maxwrk = minwrk - n +
498 $ n*ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 )
499 maxwrk = max( maxwrk, minwrk - n +
500 $ n*ilaenv( 1, 'DORMQR', ' ', n, 1, n, -1 ) )
501 IF( ilvsl ) THEN
502 maxwrk = max( maxwrk, minwrk - n +
503 $ n*ilaenv( 1, 'DORGQR', ' ', n, 1, n, -1 ) )
504 END IF
505 lwrk = maxwrk
506 IF( ijob.GE.1 )
507 $ lwrk = max( lwrk, n*n/2 )
508 ELSE
509 minwrk = 1
510 maxwrk = 1
511 lwrk = 1
512 END IF
513 work( 1 ) = lwrk
514 IF( wantsn .OR. n.EQ.0 ) THEN
515 liwmin = 1
516 ELSE
517 liwmin = n + 6
518 END IF
519 iwork( 1 ) = liwmin
520*
521 IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
522 info = -22
523 ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery ) THEN
524 info = -24
525 END IF
526 END IF
527*
528 IF( info.NE.0 ) THEN
529 CALL xerbla( 'DGGESX', -info )
530 RETURN
531 ELSE IF (lquery) THEN
532 RETURN
533 END IF
534*
535* Quick return if possible
536*
537 IF( n.EQ.0 ) THEN
538 sdim = 0
539 RETURN
540 END IF
541*
542* Get machine constants
543*
544 eps = dlamch( 'P' )
545 safmin = dlamch( 'S' )
546 safmax = one / safmin
547 smlnum = sqrt( safmin ) / eps
548 bignum = one / smlnum
549*
550* Scale A if max element outside range [SMLNUM,BIGNUM]
551*
552 anrm = dlange( 'M', n, n, a, lda, work )
553 ilascl = .false.
554 IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
555 anrmto = smlnum
556 ilascl = .true.
557 ELSE IF( anrm.GT.bignum ) THEN
558 anrmto = bignum
559 ilascl = .true.
560 END IF
561 IF( ilascl )
562 $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
563*
564* Scale B if max element outside range [SMLNUM,BIGNUM]
565*
566 bnrm = dlange( 'M', n, n, b, ldb, work )
567 ilbscl = .false.
568 IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
569 bnrmto = smlnum
570 ilbscl = .true.
571 ELSE IF( bnrm.GT.bignum ) THEN
572 bnrmto = bignum
573 ilbscl = .true.
574 END IF
575 IF( ilbscl )
576 $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
577*
578* Permute the matrix to make it more nearly triangular
579* (Workspace: need 6*N + 2*N for permutation parameters)
580*
581 ileft = 1
582 iright = n + 1
583 iwrk = iright + n
584 CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
585 $ work( iright ), work( iwrk ), ierr )
586*
587* Reduce B to triangular form (QR decomposition of B)
588* (Workspace: need N, prefer N*NB)
589*
590 irows = ihi + 1 - ilo
591 icols = n + 1 - ilo
592 itau = iwrk
593 iwrk = itau + irows
594 CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
595 $ work( iwrk ), lwork+1-iwrk, ierr )
596*
597* Apply the orthogonal transformation to matrix A
598* (Workspace: need N, prefer N*NB)
599*
600 CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
601 $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
602 $ lwork+1-iwrk, ierr )
603*
604* Initialize VSL
605* (Workspace: need N, prefer N*NB)
606*
607 IF( ilvsl ) THEN
608 CALL dlaset( 'Full', n, n, zero, one, vsl, ldvsl )
609 IF( irows.GT.1 ) THEN
610 CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
611 $ vsl( ilo+1, ilo ), ldvsl )
612 END IF
613 CALL dorgqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
614 $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
615 END IF
616*
617* Initialize VSR
618*
619 IF( ilvsr )
620 $ CALL dlaset( 'Full', n, n, zero, one, vsr, ldvsr )
621*
622* Reduce to generalized Hessenberg form
623* (Workspace: none needed)
624*
625 CALL dgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
626 $ ldvsl, vsr, ldvsr, ierr )
627*
628 sdim = 0
629*
630* Perform QZ algorithm, computing Schur vectors if desired
631* (Workspace: need N)
632*
633 iwrk = itau
634 CALL dhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
635 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
636 $ work( iwrk ), lwork+1-iwrk, ierr )
637 IF( ierr.NE.0 ) THEN
638 IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
639 info = ierr
640 ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
641 info = ierr - n
642 ELSE
643 info = n + 1
644 END IF
645 GO TO 60
646 END IF
647*
648* Sort eigenvalues ALPHA/BETA and compute the reciprocal of
649* condition number(s)
650* (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
651* otherwise, need 8*(N+1) )
652*
653 IF( wantst ) THEN
654*
655* Undo scaling on eigenvalues before SELCTGing
656*
657 IF( ilascl ) THEN
658 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
659 $ ierr )
660 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
661 $ ierr )
662 END IF
663 IF( ilbscl )
664 $ CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n,
665 $ ierr )
666*
667* Select eigenvalues
668*
669 DO 10 i = 1, n
670 bwork( i ) = selctg( alphar( i ), alphai( i ),
671 $ beta( i ) )
672 10 CONTINUE
673*
674* Reorder eigenvalues, transform Generalized Schur vectors, and
675* compute reciprocal condition numbers
676*
677 CALL dtgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
678 $ alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr,
679 $ sdim, pl, pr, dif, work( iwrk ), lwork-iwrk+1,
680 $ iwork, liwork, ierr )
681*
682 IF( ijob.GE.1 )
683 $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
684 IF( ierr.EQ.-22 ) THEN
685*
686* not enough real workspace
687*
688 info = -22
689 ELSE
690 IF( ijob.EQ.1 .OR. ijob.EQ.4 ) THEN
691 rconde( 1 ) = pl
692 rconde( 2 ) = pr
693 END IF
694 IF( ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
695 rcondv( 1 ) = dif( 1 )
696 rcondv( 2 ) = dif( 2 )
697 END IF
698 IF( ierr.EQ.1 )
699 $ info = n + 3
700 END IF
701*
702 END IF
703*
704* Apply permutation to VSL and VSR
705* (Workspace: none needed)
706*
707 IF( ilvsl )
708 $ CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
709 $ work( iright ), n, vsl, ldvsl, ierr )
710*
711 IF( ilvsr )
712 $ CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
713 $ work( iright ), n, vsr, ldvsr, ierr )
714*
715* Check if unscaling would cause over/underflow, if so, rescale
716* (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
717* B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
718*
719 IF( ilascl ) THEN
720 DO 20 i = 1, n
721 IF( alphai( i ).NE.zero ) THEN
722 IF( ( alphar( i ) / safmax ).GT.( anrmto / anrm ) .OR.
723 $ ( safmin / alphar( i ) ).GT.( anrm / anrmto ) ) THEN
724 work( 1 ) = abs( a( i, i ) / alphar( i ) )
725 beta( i ) = beta( i )*work( 1 )
726 alphar( i ) = alphar( i )*work( 1 )
727 alphai( i ) = alphai( i )*work( 1 )
728 ELSE IF( ( alphai( i ) / safmax ).GT.
729 $ ( anrmto / anrm ) .OR.
730 $ ( safmin / alphai( i ) ).GT.( anrm / anrmto ) )
731 $ THEN
732 work( 1 ) = abs( a( i, i+1 ) / alphai( i ) )
733 beta( i ) = beta( i )*work( 1 )
734 alphar( i ) = alphar( i )*work( 1 )
735 alphai( i ) = alphai( i )*work( 1 )
736 END IF
737 END IF
738 20 CONTINUE
739 END IF
740*
741 IF( ilbscl ) THEN
742 DO 30 i = 1, n
743 IF( alphai( i ).NE.zero ) THEN
744 IF( ( beta( i ) / safmax ).GT.( bnrmto / bnrm ) .OR.
745 $ ( safmin / beta( i ) ).GT.( bnrm / bnrmto ) ) THEN
746 work( 1 ) = abs( b( i, i ) / beta( i ) )
747 beta( i ) = beta( i )*work( 1 )
748 alphar( i ) = alphar( i )*work( 1 )
749 alphai( i ) = alphai( i )*work( 1 )
750 END IF
751 END IF
752 30 CONTINUE
753 END IF
754*
755* Undo scaling
756*
757 IF( ilascl ) THEN
758 CALL dlascl( 'H', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
759 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n,
760 $ ierr )
761 CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n,
762 $ ierr )
763 END IF
764*
765 IF( ilbscl ) THEN
766 CALL dlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
767 CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
768 END IF
769*
770 IF( wantst ) THEN
771*
772* Check if reordering is correct
773*
774 lastsl = .true.
775 lst2sl = .true.
776 sdim = 0
777 ip = 0
778 DO 50 i = 1, n
779 cursl = selctg( alphar( i ), alphai( i ), beta( i ) )
780 IF( alphai( i ).EQ.zero ) THEN
781 IF( cursl )
782 $ sdim = sdim + 1
783 ip = 0
784 IF( cursl .AND. .NOT.lastsl )
785 $ info = n + 2
786 ELSE
787 IF( ip.EQ.1 ) THEN
788*
789* Last eigenvalue of conjugate pair
790*
791 cursl = cursl .OR. lastsl
792 lastsl = cursl
793 IF( cursl )
794 $ sdim = sdim + 2
795 ip = -1
796 IF( cursl .AND. .NOT.lst2sl )
797 $ info = n + 2
798 ELSE
799*
800* First eigenvalue of conjugate pair
801*
802 ip = 1
803 END IF
804 END IF
805 lst2sl = lastsl
806 lastsl = cursl
807 50 CONTINUE
808*
809 END IF
810*
811 60 CONTINUE
812*
813 work( 1 ) = maxwrk
814 iwork( 1 ) = liwmin
815*
816 RETURN
817*
818* End of DGGESX
819*
820 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine dggesx(jobvsl, jobvsr, sort, selctg, sense, n, a, lda, b, ldb, sdim, alphar, alphai, beta, vsl, ldvsl, vsr, ldvsr, rconde, rcondv, work, lwork, iwork, liwork, bwork, info)
DGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE ...
Definition dggesx.f:364
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
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:142
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:108
subroutine dtgsen(ijob, wantq, wantz, select, n, a, lda, b, ldb, alphar, alphai, beta, q, ldq, z, ldz, m, pl, pr, dif, work, lwork, iwork, liwork, info)
DTGSEN
Definition dtgsen.f:450
subroutine dorgqr(m, n, k, a, lda, tau, work, lwork, info)
DORGQR
Definition dorgqr.f:126
subroutine dormqr(side, trans, m, n, k, a, lda, tau, c, ldc, work, lwork, info)
DORMQR
Definition dormqr.f:165