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