LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
dlaqz0.f
Go to the documentation of this file.
1 *> \brief \b DLAQZ0
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DLAQZ0 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaqz0.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaqz0.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaqz0.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * RECURSIVE SUBROUTINE DLAQZ0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
22 * $ LDA, B, LDB, ALPHAR, ALPHAI, BETA,
23 * $ Q, LDQ, Z, LDZ, WORK, LWORK, REC,
24 * $ INFO )
25 * IMPLICIT NONE
26 *
27 * Arguments
28 * CHARACTER, INTENT( IN ) :: WANTS, WANTQ, WANTZ
29 * INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
30 * $ REC
31 *
32 * INTEGER, INTENT( OUT ) :: INFO
33 *
34 * DOUBLE PRECISION, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ),
35 * $ Q( LDQ, * ), Z( LDZ, * ), ALPHAR( * ),
36 * $ ALPHAI( * ), BETA( * ), WORK( * )
37 * ..
38 *
39 *
40 *> \par Purpose:
41 * =============
42 *>
43 *> \verbatim
44 *>
45 *> DLAQZ0 computes the eigenvalues of a real matrix pair (H,T),
46 *> where H is an upper Hessenberg matrix and T is upper triangular,
47 *> using the double-shift QZ method.
48 *> Matrix pairs of this type are produced by the reduction to
49 *> generalized upper Hessenberg form of a real matrix pair (A,B):
50 *>
51 *> A = Q1*H*Z1**T, B = Q1*T*Z1**T,
52 *>
53 *> as computed by DGGHRD.
54 *>
55 *> If JOB='S', then the Hessenberg-triangular pair (H,T) is
56 *> also reduced to generalized Schur form,
57 *>
58 *> H = Q*S*Z**T, T = Q*P*Z**T,
59 *>
60 *> where Q and Z are orthogonal matrices, P is an upper triangular
61 *> matrix, and S is a quasi-triangular matrix with 1-by-1 and 2-by-2
62 *> diagonal blocks.
63 *>
64 *> The 1-by-1 blocks correspond to real eigenvalues of the matrix pair
65 *> (H,T) and the 2-by-2 blocks correspond to complex conjugate pairs of
66 *> eigenvalues.
67 *>
68 *> Additionally, the 2-by-2 upper triangular diagonal blocks of P
69 *> corresponding to 2-by-2 blocks of S are reduced to positive diagonal
70 *> form, i.e., if S(j+1,j) is non-zero, then P(j+1,j) = P(j,j+1) = 0,
71 *> P(j,j) > 0, and P(j+1,j+1) > 0.
72 *>
73 *> Optionally, the orthogonal matrix Q from the generalized Schur
74 *> factorization may be postmultiplied into an input matrix Q1, and the
75 *> orthogonal matrix Z may be postmultiplied into an input matrix Z1.
76 *> If Q1 and Z1 are the orthogonal matrices from DGGHRD that reduced
77 *> the matrix pair (A,B) to generalized upper Hessenberg form, then the
78 *> output matrices Q1*Q and Z1*Z are the orthogonal factors from the
79 *> generalized Schur factorization of (A,B):
80 *>
81 *> A = (Q1*Q)*S*(Z1*Z)**T, B = (Q1*Q)*P*(Z1*Z)**T.
82 *>
83 *> To avoid overflow, eigenvalues of the matrix pair (H,T) (equivalently,
84 *> of (A,B)) are computed as a pair of values (alpha,beta), where alpha is
85 *> complex and beta real.
86 *> If beta is nonzero, lambda = alpha / beta is an eigenvalue of the
87 *> generalized nonsymmetric eigenvalue problem (GNEP)
88 *> A*x = lambda*B*x
89 *> and if alpha is nonzero, mu = beta / alpha is an eigenvalue of the
90 *> alternate form of the GNEP
91 *> mu*A*y = B*y.
92 *> Real eigenvalues can be read directly from the generalized Schur
93 *> form:
94 *> alpha = S(i,i), beta = P(i,i).
95 *>
96 *> Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
97 *> Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
98 *> pp. 241--256.
99 *>
100 *> Ref: B. Kagstrom, D. Kressner, "Multishift Variants of the QZ
101 *> Algorithm with Aggressive Early Deflation", SIAM J. Numer.
102 *> Anal., 29(2006), pp. 199--227.
103 *>
104 *> Ref: T. Steel, D. Camps, K. Meerbergen, R. Vandebril "A multishift,
105 *> multipole rational QZ method with agressive early deflation"
106 *> \endverbatim
107 *
108 * Arguments:
109 * ==========
110 *
111 *> \param[in] WANTS
112 *> \verbatim
113 *> WANTS is CHARACTER*1
114 *> = 'E': Compute eigenvalues only;
115 *> = 'S': Compute eigenvalues and the Schur form.
116 *> \endverbatim
117 *>
118 *> \param[in] WANTQ
119 *> \verbatim
120 *> WANTQ is CHARACTER*1
121 *> = 'N': Left Schur vectors (Q) are not computed;
122 *> = 'I': Q is initialized to the unit matrix and the matrix Q
123 *> of left Schur vectors of (A,B) is returned;
124 *> = 'V': Q must contain an orthogonal matrix Q1 on entry and
125 *> the product Q1*Q is returned.
126 *> \endverbatim
127 *>
128 *> \param[in] WANTZ
129 *> \verbatim
130 *> WANTZ is CHARACTER*1
131 *> = 'N': Right Schur vectors (Z) are not computed;
132 *> = 'I': Z is initialized to the unit matrix and the matrix Z
133 *> of right Schur vectors of (A,B) is returned;
134 *> = 'V': Z must contain an orthogonal matrix Z1 on entry and
135 *> the product Z1*Z is returned.
136 *> \endverbatim
137 *>
138 *> \param[in] N
139 *> \verbatim
140 *> N is INTEGER
141 *> The order of the matrices A, B, Q, and Z. N >= 0.
142 *> \endverbatim
143 *>
144 *> \param[in] ILO
145 *> \verbatim
146 *> ILO is INTEGER
147 *> \endverbatim
148 *>
149 *> \param[in] IHI
150 *> \verbatim
151 *> IHI is INTEGER
152 *> ILO and IHI mark the rows and columns of A which are in
153 *> Hessenberg form. It is assumed that A is already upper
154 *> triangular in rows and columns 1:ILO-1 and IHI+1:N.
155 *> If N > 0, 1 <= ILO <= IHI <= N; if N = 0, ILO=1 and IHI=0.
156 *> \endverbatim
157 *>
158 *> \param[in,out] A
159 *> \verbatim
160 *> A is DOUBLE PRECISION array, dimension (LDA, N)
161 *> On entry, the N-by-N upper Hessenberg matrix A.
162 *> On exit, if JOB = 'S', A contains the upper quasi-triangular
163 *> matrix S from the generalized Schur factorization.
164 *> If JOB = 'E', the diagonal blocks of A match those of S, but
165 *> the rest of A is unspecified.
166 *> \endverbatim
167 *>
168 *> \param[in] LDA
169 *> \verbatim
170 *> LDA is INTEGER
171 *> The leading dimension of the array A. LDA >= max( 1, N ).
172 *> \endverbatim
173 *>
174 *> \param[in,out] B
175 *> \verbatim
176 *> B is DOUBLE PRECISION array, dimension (LDB, N)
177 *> On entry, the N-by-N upper triangular matrix B.
178 *> On exit, if JOB = 'S', B contains the upper triangular
179 *> matrix P from the generalized Schur factorization;
180 *> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks of S
181 *> are reduced to positive diagonal form, i.e., if A(j+1,j) is
182 *> non-zero, then B(j+1,j) = B(j,j+1) = 0, B(j,j) > 0, and
183 *> B(j+1,j+1) > 0.
184 *> If JOB = 'E', the diagonal blocks of B match those of P, but
185 *> the rest of B is unspecified.
186 *> \endverbatim
187 *>
188 *> \param[in] LDB
189 *> \verbatim
190 *> LDB is INTEGER
191 *> The leading dimension of the array B. LDB >= max( 1, N ).
192 *> \endverbatim
193 *>
194 *> \param[out] ALPHAR
195 *> \verbatim
196 *> ALPHAR is DOUBLE PRECISION array, dimension (N)
197 *> The real parts of each scalar alpha defining an eigenvalue
198 *> of GNEP.
199 *> \endverbatim
200 *>
201 *> \param[out] ALPHAI
202 *> \verbatim
203 *> ALPHAI is DOUBLE PRECISION array, dimension (N)
204 *> The imaginary parts of each scalar alpha defining an
205 *> eigenvalue of GNEP.
206 *> If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
207 *> positive, then the j-th and (j+1)-st eigenvalues are a
208 *> complex conjugate pair, with ALPHAI(j+1) = -ALPHAI(j).
209 *> \endverbatim
210 *>
211 *> \param[out] BETA
212 *> \verbatim
213 *> BETA is DOUBLE PRECISION array, dimension (N)
214 *> The scalars beta that define the eigenvalues of GNEP.
215 *> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
216 *> beta = BETA(j) represent the j-th eigenvalue of the matrix
217 *> pair (A,B), in one of the forms lambda = alpha/beta or
218 *> mu = beta/alpha. Since either lambda or mu may overflow,
219 *> they should not, in general, be computed.
220 *> \endverbatim
221 *>
222 *> \param[in,out] Q
223 *> \verbatim
224 *> Q is DOUBLE PRECISION array, dimension (LDQ, N)
225 *> On entry, if COMPQ = 'V', the orthogonal matrix Q1 used in
226 *> the reduction of (A,B) to generalized Hessenberg form.
227 *> On exit, if COMPQ = 'I', the orthogonal matrix of left Schur
228 *> vectors of (A,B), and if COMPQ = 'V', the orthogonal matrix
229 *> of left Schur vectors of (A,B).
230 *> Not referenced if COMPQ = 'N'.
231 *> \endverbatim
232 *>
233 *> \param[in] LDQ
234 *> \verbatim
235 *> LDQ is INTEGER
236 *> The leading dimension of the array Q. LDQ >= 1.
237 *> If COMPQ='V' or 'I', then LDQ >= N.
238 *> \endverbatim
239 *>
240 *> \param[in,out] Z
241 *> \verbatim
242 *> Z is DOUBLE PRECISION array, dimension (LDZ, N)
243 *> On entry, if COMPZ = 'V', the orthogonal matrix Z1 used in
244 *> the reduction of (A,B) to generalized Hessenberg form.
245 *> On exit, if COMPZ = 'I', the orthogonal matrix of
246 *> right Schur vectors of (H,T), and if COMPZ = 'V', the
247 *> orthogonal matrix of right Schur vectors of (A,B).
248 *> Not referenced if COMPZ = 'N'.
249 *> \endverbatim
250 *>
251 *> \param[in] LDZ
252 *> \verbatim
253 *> LDZ is INTEGER
254 *> The leading dimension of the array Z. LDZ >= 1.
255 *> If COMPZ='V' or 'I', then LDZ >= N.
256 *> \endverbatim
257 *>
258 *> \param[out] WORK
259 *> \verbatim
260 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
261 *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
262 *> \endverbatim
263 *>
264 *> \param[in] LWORK
265 *> \verbatim
266 *> LWORK is INTEGER
267 *> The dimension of the array WORK. LWORK >= max(1,N).
268 *>
269 *> If LWORK = -1, then a workspace query is assumed; the routine
270 *> only calculates the optimal size of the WORK array, returns
271 *> this value as the first entry of the WORK array, and no error
272 *> message related to LWORK is issued by XERBLA.
273 *> \endverbatim
274 *>
275 *> \param[in] REC
276 *> \verbatim
277 *> REC is INTEGER
278 *> REC indicates the current recursion level. Should be set
279 *> to 0 on first call.
280 *> \endverbatim
281 *>
282 *> \param[out] INFO
283 *> \verbatim
284 *> INFO is INTEGER
285 *> = 0: successful exit
286 *> < 0: if INFO = -i, the i-th argument had an illegal value
287 *> = 1,...,N: the QZ iteration did not converge. (A,B) is not
288 *> in Schur form, but ALPHAR(i), ALPHAI(i), and
289 *> BETA(i), i=INFO+1,...,N should be correct.
290 *> \endverbatim
291 *
292 * Authors:
293 * ========
294 *
295 *> \author Thijs Steel, KU Leuven
296 *
297 *> \date May 2020
298 *
299 *> \ingroup doubleGEcomputational
300 *>
301 * =====================================================================
302  RECURSIVE SUBROUTINE dlaqz0( WANTS, WANTQ, WANTZ, N, ILO, IHI, A,
303  $ LDA, B, LDB, ALPHAR, ALPHAI, BETA,
304  $ Q, LDQ, Z, LDZ, WORK, LWORK, REC,
305  $ INFO )
306  IMPLICIT NONE
307 
308 * Arguments
309  CHARACTER, INTENT( IN ) :: wants, wantq, wantz
310  INTEGER, INTENT( IN ) :: n, ilo, ihi, lda, ldb, ldq, ldz, lwork,
311  $ rec
312 
313  INTEGER, INTENT( OUT ) :: info
314 
315  DOUBLE PRECISION, INTENT( INOUT ) :: a( lda, * ), b( ldb, * ),
316  $ q( ldq, * ), z( ldz, * ), alphar( * ),
317  $ alphai( * ), beta( * ), work( * )
318 
319 * Parameters
320  DOUBLE PRECISION :: zero, one, half
321  PARAMETER( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
322 
323 * Local scalars
324  DOUBLE PRECISION :: smlnum, ulp, eshift, safmin, safmax, c1, s1,
325  $ temp, swap
326  INTEGER :: istart, istop, iiter, maxit, istart2, k, ld, nshifts,
327  $ nblock, nw, nmin, nibble, n_undeflated, n_deflated,
328  $ ns, sweep_info, shiftpos, lworkreq, k2, istartm,
329  $ istopm, iwants, iwantq, iwantz, norm_info, aed_info,
330  $ nwr, nbr, nsr, itemp1, itemp2, rcost, i
331  LOGICAL :: ilschur, ilq, ilz
332  CHARACTER :: jbcmpz*3
333 
334 * External Functions
335  EXTERNAL :: xerbla, dhgeqz, dlaset, dlaqz3, dlaqz4, dlabad,
336  $ dlartg, drot
337  DOUBLE PRECISION, EXTERNAL :: dlamch
338  LOGICAL, EXTERNAL :: lsame
339  INTEGER, EXTERNAL :: ilaenv
340 
341 *
342 * Decode wantS,wantQ,wantZ
343 *
344  IF( lsame( wants, 'E' ) ) THEN
345  ilschur = .false.
346  iwants = 1
347  ELSE IF( lsame( wants, 'S' ) ) THEN
348  ilschur = .true.
349  iwants = 2
350  ELSE
351  iwants = 0
352  END IF
353 
354  IF( lsame( wantq, 'N' ) ) THEN
355  ilq = .false.
356  iwantq = 1
357  ELSE IF( lsame( wantq, 'V' ) ) THEN
358  ilq = .true.
359  iwantq = 2
360  ELSE IF( lsame( wantq, 'I' ) ) THEN
361  ilq = .true.
362  iwantq = 3
363  ELSE
364  iwantq = 0
365  END IF
366 
367  IF( lsame( wantz, 'N' ) ) THEN
368  ilz = .false.
369  iwantz = 1
370  ELSE IF( lsame( wantz, 'V' ) ) THEN
371  ilz = .true.
372  iwantz = 2
373  ELSE IF( lsame( wantz, 'I' ) ) THEN
374  ilz = .true.
375  iwantz = 3
376  ELSE
377  iwantz = 0
378  END IF
379 *
380 * Check Argument Values
381 *
382  info = 0
383  IF( iwants.EQ.0 ) THEN
384  info = -1
385  ELSE IF( iwantq.EQ.0 ) THEN
386  info = -2
387  ELSE IF( iwantz.EQ.0 ) THEN
388  info = -3
389  ELSE IF( n.LT.0 ) THEN
390  info = -4
391  ELSE IF( ilo.LT.1 ) THEN
392  info = -5
393  ELSE IF( ihi.GT.n .OR. ihi.LT.ilo-1 ) THEN
394  info = -6
395  ELSE IF( lda.LT.n ) THEN
396  info = -8
397  ELSE IF( ldb.LT.n ) THEN
398  info = -10
399  ELSE IF( ldq.LT.1 .OR. ( ilq .AND. ldq.LT.n ) ) THEN
400  info = -15
401  ELSE IF( ldz.LT.1 .OR. ( ilz .AND. ldz.LT.n ) ) THEN
402  info = -17
403  END IF
404  IF( info.NE.0 ) THEN
405  CALL xerbla( 'DLAQZ0', -info )
406  RETURN
407  END IF
408 
409 *
410 * Quick return if possible
411 *
412  IF( n.LE.0 ) THEN
413  work( 1 ) = dble( 1 )
414  RETURN
415  END IF
416 
417 *
418 * Get the parameters
419 *
420  jbcmpz( 1:1 ) = wants
421  jbcmpz( 2:2 ) = wantq
422  jbcmpz( 3:3 ) = wantz
423 
424  nmin = ilaenv( 12, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
425 
426  nwr = ilaenv( 13, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
427  nwr = max( 2, nwr )
428  nwr = min( ihi-ilo+1, ( n-1 ) / 3, nwr )
429 
430  nibble = ilaenv( 14, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
431 
432  nsr = ilaenv( 15, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
433  nsr = min( nsr, ( n+6 ) / 9, ihi-ilo )
434  nsr = max( 2, nsr-mod( nsr, 2 ) )
435 
436  rcost = ilaenv( 17, 'DLAQZ0', jbcmpz, n, ilo, ihi, lwork )
437  itemp1 = int( nsr/sqrt( 1+2*nsr/( dble( rcost )/100*n ) ) )
438  itemp1 = ( ( itemp1-1 )/4 )*4+4
439  nbr = nsr+itemp1
440 
441  IF( n .LT. nmin .OR. rec .GE. 2 ) THEN
442  CALL dhgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
443  $ alphar, alphai, beta, q, ldq, z, ldz, work,
444  $ lwork, info )
445  RETURN
446  END IF
447 
448 *
449 * Find out required workspace
450 *
451 
452 * Workspace query to dlaqz3
453  nw = max( nwr, nmin )
454  CALL dlaqz3( ilschur, ilq, ilz, n, ilo, ihi, nw, a, lda, b, ldb,
455  $ q, ldq, z, ldz, n_undeflated, n_deflated, alphar,
456  $ alphai, beta, work, nw, work, nw, work, -1, rec,
457  $ aed_info )
458  itemp1 = int( work( 1 ) )
459 * Workspace query to dlaqz4
460  CALL dlaqz4( ilschur, ilq, ilz, n, ilo, ihi, nsr, nbr, alphar,
461  $ alphai, beta, a, lda, b, ldb, q, ldq, z, ldz, work,
462  $ nbr, work, nbr, work, -1, sweep_info )
463  itemp2 = int( work( 1 ) )
464 
465  lworkreq = max( itemp1+2*nw**2, itemp2+2*nbr**2 )
466  IF ( lwork .EQ.-1 ) THEN
467  work( 1 ) = dble( lworkreq )
468  RETURN
469  ELSE IF ( lwork .LT. lworkreq ) THEN
470  info = -19
471  END IF
472  IF( info.NE.0 ) THEN
473  CALL xerbla( 'DLAQZ0', info )
474  RETURN
475  END IF
476 *
477 * Initialize Q and Z
478 *
479  IF( iwantq.EQ.3 ) CALL dlaset( 'FULL', n, n, zero, one, q, ldq )
480  IF( iwantz.EQ.3 ) CALL dlaset( 'FULL', n, n, zero, one, z, ldz )
481 
482 * Get machine constants
483  safmin = dlamch( 'SAFE MINIMUM' )
484  safmax = one/safmin
485  CALL dlabad( safmin, safmax )
486  ulp = dlamch( 'PRECISION' )
487  smlnum = safmin*( dble( n )/ulp )
488 
489  istart = ilo
490  istop = ihi
491  maxit = 3*( ihi-ilo+1 )
492  ld = 0
493 
494  DO iiter = 1, maxit
495  IF( iiter .GE. maxit ) THEN
496  info = istop+1
497  GOTO 80
498  END IF
499  IF ( istart+1 .GE. istop ) THEN
500  istop = istart
501  EXIT
502  END IF
503 
504 * Check deflations at the end
505  IF ( abs( a( istop-1, istop-2 ) ) .LE. max( smlnum,
506  $ ulp*( abs( a( istop-1, istop-1 ) )+abs( a( istop-2,
507  $ istop-2 ) ) ) ) ) THEN
508  a( istop-1, istop-2 ) = zero
509  istop = istop-2
510  ld = 0
511  eshift = zero
512  ELSE IF ( abs( a( istop, istop-1 ) ) .LE. max( smlnum,
513  $ ulp*( abs( a( istop, istop ) )+abs( a( istop-1,
514  $ istop-1 ) ) ) ) ) THEN
515  a( istop, istop-1 ) = zero
516  istop = istop-1
517  ld = 0
518  eshift = zero
519  END IF
520 * Check deflations at the start
521  IF ( abs( a( istart+2, istart+1 ) ) .LE. max( smlnum,
522  $ ulp*( abs( a( istart+1, istart+1 ) )+abs( a( istart+2,
523  $ istart+2 ) ) ) ) ) THEN
524  a( istart+2, istart+1 ) = zero
525  istart = istart+2
526  ld = 0
527  eshift = zero
528  ELSE IF ( abs( a( istart+1, istart ) ) .LE. max( smlnum,
529  $ ulp*( abs( a( istart, istart ) )+abs( a( istart+1,
530  $ istart+1 ) ) ) ) ) THEN
531  a( istart+1, istart ) = zero
532  istart = istart+1
533  ld = 0
534  eshift = zero
535  END IF
536 
537  IF ( istart+1 .GE. istop ) THEN
538  EXIT
539  END IF
540 
541 * Check interior deflations
542  istart2 = istart
543  DO k = istop, istart+1, -1
544  IF ( abs( a( k, k-1 ) ) .LE. max( smlnum, ulp*( abs( a( k,
545  $ k ) )+abs( a( k-1, k-1 ) ) ) ) ) THEN
546  a( k, k-1 ) = zero
547  istart2 = k
548  EXIT
549  END IF
550  END DO
551 
552 * Get range to apply rotations to
553  IF ( ilschur ) THEN
554  istartm = 1
555  istopm = n
556  ELSE
557  istartm = istart2
558  istopm = istop
559  END IF
560 
561 * Check infinite eigenvalues, this is done without blocking so might
562 * slow down the method when many infinite eigenvalues are present
563  k = istop
564  DO WHILE ( k.GE.istart2 )
565  temp = zero
566  IF( k .LT. istop ) THEN
567  temp = temp+abs( b( k, k+1 ) )
568  END IF
569  IF( k .GT. istart2 ) THEN
570  temp = temp+abs( b( k-1, k ) )
571  END IF
572 
573  IF( abs( b( k, k ) ) .LT. max( smlnum, ulp*temp ) ) THEN
574 * A diagonal element of B is negligable, move it
575 * to the top and deflate it
576 
577  DO k2 = k, istart2+1, -1
578  CALL dlartg( b( k2-1, k2 ), b( k2-1, k2-1 ), c1, s1,
579  $ temp )
580  b( k2-1, k2 ) = temp
581  b( k2-1, k2-1 ) = zero
582 
583  CALL drot( k2-2-istartm+1, b( istartm, k2 ), 1,
584  $ b( istartm, k2-1 ), 1, c1, s1 )
585  CALL drot( min( k2+1, istop )-istartm+1, a( istartm,
586  $ k2 ), 1, a( istartm, k2-1 ), 1, c1, s1 )
587  IF ( ilz ) THEN
588  CALL drot( n, z( 1, k2 ), 1, z( 1, k2-1 ), 1, c1,
589  $ s1 )
590  END IF
591 
592  IF( k2.LT.istop ) THEN
593  CALL dlartg( a( k2, k2-1 ), a( k2+1, k2-1 ), c1,
594  $ s1, temp )
595  a( k2, k2-1 ) = temp
596  a( k2+1, k2-1 ) = zero
597 
598  CALL drot( istopm-k2+1, a( k2, k2 ), lda, a( k2+1,
599  $ k2 ), lda, c1, s1 )
600  CALL drot( istopm-k2+1, b( k2, k2 ), ldb, b( k2+1,
601  $ k2 ), ldb, c1, s1 )
602  IF( ilq ) THEN
603  CALL drot( n, q( 1, k2 ), 1, q( 1, k2+1 ), 1,
604  $ c1, s1 )
605  END IF
606  END IF
607 
608  END DO
609 
610  IF( istart2.LT.istop )THEN
611  CALL dlartg( a( istart2, istart2 ), a( istart2+1,
612  $ istart2 ), c1, s1, temp )
613  a( istart2, istart2 ) = temp
614  a( istart2+1, istart2 ) = zero
615 
616  CALL drot( istopm-( istart2+1 )+1, a( istart2,
617  $ istart2+1 ), lda, a( istart2+1,
618  $ istart2+1 ), lda, c1, s1 )
619  CALL drot( istopm-( istart2+1 )+1, b( istart2,
620  $ istart2+1 ), ldb, b( istart2+1,
621  $ istart2+1 ), ldb, c1, s1 )
622  IF( ilq ) THEN
623  CALL drot( n, q( 1, istart2 ), 1, q( 1,
624  $ istart2+1 ), 1, c1, s1 )
625  END IF
626  END IF
627 
628  istart2 = istart2+1
629 
630  END IF
631  k = k-1
632  END DO
633 
634 * istart2 now points to the top of the bottom right
635 * unreduced Hessenberg block
636  IF ( istart2 .GE. istop ) THEN
637  istop = istart2-1
638  ld = 0
639  eshift = zero
640  cycle
641  END IF
642 
643  nw = nwr
644  nshifts = nsr
645  nblock = nbr
646 
647  IF ( istop-istart2+1 .LT. nmin ) THEN
648 * Setting nw to the size of the subblock will make AED deflate
649 * all the eigenvalues. This is slightly more efficient than just
650 * using DHGEQZ because the off diagonal part gets updated via BLAS.
651  IF ( istop-istart+1 .LT. nmin ) THEN
652  nw = istop-istart+1
653  istart2 = istart
654  ELSE
655  nw = istop-istart2+1
656  END IF
657  END IF
658 
659 *
660 * Time for AED
661 *
662  CALL dlaqz3( ilschur, ilq, ilz, n, istart2, istop, nw, a, lda,
663  $ b, ldb, q, ldq, z, ldz, n_undeflated, n_deflated,
664  $ alphar, alphai, beta, work, nw, work( nw**2+1 ),
665  $ nw, work( 2*nw**2+1 ), lwork-2*nw**2, rec,
666  $ aed_info )
667 
668  IF ( n_deflated > 0 ) THEN
669  istop = istop-n_deflated
670  ld = 0
671  eshift = zero
672  END IF
673 
674  IF ( 100*n_deflated > nibble*( n_deflated+n_undeflated ) .OR.
675  $ istop-istart2+1 .LT. nmin ) THEN
676 * AED has uncovered many eigenvalues. Skip a QZ sweep and run
677 * AED again.
678  cycle
679  END IF
680 
681  ld = ld+1
682 
683  ns = min( nshifts, istop-istart2 )
684  ns = min( ns, n_undeflated )
685  shiftpos = istop-n_deflated-n_undeflated+1
686 *
687 * Shuffle shifts to put double shifts in front
688 * This ensures that we don't split up a double shift
689 *
690  DO i = shiftpos, shiftpos+n_undeflated-1, 2
691  IF( alphai( i ).NE.-alphai( i+1 ) ) THEN
692 *
693  swap = alphar( i )
694  alphar( i ) = alphar( i+1 )
695  alphar( i+1 ) = alphar( i+2 )
696  alphar( i+2 ) = swap
697 
698  swap = alphai( i )
699  alphai( i ) = alphai( i+1 )
700  alphai( i+1 ) = alphai( i+2 )
701  alphai( i+2 ) = swap
702 
703  swap = beta( i )
704  beta( i ) = beta( i+1 )
705  beta( i+1 ) = beta( i+2 )
706  beta( i+2 ) = swap
707  END IF
708  END DO
709 
710  IF ( mod( ld, 6 ) .EQ. 0 ) THEN
711 *
712 * Exceptional shift. Chosen for no particularly good reason.
713 *
714  IF( ( dble( maxit )*safmin )*abs( a( istop,
715  $ istop-1 ) ).LT.abs( a( istop-1, istop-1 ) ) ) THEN
716  eshift = a( istop, istop-1 )/b( istop-1, istop-1 )
717  ELSE
718  eshift = eshift+one/( safmin*dble( maxit ) )
719  END IF
720  alphar( shiftpos ) = one
721  alphar( shiftpos+1 ) = zero
722  alphai( shiftpos ) = zero
723  alphai( shiftpos+1 ) = zero
724  beta( shiftpos ) = eshift
725  beta( shiftpos+1 ) = eshift
726  ns = 2
727  END IF
728 
729 *
730 * Time for a QZ sweep
731 *
732  CALL dlaqz4( ilschur, ilq, ilz, n, istart2, istop, ns, nblock,
733  $ alphar( shiftpos ), alphai( shiftpos ),
734  $ beta( shiftpos ), a, lda, b, ldb, q, ldq, z, ldz,
735  $ work, nblock, work( nblock**2+1 ), nblock,
736  $ work( 2*nblock**2+1 ), lwork-2*nblock**2,
737  $ sweep_info )
738 
739  END DO
740 
741 *
742 * Call DHGEQZ to normalize the eigenvalue blocks and set the eigenvalues
743 * If all the eigenvalues have been found, DHGEQZ will not do any iterations
744 * and only normalize the blocks. In case of a rare convergence failure,
745 * the single shift might perform better.
746 *
747  80 CALL dhgeqz( wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb,
748  $ alphar, alphai, beta, q, ldq, z, ldz, work, lwork,
749  $ norm_info )
750 
751  info = norm_info
752 
753  END SUBROUTINE
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine dlartg(f, g, c, s, r)
DLARTG generates a plane rotation with real cosine and real sine.
Definition: dlartg.f90:113
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
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
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
recursive subroutine dlaqz0(WANTS, WANTQ, WANTZ, N, ILO, IHI, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, REC, INFO)
DLAQZ0
Definition: dlaqz0.f:306
subroutine dlaqz4(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS, NBLOCK_DESIRED, SR, SI, SS, A, LDA, B, LDB, Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK, LWORK, INFO)
DLAQZ4
Definition: dlaqz4.f:213
recursive subroutine dlaqz3(ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHAR, ALPHAI, BETA, QC, LDQC, ZC, LDZC, WORK, LWORK, REC, INFO)
DLAQZ3
Definition: dlaqz3.f:239