LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
zlaqz2.f
Go to the documentation of this file.
1*> \brief \b ZLAQZ2
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download ZLAQZ2 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ZLAQZ2.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ZLAQZ2.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ZLAQZ2.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B,
22* $ LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
23* $ WORK, LWORK, RWORK, REC, INFO )
24* IMPLICIT NONE
25*
26* Arguments
27* LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
28* INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
29* $ LDQC, LDZC, LWORK, REC
30*
31* COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
32* $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
33* INTEGER, INTENT( OUT ) :: NS, ND, INFO
34* COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
35* DOUBLE PRECISION :: RWORK( * )
36* ..
37*
38*
39*> \par Purpose:
40* =============
41*>
42*> \verbatim
43*>
44*> ZLAQZ2 performs AED
45*> \endverbatim
46*
47* Arguments:
48* ==========
49*
50*> \param[in] ILSCHUR
51*> \verbatim
52*> ILSCHUR is LOGICAL
53*> Determines whether or not to update the full Schur form
54*> \endverbatim
55*> \param[in] ILQ
56*> \verbatim
57*> ILQ is LOGICAL
58*> Determines whether or not to update the matrix Q
59*> \endverbatim
60*>
61*> \param[in] ILZ
62*> \verbatim
63*> ILZ is LOGICAL
64*> Determines whether or not to update the matrix Z
65*> \endverbatim
66*>
67*> \param[in] N
68*> \verbatim
69*> N is INTEGER
70*> The order of the matrices A, B, Q, and Z. N >= 0.
71*> \endverbatim
72*>
73*> \param[in] ILO
74*> \verbatim
75*> ILO is INTEGER
76*> \endverbatim
77*>
78*> \param[in] IHI
79*> \verbatim
80*> IHI is INTEGER
81*> ILO and IHI mark the rows and columns of (A,B) which
82*> are to be normalized
83*> \endverbatim
84*>
85*> \param[in] NW
86*> \verbatim
87*> NW is INTEGER
88*> The desired size of the deflation window.
89*> \endverbatim
90*>
91*> \param[in,out] A
92*> \verbatim
93*> A is COMPLEX*16 array, dimension (LDA, N)
94*> \endverbatim
95*>
96*> \param[in] LDA
97*> \verbatim
98*> LDA is INTEGER
99*> The leading dimension of the array A. LDA >= max( 1, N ).
100*> \endverbatim
101*>
102*> \param[in,out] B
103*> \verbatim
104*> B is COMPLEX*16 array, dimension (LDB, N)
105*> \endverbatim
106*>
107*> \param[in] LDB
108*> \verbatim
109*> LDB is INTEGER
110*> The leading dimension of the array B. LDB >= max( 1, N ).
111*> \endverbatim
112*>
113*> \param[in,out] Q
114*> \verbatim
115*> Q is COMPLEX*16 array, dimension (LDQ, N)
116*> \endverbatim
117*>
118*> \param[in] LDQ
119*> \verbatim
120*> LDQ is INTEGER
121*> \endverbatim
122*>
123*> \param[in,out] Z
124*> \verbatim
125*> Z is COMPLEX*16 array, dimension (LDZ, N)
126*> \endverbatim
127*>
128*> \param[in] LDZ
129*> \verbatim
130*> LDZ is INTEGER
131*> \endverbatim
132*>
133*> \param[out] NS
134*> \verbatim
135*> NS is INTEGER
136*> The number of unconverged eigenvalues available to
137*> use as shifts.
138*> \endverbatim
139*>
140*> \param[out] ND
141*> \verbatim
142*> ND is INTEGER
143*> The number of converged eigenvalues found.
144*> \endverbatim
145*>
146*> \param[out] ALPHA
147*> \verbatim
148*> ALPHA is COMPLEX*16 array, dimension (N)
149*> Each scalar alpha defining an eigenvalue
150*> of GNEP.
151*> \endverbatim
152*>
153*> \param[out] BETA
154*> \verbatim
155*> BETA is COMPLEX*16 array, dimension (N)
156*> The scalars beta that define the eigenvalues of GNEP.
157*> Together, the quantities alpha = ALPHA(j) and
158*> beta = BETA(j) represent the j-th eigenvalue of the matrix
159*> pair (A,B), in one of the forms lambda = alpha/beta or
160*> mu = beta/alpha. Since either lambda or mu may overflow,
161*> they should not, in general, be computed.
162*> \endverbatim
163*>
164*> \param[in,out] QC
165*> \verbatim
166*> QC is COMPLEX*16 array, dimension (LDQC, NW)
167*> \endverbatim
168*>
169*> \param[in] LDQC
170*> \verbatim
171*> LDQC is INTEGER
172*> \endverbatim
173*>
174*> \param[in,out] ZC
175*> \verbatim
176*> ZC is COMPLEX*16 array, dimension (LDZC, NW)
177*> \endverbatim
178*>
179*> \param[in] LDZC
180*> \verbatim
181*> LDZ is INTEGER
182*> \endverbatim
183*>
184*> \param[out] WORK
185*> \verbatim
186*> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
187*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
188*> \endverbatim
189*>
190*> \param[in] LWORK
191*> \verbatim
192*> LWORK is INTEGER
193*> The dimension of the array WORK. LWORK >= max(1,N).
194*>
195*> If LWORK = -1, then a workspace query is assumed; the routine
196*> only calculates the optimal size of the WORK array, returns
197*> this value as the first entry of the WORK array, and no error
198*> message related to LWORK is issued by XERBLA.
199*> \endverbatim
200*>
201*> \param[out] RWORK
202*> \verbatim
203*> RWORK is DOUBLE PRECISION array, dimension (N)
204*> \endverbatim
205*>
206*> \param[in] REC
207*> \verbatim
208*> REC is INTEGER
209*> REC indicates the current recursion level. Should be set
210*> to 0 on first call.
211*> \endverbatim
212*>
213*> \param[out] INFO
214*> \verbatim
215*> INFO is INTEGER
216*> = 0: successful exit
217*> < 0: if INFO = -i, the i-th argument had an illegal value
218*> \endverbatim
219*
220* Authors:
221* ========
222*
223*> \author Thijs Steel, KU Leuven
224*
225*> \date May 2020
226*
227*> \ingroup laqz2
228*>
229* =====================================================================
230 RECURSIVE SUBROUTINE zlaqz2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
231 $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
232 $ ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
233 $ WORK, LWORK, RWORK, REC, INFO )
234 IMPLICIT NONE
235
236* Arguments
237 LOGICAL, INTENT( IN ) :: ilschur, ilq, ilz
238 INTEGER, INTENT( IN ) :: n, ilo, ihi, nw, lda, ldb, ldq, ldz,
239 $ ldqc, ldzc, lwork, rec
240
241 COMPLEX*16, INTENT( INOUT ) :: a( lda, * ), b( ldb, * ), q( ldq,
242 $ * ), z( ldz, * ), alpha( * ), beta( * )
243 INTEGER, INTENT( OUT ) :: ns, nd, info
244 COMPLEX*16 :: qc( ldqc, * ), zc( ldzc, * ), work( * )
245 DOUBLE PRECISION :: rwork( * )
246
247* Parameters
248 COMPLEX*16 czero, cone
249 PARAMETER ( czero = ( 0.0d+0, 0.0d+0 ), cone = ( 1.0d+0,
250 $ 0.0d+0 ) )
251 DOUBLE PRECISION :: zero, one, half
252 parameter( zero = 0.0d0, one = 1.0d0, half = 0.5d0 )
253
254* Local Scalars
255 INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, ztgexc_info,
256 $ ifst, ilst, lworkreq, qz_small_info
257 DOUBLE PRECISION ::smlnum, ulp, safmin, safmax, c1, tempr
258 COMPLEX*16 :: s, s1, temp
259
260* External Functions
261 EXTERNAL :: xerbla, zlaqz0, zlaqz1, zlacpy, zlaset, zgemm,
262 $ ztgexc, zlartg, zrot
263 DOUBLE PRECISION, EXTERNAL :: dlamch
264
265 info = 0
266
267* Set up deflation window
268 jw = min( nw, ihi-ilo+1 )
269 kwtop = ihi-jw+1
270 IF ( kwtop .EQ. ilo ) THEN
271 s = czero
272 ELSE
273 s = a( kwtop, kwtop-1 )
274 END IF
275
276* Determine required workspace
277 ifst = 1
278 ilst = jw
279 CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
280 $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
281 $ ldzc, work, -1, rwork, rec+1, qz_small_info )
282 lworkreq = int( work( 1 ) )+2*jw**2
283 lworkreq = max( lworkreq, n*nw, 2*nw**2+n )
284 IF ( lwork .EQ.-1 ) THEN
285* workspace query, quick return
286 work( 1 ) = lworkreq
287 RETURN
288 ELSE IF ( lwork .LT. lworkreq ) THEN
289 info = -26
290 END IF
291
292 IF( info.NE.0 ) THEN
293 CALL xerbla( 'ZLAQZ2', -info )
294 RETURN
295 END IF
296
297* Get machine constants
298 safmin = dlamch( 'SAFE MINIMUM' )
299 safmax = one/safmin
300 ulp = dlamch( 'PRECISION' )
301 smlnum = safmin*( dble( n )/ulp )
302
303 IF ( ihi .EQ. kwtop ) THEN
304* 1 by 1 deflation window, just try a regular deflation
305 alpha( kwtop ) = a( kwtop, kwtop )
306 beta( kwtop ) = b( kwtop, kwtop )
307 ns = 1
308 nd = 0
309 IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
310 $ kwtop ) ) ) ) THEN
311 ns = 0
312 nd = 1
313 IF ( kwtop .GT. ilo ) THEN
314 a( kwtop, kwtop-1 ) = czero
315 END IF
316 END IF
317 END IF
318
319
320* Store window in case of convergence failure
321 CALL zlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
322 CALL zlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
323 $ 1 ), jw )
324
325* Transform window to real schur form
326 CALL zlaset( 'FULL', jw, jw, czero, cone, qc, ldqc )
327 CALL zlaset( 'FULL', jw, jw, czero, cone, zc, ldzc )
328 CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
329 $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
330 $ ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
331 $ rec+1, qz_small_info )
332
333 IF( qz_small_info .NE. 0 ) THEN
334* Convergence failure, restore the window and exit
335 nd = 0
336 ns = jw-qz_small_info
337 CALL zlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
338 CALL zlacpy( 'ALL', jw, jw, work( jw**2+1 ), jw, b( kwtop,
339 $ kwtop ), ldb )
340 RETURN
341 END IF
342
343* Deflation detection loop
344 IF ( kwtop .EQ. ilo .OR. s .EQ. czero ) THEN
345 kwbot = kwtop-1
346 ELSE
347 kwbot = ihi
348 k = 1
349 k2 = 1
350 DO WHILE ( k .LE. jw )
351* Try to deflate eigenvalue
352 tempr = abs( a( kwbot, kwbot ) )
353 IF( tempr .EQ. zero ) THEN
354 tempr = abs( s )
355 END IF
356 IF ( ( abs( s*qc( 1, kwbot-kwtop+1 ) ) ) .LE. max( ulp*
357 $ tempr, smlnum ) ) THEN
358* Deflatable
359 kwbot = kwbot-1
360 ELSE
361* Not deflatable, move out of the way
362 ifst = kwbot-kwtop+1
363 ilst = k2
364 CALL ztgexc( .true., .true., jw, a( kwtop, kwtop ),
365 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
366 $ zc, ldzc, ifst, ilst, ztgexc_info )
367 k2 = k2+1
368 END IF
369
370 k = k+1
371 END DO
372 END IF
373
374* Store eigenvalues
375 nd = ihi-kwbot
376 ns = jw-nd
377 k = kwtop
378 DO WHILE ( k .LE. ihi )
379 alpha( k ) = a( k, k )
380 beta( k ) = b( k, k )
381 k = k+1
382 END DO
383
384 IF ( kwtop .NE. ilo .AND. s .NE. czero ) THEN
385* Reflect spike back, this will create optimally packed bulges
386 a( kwtop:kwbot, kwtop-1 ) = a( kwtop, kwtop-1 ) *dconjg( qc( 1,
387 $ 1:jw-nd ) )
388 DO k = kwbot-1, kwtop, -1
389 CALL zlartg( a( k, kwtop-1 ), a( k+1, kwtop-1 ), c1, s1,
390 $ temp )
391 a( k, kwtop-1 ) = temp
392 a( k+1, kwtop-1 ) = czero
393 k2 = max( kwtop, k-1 )
394 CALL zrot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
395 $ s1 )
396 CALL zrot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
397 $ ldb, c1, s1 )
398 CALL zrot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
399 $ 1, c1, dconjg( s1 ) )
400 END DO
401
402* Chase bulges down
403 istartm = kwtop
404 istopm = ihi
405 k = kwbot-1
406 DO WHILE ( k .GE. kwtop )
407
408* Move bulge down and remove it
409 DO k2 = k, kwbot-1
410 CALL zlaqz1( .true., .true., k2, kwtop, kwtop+jw-1,
411 $ kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
412 $ jw, kwtop, zc, ldzc )
413 END DO
414
415 k = k-1
416 END DO
417
418 END IF
419
420* Apply Qc and Zc to rest of the matrix
421 IF ( ilschur ) THEN
422 istartm = 1
423 istopm = n
424 ELSE
425 istartm = ilo
426 istopm = ihi
427 END IF
428
429 IF ( istopm-ihi > 0 ) THEN
430 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
431 $ a( kwtop, ihi+1 ), lda, czero, work, jw )
432 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
433 $ ihi+1 ), lda )
434 CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
435 $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
436 CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
437 $ ihi+1 ), ldb )
438 END IF
439 IF ( ilq ) THEN
440 CALL zgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
441 $ ldqc, czero, work, n )
442 CALL zlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
443 END IF
444
445 IF ( kwtop-1-istartm+1 > 0 ) THEN
446 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
447 $ kwtop ), lda, zc, ldzc, czero, work,
448 $ kwtop-istartm )
449 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
450 $ a( istartm, kwtop ), lda )
451 CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
452 $ kwtop ), ldb, zc, ldzc, czero, work,
453 $ kwtop-istartm )
454 CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
455 $ b( istartm, kwtop ), ldb )
456 END IF
457 IF ( ilz ) THEN
458 CALL zgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
459 $ ldzc, czero, work, n )
460 CALL zlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
461 END IF
462
463 END SUBROUTINE
subroutine xerbla(srname, info)
Definition cblat2.f:3285
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
Definition zgemm.f:188
subroutine zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
Definition zlacpy.f:103
double precision function dlamch(cmach)
DLAMCH
Definition dlamch.f:69
recursive subroutine zlaqz0(wants, wantq, wantz, n, ilo, ihi, a, lda, b, ldb, alpha, beta, q, ldq, z, ldz, work, lwork, rwork, rec, info)
ZLAQZ0
Definition zlaqz0.f:284
subroutine zlaqz1(ilq, ilz, k, istartm, istopm, ihi, a, lda, b, ldb, nq, qstart, q, ldq, nz, zstart, z, ldz)
ZLAQZ1
Definition zlaqz1.f:173
recursive subroutine zlaqz2(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)
ZLAQZ2
Definition zlaqz2.f:234
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition zlartg.f90:116
subroutine zlaset(uplo, m, n, alpha, beta, a, lda)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition zlaset.f:106
subroutine zrot(n, cx, incx, cy, incy, c, s)
ZROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors.
Definition zrot.f:103
subroutine ztgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, info)
ZTGEXC
Definition ztgexc.f:200