LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dhgeqz.f
Go to the documentation of this file.
1*> \brief \b DHGEQZ
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DHGEQZ + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dhgeqz.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dhgeqz.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dhgeqz.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DHGEQZ( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
22* ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
23* LWORK, INFO )
24*
25* .. Scalar Arguments ..
26* CHARACTER COMPQ, COMPZ, JOB
27* INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
28* ..
29* .. Array Arguments ..
30* DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ),
31* $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
32* $ WORK( * ), Z( LDZ, * )
33* ..
34*
35*
36*> \par Purpose:
37* =============
38*>
39*> \verbatim
40*>
41*> DHGEQZ computes the eigenvalues of a real matrix pair (H,T),
42*> where H is an upper Hessenberg matrix and T is upper triangular,
43*> using the double-shift QZ method.
44*> Matrix pairs of this type are produced by the reduction to
45*> generalized upper Hessenberg form of a real matrix pair (A,B):
46*>
47*> A = Q1*H*Z1**T, B = Q1*T*Z1**T,
48*>
49*> as computed by DGGHRD.
50*>
51*> If JOB='S', then the Hessenberg-triangular pair (H,T) is
52*> also reduced to generalized Schur form,
53*>
54*> H = Q*S*Z**T, T = Q*P*Z**T,
55*>
56*> where Q and Z are orthogonal matrices, P is an upper triangular
57*> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
58*> diagonal blocks.
59*>
60*> The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
61*> (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
62*> eigenvalues.
63*>
64*> Additionally, the 2-by-2 upper triangular diagonal blocks of P
65*> corresponding to 2-by-2 blocks of S are reduced to positive diagonal
66*> form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
67*> P(j,j) > 0, and P(j+1,j+1) > 0.
68*>
69*> Optionally, the orthogonal matrix Q from the generalized Schur
70*> factorization may be postmultiplied into an input matrix Q1, and the
71*> orthogonal matrix Z may be postmultiplied into an input matrix Z1.
72*> If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
73*> the matrix pair (A,B) to generalized upper Hessenberg form, then the
74*> output matrices Q1*Q and Z1*Z are the orthogonal factors from the
75*> generalized Schur factorization of (A,B):
76*>
77*> A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T.
78*>
79*> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
80*> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
81*> complex and beta real.
82*> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
83*> generalized nonsymmetric eigenvalue problem (GNEP)
84*> A*x = lambda*B*x
85*> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
86*> alternate form of the GNEP
87*> mu*A*y = B*y.
88*> Real eigenvalues can be read directly from the generalized Schur
89*> form:
90*> alpha = S(i,i), beta = P(i,i).
91*>
92*> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
93*> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
94*> pp. 241--256.
95*> \endverbatim
96*
97* Arguments:
98* ==========
99*
100*> \param[in] JOB
101*> \verbatim
102*> JOB is CHARACTER*1
103*> = 'E': Compute eigenvalues only;
104*> = 'S': Compute eigenvalues and the Schur form.
105*> \endverbatim
106*>
107*> \param[in] COMPQ
108*> \verbatim
109*> COMPQ is CHARACTER*1
110*> = 'N': Left Schur vectors (Q) are not computed;
111*> = 'I': Q is initialized to the unit matrix and the matrix Q
112*> of left Schur vectors of (H,T) is returned;
113*> = 'V': Q must contain an orthogonal matrix Q1 on entry and
114*> the product Q1*Q is returned.
115*> \endverbatim
116*>
117*> \param[in] COMPZ
118*> \verbatim
119*> COMPZ is CHARACTER*1
120*> = 'N': Right Schur vectors (Z) are not computed;
121*> = 'I': Z is initialized to the unit matrix and the matrix Z
122*> of right Schur vectors of (H,T) is returned;
123*> = 'V': Z must contain an orthogonal matrix Z1 on entry and
124*> the product Z1*Z is returned.
125*> \endverbatim
126*>
127*> \param[in] N
128*> \verbatim
129*> N is INTEGER
130*> The order of the matrices H, T, Q, and Z. N >= 0.
131*> \endverbatim
132*>
133*> \param[in] ILO
134*> \verbatim
135*> ILO is INTEGER
136*> \endverbatim
137*>
138*> \param[in] IHI
139*> \verbatim
140*> IHI is INTEGER
141*> ILO and IHI mark the rows and columns of H which are in
142*> Hessenberg form. It is assumed that A is already upper
143*> triangular in rows and columns 1:ILO-1 and IHI+1:N.
144*> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
145*> \endverbatim
146*>
147*> \param[in,out] H
148*> \verbatim
149*> H is DOUBLE PRECISION array, dimension (LDH, N)
150*> On entry, the N-by-N upper Hessenberg matrix H.
151*> On exit, if JOB = 'S', H contains the upper quasi-triangular
152*> matrix S from the generalized Schur factorization.
153*> If JOB = 'E', the diagonal blocks of H match those of S, but
154*> the rest of H is unspecified.
155*> \endverbatim
156*>
157*> \param[in] LDH
158*> \verbatim
159*> LDH is INTEGER
160*> The leading dimension of the array H. LDH >= max( 1, N ).
161*> \endverbatim
162*>
163*> \param[in,out] T
164*> \verbatim
165*> T is DOUBLE PRECISION array, dimension (LDT, N)
166*> On entry, the N-by-N upper triangular matrix T.
167*> On exit, if JOB = 'S', T contains the upper triangular
168*> matrix P from the generalized Schur factorization;
169*> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
170*> are reduced to positive diagonal form, i.e., if H(j+1,j) is
171*> non-zero, then T(j+1,j) = T(j,j+1) = 0, T(j,j) > 0, and
172*> T(j+1,j+1) > 0.
173*> If JOB = 'E', the diagonal blocks of T match those of P, but
174*> the rest of T is unspecified.
175*> \endverbatim
176*>
177*> \param[in] LDT
178*> \verbatim
179*> LDT is INTEGER
180*> The leading dimension of the array T. LDT >= max( 1, N ).
181*> \endverbatim
182*>
183*> \param[out] ALPHAR
184*> \verbatim
185*> ALPHAR is DOUBLE PRECISION array, dimension (N)
186*> The real parts of each scalar alpha defining an eigenvalue
187*> of GNEP.
188*> \endverbatim
189*>
190*> \param[out] ALPHAI
191*> \verbatim
192*> ALPHAI is DOUBLE PRECISION array, dimension (N)
193*> The imaginary parts of each scalar alpha defining an
194*> eigenvalue of GNEP.
195*> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
196*> positive, then the j-th and (j+1)-st eigenvalues are a
197*> complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
198*> \endverbatim
199*>
200*> \param[out] BETA
201*> \verbatim
202*> BETA is DOUBLE PRECISION array, dimension (N)
203*> The scalars beta that define the eigenvalues of GNEP.
204*> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
205*> beta = BETA(j) represent the j-th eigenvalue of the matrix
206*> pair (A,B), in one of the forms lambda = alpha/beta or
207*> mu = beta/alpha. Since either lambda or mu may overflow,
208*> they should not, in general, be computed.
209*> \endverbatim
210*>
211*> \param[in,out] Q
212*> \verbatim
213*> Q is DOUBLE PRECISION array, dimension (LDQ, N)
214*> On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
215*> the reduction of (A,B) to generalized Hessenberg form.
216*> On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
217*> vectors of (H,T), and if COMPQ = 'V', the orthogonal matrix
218*> of left Schur vectors of (A,B).
219*> Not referenced if COMPQ = 'N'.
220*> \endverbatim
221*>
222*> \param[in] LDQ
223*> \verbatim
224*> LDQ is INTEGER
225*> The leading dimension of the array Q. LDQ >= 1.
226*> If COMPQ='V' or 'I', then LDQ >= N.
227*> \endverbatim
228*>
229*> \param[in,out] Z
230*> \verbatim
231*> Z is DOUBLE PRECISION array, dimension (LDZ, N)
232*> On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
233*> the reduction of (A,B) to generalized Hessenberg form.
234*> On exit, if COMPZ = 'I', the orthogonal matrix of
235*> right Schur vectors of (H,T), and if COMPZ = 'V', the
236*> orthogonal matrix of right Schur vectors of (A,B).
237*> Not referenced if COMPZ = 'N'.
238*> \endverbatim
239*>
240*> \param[in] LDZ
241*> \verbatim
242*> LDZ is INTEGER
243*> The leading dimension of the array Z. LDZ >= 1.
244*> If COMPZ='V' or 'I', then LDZ >= N.
245*> \endverbatim
246*>
247*> \param[out] WORK
248*> \verbatim
249*> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
250*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
251*> \endverbatim
252*>
253*> \param[in] LWORK
254*> \verbatim
255*> LWORK is INTEGER
256*> The dimension of the array WORK. LWORK >= max(1,N).
257*>
258*> If LWORK = -1, then a workspace query is assumed; the routine
259*> only calculates the optimal size of the WORK array, returns
260*> this value as the first entry of the WORK array, and no error
261*> message related to LWORK is issued by XERBLA.
262*> \endverbatim
263*>
264*> \param[out] INFO
265*> \verbatim
266*> INFO is INTEGER
267*> = 0: successful exit
268*> < 0: if INFO = -i, the i-th argument had an illegal value
269*> = 1,...,N: the QZ iteration did not converge. (H,T) is not
270*> in Schur form, but ALPHAR(i), ALPHAI(i), and
271*> BETA(i), i=INFO+1,...,N should be correct.
272*> = N+1,...,2*N: the shift calculation failed. (H,T) is not
273*> in Schur form, but ALPHAR(i), ALPHAI(i), and
274*> BETA(i), i=INFO-N+1,...,N should be correct.
275*> \endverbatim
276*
277* Authors:
278* ========
279*
280*> \author Univ. of Tennessee
281*> \author Univ. of California Berkeley
282*> \author Univ. of Colorado Denver
283*> \author NAG Ltd.
284*
285*> \ingroup hgeqz
286*
287*> \par Further Details:
288* =====================
289*>
290*> \verbatim
291*>
292*> Iteration counters:
293*>
294*> JITER -- counts iterations.
295*> IITER -- counts iterations run since ILAST was last
296*> changed. This is therefore reset only when a 1-by-1 or
297*> 2-by-2 block deflates off the bottom.
298*> \endverbatim
299*>
300* =====================================================================
301 SUBROUTINE dhgeqz( JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT,
302 $ ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK,
303 $ LWORK, INFO )
304*
305* -- LAPACK computational routine --
306* -- LAPACK is a software package provided by Univ. of Tennessee, --
307* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
308*
309* .. Scalar Arguments ..
310 CHARACTER COMPQ, COMPZ, JOB
311 INTEGER IHI, ILO, INFO, LDH, LDQ, LDT, LDZ, LWORK, N
312* ..
313* .. Array Arguments ..
314 DOUBLE PRECISION ALPHAI( * ), ALPHAR( * ), BETA( * ),
315 $ H( LDH, * ), Q( LDQ, * ), T( LDT, * ),
316 $ work( * ), z( ldz, * )
317* ..
318*
319* =====================================================================
320*
321* .. Parameters ..
322* $ SAFETY = 1.0E+0 )
323 DOUBLE PRECISION HALF, ZERO, ONE, SAFETY
324 PARAMETER ( HALF = 0.5d+0, zero = 0.0d+0, one = 1.0d+0,
325 $ safety = 1.0d+2 )
326* ..
327* .. Local Scalars ..
328 LOGICAL ILAZR2, ILAZRO, ILPIVT, ILQ, ILSCHR, ILZ,
329 $ LQUERY
330 INTEGER ICOMPQ, ICOMPZ, IFIRST, IFRSTM, IITER, ILAST,
331 $ ILASTM, IN, ISCHUR, ISTART, J, JC, JCH, JITER,
332 $ jr, maxit
333 DOUBLE PRECISION A11, A12, A1I, A1R, A21, A22, A2I, A2R, AD11,
334 $ AD11L, AD12, AD12L, AD21, AD21L, AD22, AD22L,
335 $ ad32l, an, anorm, ascale, atol, b11, b1a, b1i,
336 $ b1r, b22, b2a, b2i, b2r, bn, bnorm, bscale,
337 $ btol, c, c11i, c11r, c12, c21, c22i, c22r, cl,
338 $ cq, cr, cz, eshift, s, s1, s1inv, s2, safmax,
339 $ safmin, scale, sl, sqi, sqr, sr, szi, szr, t1,
340 $ t2, t3, tau, temp, temp2, tempi, tempr, u1,
341 $ u12, u12l, u2, ulp, vs, w11, w12, w21, w22,
342 $ wabs, wi, wr, wr2
343* ..
344* .. Local Arrays ..
345 DOUBLE PRECISION V( 3 )
346* ..
347* .. External Functions ..
348 LOGICAL LSAME
349 DOUBLE PRECISION DLAMCH, DLANHS, DLAPY2, DLAPY3
350 EXTERNAL lsame, dlamch, dlanhs, dlapy2, dlapy3
351* ..
352* .. External Subroutines ..
353 EXTERNAL dlag2, dlarfg, dlartg, dlaset, dlasv2, drot,
354 $ xerbla
355* ..
356* .. Intrinsic Functions ..
357 INTRINSIC abs, dble, max, min, sqrt
358* ..
359* .. Executable Statements ..
360*
361* Decode JOB, COMPQ, COMPZ
362*
363 IF( lsame( job, 'E' ) ) THEN
364 ilschr = .false.
365 ischur = 1
366 ELSE IF( lsame( job, 'S' ) ) THEN
367 ilschr = .true.
368 ischur = 2
369 ELSE
370 ischur = 0
371 END IF
372*
373 IF( lsame( compq, 'N' ) ) THEN
374 ilq = .false.
375 icompq = 1
376 ELSE IF( lsame( compq, 'V' ) ) THEN
377 ilq = .true.
378 icompq = 2
379 ELSE IF( lsame( compq, 'I' ) ) THEN
380 ilq = .true.
381 icompq = 3
382 ELSE
383 icompq = 0
384 END IF
385*
386 IF( lsame( compz, 'N' ) ) THEN
387 ilz = .false.
388 icompz = 1
389 ELSE IF( lsame( compz, 'V' ) ) THEN
390 ilz = .true.
391 icompz = 2
392 ELSE IF( lsame( compz, 'I' ) ) THEN
393 ilz = .true.
394 icompz = 3
395 ELSE
396 icompz = 0
397 END IF
398*
399* Check Argument Values
400*
401 info = 0
402 work( 1 ) = max( 1, n )
403 lquery = ( lwork.EQ.-1 )
404 IF( ischur.EQ.0 ) THEN
405 info = -1
406 ELSE IF( icompq.EQ.0 ) THEN
407 info = -2
408 ELSE IF( icompz.EQ.0 ) THEN
409 info = -3
410 ELSE IF( n.LT.0 ) THEN
411 info = -4
412 ELSE IF( ilo.LT.1 ) THEN
413 info = -5
414 ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
415 info = -6
416 ELSE IF( ldh.LT.n ) THEN
417 info = -8
418 ELSE IF( ldt.LT.n ) THEN
419 info = -10
420 ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
421 info = -15
422 ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
423 info = -17
424 ELSE IF( lwork.LT.max( 1, n ) .AND. .NOT.lquery ) THEN
425 info = -19
426 END IF
427 IF( info.NE.0 ) THEN
428 CALL xerbla( 'DHGEQZ', -info )
429 RETURN
430 ELSE IF( lquery ) THEN
431 RETURN
432 END IF
433*
434* Quick return if possible
435*
436 IF( n.LE.0 ) THEN
437 work( 1 ) = dble( 1 )
438 RETURN
439 END IF
440*
441* Initialize Q and Z
442*
443 IF( icompq.EQ.3 )
444 $ CALL dlaset( 'Full', n, n, zero, one, q, ldq )
445 IF( icompz.EQ.3 )
446 $ CALL dlaset( 'Full', n, n, zero, one, z, ldz )
447*
448* Machine Constants
449*
450 in = ihi + 1 - ilo
451 safmin = dlamch( 'S' )
452 safmax = one / safmin
453 ulp = dlamch( 'E' )*dlamch( 'B' )
454 anorm = dlanhs( 'F', in, h( ilo, ilo ), ldh, work )
455 bnorm = dlanhs( 'F', in, t( ilo, ilo ), ldt, work )
456 atol = max( safmin, ulp*anorm )
457 btol = max( safmin, ulp*bnorm )
458 ascale = one / max( safmin, anorm )
459 bscale = one / max( safmin, bnorm )
460*
461* Set Eigenvalues IHI+1:N
462*
463 DO 30 j = ihi + 1, n
464 IF( t( j, j ).LT.zero ) THEN
465 IF( ilschr ) THEN
466 DO 10 jr = 1, j
467 h( jr, j ) = -h( jr, j )
468 t( jr, j ) = -t( jr, j )
469 10 CONTINUE
470 ELSE
471 h( j, j ) = -h( j, j )
472 t( j, j ) = -t( j, j )
473 END IF
474 IF( ilz ) THEN
475 DO 20 jr = 1, n
476 z( jr, j ) = -z( jr, j )
477 20 CONTINUE
478 END IF
479 END IF
480 alphar( j ) = h( j, j )
481 alphai( j ) = zero
482 beta( j ) = t( j, j )
483 30 CONTINUE
484*
485* If IHI < ILO, skip QZ steps
486*
487 IF( ihi.LT.ilo )
488 $ GO TO 380
489*
490* MAIN QZ ITERATION LOOP
491*
492* Initialize dynamic indices
493*
494* Eigenvalues ILAST+1:N have been found.
495* Column operations modify rows IFRSTM:whatever.
496* Row operations modify columns whatever:ILASTM.
497*
498* If only eigenvalues are being computed, then
499* IFRSTM is the row of the last splitting row above row ILAST;
500* this is always at least ILO.
501* IITER counts iterations since the last eigenvalue was found,
502* to tell when to use an extraordinary shift.
503* MAXIT is the maximum number of QZ sweeps allowed.
504*
505 ilast = ihi
506 IF( ilschr ) THEN
507 ifrstm = 1
508 ilastm = n
509 ELSE
510 ifrstm = ilo
511 ilastm = ihi
512 END IF
513 iiter = 0
514 eshift = zero
515 maxit = 30*( ihi-ilo+1 )
516*
517 DO 360 jiter = 1, maxit
518*
519* Split the matrix if possible.
520*
521* Two tests:
522* 1: H(j,j-1)=0 or j=ILO
523* 2: T(j,j)=0
524*
525 IF( ilast.EQ.ilo ) THEN
526*
527* Special case: j=ILAST
528*
529 GO TO 80
530 ELSE
531 IF( abs( h( ilast, ilast-1 ) ).LE.max( safmin, ulp*(
532 $ abs( h( ilast, ilast ) ) + abs( h( ilast-1, ilast-1 ) )
533 $ ) ) ) THEN
534 h( ilast, ilast-1 ) = zero
535 GO TO 80
536 END IF
537 END IF
538*
539 IF( abs( t( ilast, ilast ) ).LE.btol ) THEN
540 t( ilast, ilast ) = zero
541 GO TO 70
542 END IF
543*
544* General case: j<ILAST
545*
546 DO 60 j = ilast - 1, ilo, -1
547*
548* Test 1: for H(j,j-1)=0 or j=ILO
549*
550 IF( j.EQ.ilo ) THEN
551 ilazro = .true.
552 ELSE
553 IF( abs( h( j, j-1 ) ).LE.max( safmin, ulp*(
554 $ abs( h( j, j ) ) + abs( h( j-1, j-1 ) )
555 $ ) ) ) THEN
556 h( j, j-1 ) = zero
557 ilazro = .true.
558 ELSE
559 ilazro = .false.
560 END IF
561 END IF
562*
563* Test 2: for T(j,j)=0
564*
565 IF( abs( t( j, j ) ).LT.btol ) THEN
566 t( j, j ) = zero
567*
568* Test 1a: Check for 2 consecutive small subdiagonals in A
569*
570 ilazr2 = .false.
571 IF( .NOT.ilazro ) THEN
572 temp = abs( h( j, j-1 ) )
573 temp2 = abs( h( j, j ) )
574 tempr = max( temp, temp2 )
575 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
576 temp = temp / tempr
577 temp2 = temp2 / tempr
578 END IF
579 IF( temp*( ascale*abs( h( j+1, j ) ) ).LE.temp2*
580 $ ( ascale*atol ) )ilazr2 = .true.
581 END IF
582*
583* If both tests pass (1 & 2), i.e., the leading diagonal
584* element of B in the block is zero, split a 1x1 block off
585* at the top. (I.e., at the J-th row/column) The leading
586* diagonal element of the remainder can also be zero, so
587* this may have to be done repeatedly.
588*
589 IF( ilazro .OR. ilazr2 ) THEN
590 DO 40 jch = j, ilast - 1
591 temp = h( jch, jch )
592 CALL dlartg( temp, h( jch+1, jch ), c, s,
593 $ h( jch, jch ) )
594 h( jch+1, jch ) = zero
595 CALL drot( ilastm-jch, h( jch, jch+1 ), ldh,
596 $ h( jch+1, jch+1 ), ldh, c, s )
597 CALL drot( ilastm-jch, t( jch, jch+1 ), ldt,
598 $ t( jch+1, jch+1 ), ldt, c, s )
599 IF( ilq )
600 $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
601 $ c, s )
602 IF( ilazr2 )
603 $ h( jch, jch-1 ) = h( jch, jch-1 )*c
604 ilazr2 = .false.
605 IF( abs( t( jch+1, jch+1 ) ).GE.btol ) THEN
606 IF( jch+1.GE.ilast ) THEN
607 GO TO 80
608 ELSE
609 ifirst = jch + 1
610 GO TO 110
611 END IF
612 END IF
613 t( jch+1, jch+1 ) = zero
614 40 CONTINUE
615 GO TO 70
616 ELSE
617*
618* Only test 2 passed -- chase the zero to T(ILAST,ILAST)
619* Then process as in the case T(ILAST,ILAST)=0
620*
621 DO 50 jch = j, ilast - 1
622 temp = t( jch, jch+1 )
623 CALL dlartg( temp, t( jch+1, jch+1 ), c, s,
624 $ t( jch, jch+1 ) )
625 t( jch+1, jch+1 ) = zero
626 IF( jch.LT.ilastm-1 )
627 $ CALL drot( ilastm-jch-1, t( jch, jch+2 ), ldt,
628 $ t( jch+1, jch+2 ), ldt, c, s )
629 CALL drot( ilastm-jch+2, h( jch, jch-1 ), ldh,
630 $ h( jch+1, jch-1 ), ldh, c, s )
631 IF( ilq )
632 $ CALL drot( n, q( 1, jch ), 1, q( 1, jch+1 ), 1,
633 $ c, s )
634 temp = h( jch+1, jch )
635 CALL dlartg( temp, h( jch+1, jch-1 ), c, s,
636 $ h( jch+1, jch ) )
637 h( jch+1, jch-1 ) = zero
638 CALL drot( jch+1-ifrstm, h( ifrstm, jch ), 1,
639 $ h( ifrstm, jch-1 ), 1, c, s )
640 CALL drot( jch-ifrstm, t( ifrstm, jch ), 1,
641 $ t( ifrstm, jch-1 ), 1, c, s )
642 IF( ilz )
643 $ CALL drot( n, z( 1, jch ), 1, z( 1, jch-1 ), 1,
644 $ c, s )
645 50 CONTINUE
646 GO TO 70
647 END IF
648 ELSE IF( ilazro ) THEN
649*
650* Only test 1 passed -- work on J:ILAST
651*
652 ifirst = j
653 GO TO 110
654 END IF
655*
656* Neither test passed -- try next J
657*
658 60 CONTINUE
659*
660* (Drop-through is "impossible")
661*
662 info = n + 1
663 GO TO 420
664*
665* T(ILAST,ILAST)=0 -- clear H(ILAST,ILAST-1) to split off a
666* 1x1 block.
667*
668 70 CONTINUE
669 temp = h( ilast, ilast )
670 CALL dlartg( temp, h( ilast, ilast-1 ), c, s,
671 $ h( ilast, ilast ) )
672 h( ilast, ilast-1 ) = zero
673 CALL drot( ilast-ifrstm, h( ifrstm, ilast ), 1,
674 $ h( ifrstm, ilast-1 ), 1, c, s )
675 CALL drot( ilast-ifrstm, t( ifrstm, ilast ), 1,
676 $ t( ifrstm, ilast-1 ), 1, c, s )
677 IF( ilz )
678 $ CALL drot( n, z( 1, ilast ), 1, z( 1, ilast-1 ), 1, c, s )
679*
680* H(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
681* and BETA
682*
683 80 CONTINUE
684 IF( t( ilast, ilast ).LT.zero ) THEN
685 IF( ilschr ) THEN
686 DO 90 j = ifrstm, ilast
687 h( j, ilast ) = -h( j, ilast )
688 t( j, ilast ) = -t( j, ilast )
689 90 CONTINUE
690 ELSE
691 h( ilast, ilast ) = -h( ilast, ilast )
692 t( ilast, ilast ) = -t( ilast, ilast )
693 END IF
694 IF( ilz ) THEN
695 DO 100 j = 1, n
696 z( j, ilast ) = -z( j, ilast )
697 100 CONTINUE
698 END IF
699 END IF
700 alphar( ilast ) = h( ilast, ilast )
701 alphai( ilast ) = zero
702 beta( ilast ) = t( ilast, ilast )
703*
704* Go to next block -- exit if finished.
705*
706 ilast = ilast - 1
707 IF( ilast.LT.ilo )
708 $ GO TO 380
709*
710* Reset counters
711*
712 iiter = 0
713 eshift = zero
714 IF( .NOT.ilschr ) THEN
715 ilastm = ilast
716 IF( ifrstm.GT.ilast )
717 $ ifrstm = ilo
718 END IF
719 GO TO 350
720*
721* QZ step
722*
723* This iteration only involves rows/columns IFIRST:ILAST. We
724* assume IFIRST < ILAST, and that the diagonal of B is non-zero.
725*
726 110 CONTINUE
727 iiter = iiter + 1
728 IF( .NOT.ilschr ) THEN
729 ifrstm = ifirst
730 END IF
731*
732* Compute single shifts.
733*
734* At this point, IFIRST < ILAST, and the diagonal elements of
735* T(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
736* magnitude)
737*
738 IF( ( iiter / 10 )*10.EQ.iiter ) THEN
739*
740* Exceptional shift. Chosen for no particularly good reason.
741* (Single shift only.)
742*
743 IF( ( dble( maxit )*safmin )*abs( h( ilast, ilast-1 ) ).LT.
744 $ abs( t( ilast-1, ilast-1 ) ) ) THEN
745 eshift = h( ilast, ilast-1 ) /
746 $ t( ilast-1, ilast-1 )
747 ELSE
748 eshift = eshift + one / ( safmin*dble( maxit ) )
749 END IF
750 s1 = one
751 wr = eshift
752*
753 ELSE
754*
755* Shifts based on the generalized eigenvalues of the
756* bottom-right 2x2 block of A and B. The first eigenvalue
757* returned by DLAG2 is the Wilkinson shift (AEP p.512),
758*
759 CALL dlag2( h( ilast-1, ilast-1 ), ldh,
760 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
761 $ s2, wr, wr2, wi )
762*
763 IF ( abs( (wr/s1)*t( ilast, ilast ) - h( ilast, ilast ) )
764 $ .GT. abs( (wr2/s2)*t( ilast, ilast )
765 $ - h( ilast, ilast ) ) ) THEN
766 temp = wr
767 wr = wr2
768 wr2 = temp
769 temp = s1
770 s1 = s2
771 s2 = temp
772 END IF
773 temp = max( s1, safmin*max( one, abs( wr ), abs( wi ) ) )
774 IF( wi.NE.zero )
775 $ GO TO 200
776 END IF
777*
778* Fiddle with shift to avoid overflow
779*
780 temp = min( ascale, one )*( half*safmax )
781 IF( s1.GT.temp ) THEN
782 scale = temp / s1
783 ELSE
784 scale = one
785 END IF
786*
787 temp = min( bscale, one )*( half*safmax )
788 IF( abs( wr ).GT.temp )
789 $ scale = min( scale, temp / abs( wr ) )
790 s1 = scale*s1
791 wr = scale*wr
792*
793* Now check for two consecutive small subdiagonals.
794*
795 DO 120 j = ilast - 1, ifirst + 1, -1
796 istart = j
797 temp = abs( s1*h( j, j-1 ) )
798 temp2 = abs( s1*h( j, j )-wr*t( j, j ) )
799 tempr = max( temp, temp2 )
800 IF( tempr.LT.one .AND. tempr.NE.zero ) THEN
801 temp = temp / tempr
802 temp2 = temp2 / tempr
803 END IF
804 IF( abs( ( ascale*h( j+1, j ) )*temp ).LE.( ascale*atol )*
805 $ temp2 )GO TO 130
806 120 CONTINUE
807*
808 istart = ifirst
809 130 CONTINUE
810*
811* Do an implicit single-shift QZ sweep.
812*
813* Initial Q
814*
815 temp = s1*h( istart, istart ) - wr*t( istart, istart )
816 temp2 = s1*h( istart+1, istart )
817 CALL dlartg( temp, temp2, c, s, tempr )
818*
819* Sweep
820*
821 DO 190 j = istart, ilast - 1
822 IF( j.GT.istart ) THEN
823 temp = h( j, j-1 )
824 CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
825 h( j+1, j-1 ) = zero
826 END IF
827*
828 DO 140 jc = j, ilastm
829 temp = c*h( j, jc ) + s*h( j+1, jc )
830 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
831 h( j, jc ) = temp
832 temp2 = c*t( j, jc ) + s*t( j+1, jc )
833 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
834 t( j, jc ) = temp2
835 140 CONTINUE
836 IF( ilq ) THEN
837 DO 150 jr = 1, n
838 temp = c*q( jr, j ) + s*q( jr, j+1 )
839 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
840 q( jr, j ) = temp
841 150 CONTINUE
842 END IF
843*
844 temp = t( j+1, j+1 )
845 CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
846 t( j+1, j ) = zero
847*
848 DO 160 jr = ifrstm, min( j+2, ilast )
849 temp = c*h( jr, j+1 ) + s*h( jr, j )
850 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
851 h( jr, j+1 ) = temp
852 160 CONTINUE
853 DO 170 jr = ifrstm, j
854 temp = c*t( jr, j+1 ) + s*t( jr, j )
855 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
856 t( jr, j+1 ) = temp
857 170 CONTINUE
858 IF( ilz ) THEN
859 DO 180 jr = 1, n
860 temp = c*z( jr, j+1 ) + s*z( jr, j )
861 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
862 z( jr, j+1 ) = temp
863 180 CONTINUE
864 END IF
865 190 CONTINUE
866*
867 GO TO 350
868*
869* Use Francis double-shift
870*
871* Note: the Francis double-shift should work with real shifts,
872* but only if the block is at least 3x3.
873* This code may break if this point is reached with
874* a 2x2 block with real eigenvalues.
875*
876 200 CONTINUE
877 IF( ifirst+1.EQ.ilast ) THEN
878*
879* Special case -- 2x2 block with complex eigenvectors
880*
881* Step 1: Standardize, that is, rotate so that
882*
883* ( B11 0 )
884* B = ( ) with B11 non-negative.
885* ( 0 B22 )
886*
887 CALL dlasv2( t( ilast-1, ilast-1 ), t( ilast-1, ilast ),
888 $ t( ilast, ilast ), b22, b11, sr, cr, sl, cl )
889*
890 IF( b11.LT.zero ) THEN
891 cr = -cr
892 sr = -sr
893 b11 = -b11
894 b22 = -b22
895 END IF
896*
897 CALL drot( ilastm+1-ifirst, h( ilast-1, ilast-1 ), ldh,
898 $ h( ilast, ilast-1 ), ldh, cl, sl )
899 CALL drot( ilast+1-ifrstm, h( ifrstm, ilast-1 ), 1,
900 $ h( ifrstm, ilast ), 1, cr, sr )
901*
902 IF( ilast.LT.ilastm )
903 $ CALL drot( ilastm-ilast, t( ilast-1, ilast+1 ), ldt,
904 $ t( ilast, ilast+1 ), ldt, cl, sl )
905 IF( ifrstm.LT.ilast-1 )
906 $ CALL drot( ifirst-ifrstm, t( ifrstm, ilast-1 ), 1,
907 $ t( ifrstm, ilast ), 1, cr, sr )
908*
909 IF( ilq )
910 $ CALL drot( n, q( 1, ilast-1 ), 1, q( 1, ilast ), 1, cl,
911 $ sl )
912 IF( ilz )
913 $ CALL drot( n, z( 1, ilast-1 ), 1, z( 1, ilast ), 1, cr,
914 $ sr )
915*
916 t( ilast-1, ilast-1 ) = b11
917 t( ilast-1, ilast ) = zero
918 t( ilast, ilast-1 ) = zero
919 t( ilast, ilast ) = b22
920*
921* If B22 is negative, negate column ILAST
922*
923 IF( b22.LT.zero ) THEN
924 DO 210 j = ifrstm, ilast
925 h( j, ilast ) = -h( j, ilast )
926 t( j, ilast ) = -t( j, ilast )
927 210 CONTINUE
928*
929 IF( ilz ) THEN
930 DO 220 j = 1, n
931 z( j, ilast ) = -z( j, ilast )
932 220 CONTINUE
933 END IF
934 b22 = -b22
935 END IF
936*
937* Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
938*
939* Recompute shift
940*
941 CALL dlag2( h( ilast-1, ilast-1 ), ldh,
942 $ t( ilast-1, ilast-1 ), ldt, safmin*safety, s1,
943 $ temp, wr, temp2, wi )
944*
945* If standardization has perturbed the shift onto real line,
946* do another (real single-shift) QR step.
947*
948 IF( wi.EQ.zero )
949 $ GO TO 350
950 s1inv = one / s1
951*
952* Do EISPACK (QZVAL) computation of alpha and beta
953*
954 a11 = h( ilast-1, ilast-1 )
955 a21 = h( ilast, ilast-1 )
956 a12 = h( ilast-1, ilast )
957 a22 = h( ilast, ilast )
958*
959* Compute complex Givens rotation on right
960* (Assume some element of C = (sA - wB) > unfl )
961* __
962* (sA - wB) ( CZ -SZ )
963* ( SZ CZ )
964*
965 c11r = s1*a11 - wr*b11
966 c11i = -wi*b11
967 c12 = s1*a12
968 c21 = s1*a21
969 c22r = s1*a22 - wr*b22
970 c22i = -wi*b22
971*
972 IF( abs( c11r )+abs( c11i )+abs( c12 ).GT.abs( c21 )+
973 $ abs( c22r )+abs( c22i ) ) THEN
974 t1 = dlapy3( c12, c11r, c11i )
975 cz = c12 / t1
976 szr = -c11r / t1
977 szi = -c11i / t1
978 ELSE
979 cz = dlapy2( c22r, c22i )
980 IF( cz.LE.safmin ) THEN
981 cz = zero
982 szr = one
983 szi = zero
984 ELSE
985 tempr = c22r / cz
986 tempi = c22i / cz
987 t1 = dlapy2( cz, c21 )
988 cz = cz / t1
989 szr = -c21*tempr / t1
990 szi = c21*tempi / t1
991 END IF
992 END IF
993*
994* Compute Givens rotation on left
995*
996* ( CQ SQ )
997* ( __ ) A or B
998* ( -SQ CQ )
999*
1000 an = abs( a11 ) + abs( a12 ) + abs( a21 ) + abs( a22 )
1001 bn = abs( b11 ) + abs( b22 )
1002 wabs = abs( wr ) + abs( wi )
1003 IF( s1*an.GT.wabs*bn ) THEN
1004 cq = cz*b11
1005 sqr = szr*b22
1006 sqi = -szi*b22
1007 ELSE
1008 a1r = cz*a11 + szr*a12
1009 a1i = szi*a12
1010 a2r = cz*a21 + szr*a22
1011 a2i = szi*a22
1012 cq = dlapy2( a1r, a1i )
1013 IF( cq.LE.safmin ) THEN
1014 cq = zero
1015 sqr = one
1016 sqi = zero
1017 ELSE
1018 tempr = a1r / cq
1019 tempi = a1i / cq
1020 sqr = tempr*a2r + tempi*a2i
1021 sqi = tempi*a2r - tempr*a2i
1022 END IF
1023 END IF
1024 t1 = dlapy3( cq, sqr, sqi )
1025 cq = cq / t1
1026 sqr = sqr / t1
1027 sqi = sqi / t1
1028*
1029* Compute diagonal elements of QBZ
1030*
1031 tempr = sqr*szr - sqi*szi
1032 tempi = sqr*szi + sqi*szr
1033 b1r = cq*cz*b11 + tempr*b22
1034 b1i = tempi*b22
1035 b1a = dlapy2( b1r, b1i )
1036 b2r = cq*cz*b22 + tempr*b11
1037 b2i = -tempi*b11
1038 b2a = dlapy2( b2r, b2i )
1039*
1040* Normalize so beta > 0, and Im( alpha1 ) > 0
1041*
1042 beta( ilast-1 ) = b1a
1043 beta( ilast ) = b2a
1044 alphar( ilast-1 ) = ( wr*b1a )*s1inv
1045 alphai( ilast-1 ) = ( wi*b1a )*s1inv
1046 alphar( ilast ) = ( wr*b2a )*s1inv
1047 alphai( ilast ) = -( wi*b2a )*s1inv
1048*
1049* Step 3: Go to next block -- exit if finished.
1050*
1051 ilast = ifirst - 1
1052 IF( ilast.LT.ilo )
1053 $ GO TO 380
1054*
1055* Reset counters
1056*
1057 iiter = 0
1058 eshift = zero
1059 IF( .NOT.ilschr ) THEN
1060 ilastm = ilast
1061 IF( ifrstm.GT.ilast )
1062 $ ifrstm = ilo
1063 END IF
1064 GO TO 350
1065 ELSE
1066*
1067* Usual case: 3x3 or larger block, using Francis implicit
1068* double-shift
1069*
1070* 2
1071* Eigenvalue equation is w - c w + d = 0,
1072*
1073* -1 2 -1
1074* so compute 1st column of (A B ) - c A B + d
1075* using the formula in QZIT (from EISPACK)
1076*
1077* We assume that the block is at least 3x3
1078*
1079 ad11 = ( ascale*h( ilast-1, ilast-1 ) ) /
1080 $ ( bscale*t( ilast-1, ilast-1 ) )
1081 ad21 = ( ascale*h( ilast, ilast-1 ) ) /
1082 $ ( bscale*t( ilast-1, ilast-1 ) )
1083 ad12 = ( ascale*h( ilast-1, ilast ) ) /
1084 $ ( bscale*t( ilast, ilast ) )
1085 ad22 = ( ascale*h( ilast, ilast ) ) /
1086 $ ( bscale*t( ilast, ilast ) )
1087 u12 = t( ilast-1, ilast ) / t( ilast, ilast )
1088 ad11l = ( ascale*h( ifirst, ifirst ) ) /
1089 $ ( bscale*t( ifirst, ifirst ) )
1090 ad21l = ( ascale*h( ifirst+1, ifirst ) ) /
1091 $ ( bscale*t( ifirst, ifirst ) )
1092 ad12l = ( ascale*h( ifirst, ifirst+1 ) ) /
1093 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1094 ad22l = ( ascale*h( ifirst+1, ifirst+1 ) ) /
1095 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1096 ad32l = ( ascale*h( ifirst+2, ifirst+1 ) ) /
1097 $ ( bscale*t( ifirst+1, ifirst+1 ) )
1098 u12l = t( ifirst, ifirst+1 ) / t( ifirst+1, ifirst+1 )
1099*
1100 v( 1 ) = ( ad11-ad11l )*( ad22-ad11l ) - ad12*ad21 +
1101 $ ad21*u12*ad11l + ( ad12l-ad11l*u12l )*ad21l
1102 v( 2 ) = ( ( ad22l-ad11l )-ad21l*u12l-( ad11-ad11l )-
1103 $ ( ad22-ad11l )+ad21*u12 )*ad21l
1104 v( 3 ) = ad32l*ad21l
1105*
1106 istart = ifirst
1107*
1108 CALL dlarfg( 3, v( 1 ), v( 2 ), 1, tau )
1109 v( 1 ) = one
1110*
1111* Sweep
1112*
1113 DO 290 j = istart, ilast - 2
1114*
1115* All but last elements: use 3x3 Householder transforms.
1116*
1117* Zero (j-1)st column of A
1118*
1119 IF( j.GT.istart ) THEN
1120 v( 1 ) = h( j, j-1 )
1121 v( 2 ) = h( j+1, j-1 )
1122 v( 3 ) = h( j+2, j-1 )
1123*
1124 CALL dlarfg( 3, h( j, j-1 ), v( 2 ), 1, tau )
1125 v( 1 ) = one
1126 h( j+1, j-1 ) = zero
1127 h( j+2, j-1 ) = zero
1128 END IF
1129*
1130 t2 = tau*v( 2 )
1131 t3 = tau*v( 3 )
1132 DO 230 jc = j, ilastm
1133 temp = h( j, jc )+v( 2 )*h( j+1, jc )+v( 3 )*
1134 $ h( j+2, jc )
1135 h( j, jc ) = h( j, jc ) - temp*tau
1136 h( j+1, jc ) = h( j+1, jc ) - temp*t2
1137 h( j+2, jc ) = h( j+2, jc ) - temp*t3
1138 temp2 = t( j, jc )+v( 2 )*t( j+1, jc )+v( 3 )*
1139 $ t( j+2, jc )
1140 t( j, jc ) = t( j, jc ) - temp2*tau
1141 t( j+1, jc ) = t( j+1, jc ) - temp2*t2
1142 t( j+2, jc ) = t( j+2, jc ) - temp2*t3
1143 230 CONTINUE
1144 IF( ilq ) THEN
1145 DO 240 jr = 1, n
1146 temp = q( jr, j )+v( 2 )*q( jr, j+1 )+v( 3 )*
1147 $ q( jr, j+2 )
1148 q( jr, j ) = q( jr, j ) - temp*tau
1149 q( jr, j+1 ) = q( jr, j+1 ) - temp*t2
1150 q( jr, j+2 ) = q( jr, j+2 ) - temp*t3
1151 240 CONTINUE
1152 END IF
1153*
1154* Zero j-th column of B (see DLAGBC for details)
1155*
1156* Swap rows to pivot
1157*
1158 ilpivt = .false.
1159 temp = max( abs( t( j+1, j+1 ) ), abs( t( j+1, j+2 ) ) )
1160 temp2 = max( abs( t( j+2, j+1 ) ), abs( t( j+2, j+2 ) ) )
1161 IF( max( temp, temp2 ).LT.safmin ) THEN
1162 scale = zero
1163 u1 = one
1164 u2 = zero
1165 GO TO 250
1166 ELSE IF( temp.GE.temp2 ) THEN
1167 w11 = t( j+1, j+1 )
1168 w21 = t( j+2, j+1 )
1169 w12 = t( j+1, j+2 )
1170 w22 = t( j+2, j+2 )
1171 u1 = t( j+1, j )
1172 u2 = t( j+2, j )
1173 ELSE
1174 w21 = t( j+1, j+1 )
1175 w11 = t( j+2, j+1 )
1176 w22 = t( j+1, j+2 )
1177 w12 = t( j+2, j+2 )
1178 u2 = t( j+1, j )
1179 u1 = t( j+2, j )
1180 END IF
1181*
1182* Swap columns if nec.
1183*
1184 IF( abs( w12 ).GT.abs( w11 ) ) THEN
1185 ilpivt = .true.
1186 temp = w12
1187 temp2 = w22
1188 w12 = w11
1189 w22 = w21
1190 w11 = temp
1191 w21 = temp2
1192 END IF
1193*
1194* LU-factor
1195*
1196 temp = w21 / w11
1197 u2 = u2 - temp*u1
1198 w22 = w22 - temp*w12
1199 w21 = zero
1200*
1201* Compute SCALE
1202*
1203 scale = one
1204 IF( abs( w22 ).LT.safmin ) THEN
1205 scale = zero
1206 u2 = one
1207 u1 = -w12 / w11
1208 GO TO 250
1209 END IF
1210 IF( abs( w22 ).LT.abs( u2 ) )
1211 $ scale = abs( w22 / u2 )
1212 IF( abs( w11 ).LT.abs( u1 ) )
1213 $ scale = min( scale, abs( w11 / u1 ) )
1214*
1215* Solve
1216*
1217 u2 = ( scale*u2 ) / w22
1218 u1 = ( scale*u1-w12*u2 ) / w11
1219*
1220 250 CONTINUE
1221 IF( ilpivt ) THEN
1222 temp = u2
1223 u2 = u1
1224 u1 = temp
1225 END IF
1226*
1227* Compute Householder Vector
1228*
1229 t1 = sqrt( scale**2+u1**2+u2**2 )
1230 tau = one + scale / t1
1231 vs = -one / ( scale+t1 )
1232 v( 1 ) = one
1233 v( 2 ) = vs*u1
1234 v( 3 ) = vs*u2
1235*
1236* Apply transformations from the right.
1237*
1238 t2 = tau*v(2)
1239 t3 = tau*v(3)
1240 DO 260 jr = ifrstm, min( j+3, ilast )
1241 temp = h( jr, j )+v( 2 )*h( jr, j+1 )+v( 3 )*
1242 $ h( jr, j+2 )
1243 h( jr, j ) = h( jr, j ) - temp*tau
1244 h( jr, j+1 ) = h( jr, j+1 ) - temp*t2
1245 h( jr, j+2 ) = h( jr, j+2 ) - temp*t3
1246 260 CONTINUE
1247 DO 270 jr = ifrstm, j + 2
1248 temp = t( jr, j )+v( 2 )*t( jr, j+1 )+v( 3 )*
1249 $ t( jr, j+2 )
1250 t( jr, j ) = t( jr, j ) - temp*tau
1251 t( jr, j+1 ) = t( jr, j+1 ) - temp*t2
1252 t( jr, j+2 ) = t( jr, j+2 ) - temp*t3
1253 270 CONTINUE
1254 IF( ilz ) THEN
1255 DO 280 jr = 1, n
1256 temp = z( jr, j )+v( 2 )*z( jr, j+1 )+v( 3 )*
1257 $ z( jr, j+2 )
1258 z( jr, j ) = z( jr, j ) - temp*tau
1259 z( jr, j+1 ) = z( jr, j+1 ) - temp*t2
1260 z( jr, j+2 ) = z( jr, j+2 ) - temp*t3
1261 280 CONTINUE
1262 END IF
1263 t( j+1, j ) = zero
1264 t( j+2, j ) = zero
1265 290 CONTINUE
1266*
1267* Last elements: Use Givens rotations
1268*
1269* Rotations from the left
1270*
1271 j = ilast - 1
1272 temp = h( j, j-1 )
1273 CALL dlartg( temp, h( j+1, j-1 ), c, s, h( j, j-1 ) )
1274 h( j+1, j-1 ) = zero
1275*
1276 DO 300 jc = j, ilastm
1277 temp = c*h( j, jc ) + s*h( j+1, jc )
1278 h( j+1, jc ) = -s*h( j, jc ) + c*h( j+1, jc )
1279 h( j, jc ) = temp
1280 temp2 = c*t( j, jc ) + s*t( j+1, jc )
1281 t( j+1, jc ) = -s*t( j, jc ) + c*t( j+1, jc )
1282 t( j, jc ) = temp2
1283 300 CONTINUE
1284 IF( ilq ) THEN
1285 DO 310 jr = 1, n
1286 temp = c*q( jr, j ) + s*q( jr, j+1 )
1287 q( jr, j+1 ) = -s*q( jr, j ) + c*q( jr, j+1 )
1288 q( jr, j ) = temp
1289 310 CONTINUE
1290 END IF
1291*
1292* Rotations from the right.
1293*
1294 temp = t( j+1, j+1 )
1295 CALL dlartg( temp, t( j+1, j ), c, s, t( j+1, j+1 ) )
1296 t( j+1, j ) = zero
1297*
1298 DO 320 jr = ifrstm, ilast
1299 temp = c*h( jr, j+1 ) + s*h( jr, j )
1300 h( jr, j ) = -s*h( jr, j+1 ) + c*h( jr, j )
1301 h( jr, j+1 ) = temp
1302 320 CONTINUE
1303 DO 330 jr = ifrstm, ilast - 1
1304 temp = c*t( jr, j+1 ) + s*t( jr, j )
1305 t( jr, j ) = -s*t( jr, j+1 ) + c*t( jr, j )
1306 t( jr, j+1 ) = temp
1307 330 CONTINUE
1308 IF( ilz ) THEN
1309 DO 340 jr = 1, n
1310 temp = c*z( jr, j+1 ) + s*z( jr, j )
1311 z( jr, j ) = -s*z( jr, j+1 ) + c*z( jr, j )
1312 z( jr, j+1 ) = temp
1313 340 CONTINUE
1314 END IF
1315*
1316* End of Double-Shift code
1317*
1318 END IF
1319*
1320 GO TO 350
1321*
1322* End of iteration loop
1323*
1324 350 CONTINUE
1325 360 CONTINUE
1326*
1327* Drop-through = non-convergence
1328*
1329 info = ilast
1330 GO TO 420
1331*
1332* Successful completion of all QZ steps
1333*
1334 380 CONTINUE
1335*
1336* Set Eigenvalues 1:ILO-1
1337*
1338 DO 410 j = 1, ilo - 1
1339 IF( t( j, j ).LT.zero ) THEN
1340 IF( ilschr ) THEN
1341 DO 390 jr = 1, j
1342 h( jr, j ) = -h( jr, j )
1343 t( jr, j ) = -t( jr, j )
1344 390 CONTINUE
1345 ELSE
1346 h( j, j ) = -h( j, j )
1347 t( j, j ) = -t( j, j )
1348 END IF
1349 IF( ilz ) THEN
1350 DO 400 jr = 1, n
1351 z( jr, j ) = -z( jr, j )
1352 400 CONTINUE
1353 END IF
1354 END IF
1355 alphar( j ) = h( j, j )
1356 alphai( j ) = zero
1357 beta( j ) = t( j, j )
1358 410 CONTINUE
1359*
1360* Normal Termination
1361*
1362 info = 0
1363*
1364* Exit (other than argument error) -- return optimal workspace size
1365*
1366 420 CONTINUE
1367 work( 1 ) = dble( n )
1368 RETURN
1369*
1370* End of DHGEQZ
1371*
1372 END
subroutine xerbla(srname, info)
Definition cblat2.f:3285
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 dlag2(a, lda, b, ldb, safmin, scale1, scale2, wr1, wr2, wi)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition dlag2.f:156
subroutine dlarfg(n, alpha, x, incx, tau)
DLARFG generates an elementary reflector (Householder matrix).
Definition dlarfg.f:106
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition dlartg.f90:111
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 dlasv2(f, g, h, ssmin, ssmax, snr, csr, snl, csl)
DLASV2 computes the singular value decomposition of a 2-by-2 triangular matrix.
Definition dlasv2.f:136
subroutine drot(n, dx, incx, dy, incy, c, s)
DROT
Definition drot.f:92