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