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