LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
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 complex16GEcomputational
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, dlabad, 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  CALL dlabad( safmin, safmax )
301  ulp = dlamch( 'PRECISION' )
302  smlnum = safmin*( dble( n )/ulp )
303 
304  IF ( ihi .EQ. kwtop ) THEN
305 * 1 by 1 deflation window, just try a regular deflation
306  alpha( kwtop ) = a( kwtop, kwtop )
307  beta( kwtop ) = b( kwtop, kwtop )
308  ns = 1
309  nd = 0
310  IF ( abs( s ) .LE. max( smlnum, ulp*abs( a( kwtop,
311  $ kwtop ) ) ) ) THEN
312  ns = 0
313  nd = 1
314  IF ( kwtop .GT. ilo ) THEN
315  a( kwtop, kwtop-1 ) = czero
316  END IF
317  END IF
318  END IF
319 
320 
321 * Store window in case of convergence failure
322  CALL zlacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
323  CALL zlacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
324  $ 1 ), jw )
325 
326 * Transform window to real schur form
327  CALL zlaset( 'FULL', jw, jw, czero, cone, qc, ldqc )
328  CALL zlaset( 'FULL', jw, jw, czero, cone, zc, ldzc )
329  CALL zlaqz0( 'S', 'V', 'V', jw, 1, jw, a( kwtop, kwtop ), lda,
330  $ b( kwtop, kwtop ), ldb, alpha, beta, qc, ldqc, zc,
331  $ ldzc, work( 2*jw**2+1 ), lwork-2*jw**2, rwork,
332  $ rec+1, qz_small_info )
333 
334  IF( qz_small_info .NE. 0 ) THEN
335 * Convergence failure, restore the window and exit
336  nd = 0
337  ns = jw-qz_small_info
338  CALL zlacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), 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, c1,
396  $ s1 )
397  CALL zrot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
398  $ ldb, c1, s1 )
399  CALL zrot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
400  $ 1, c1, dconjg( s1 ) )
401  END DO
402 
403 * Chase bulges down
404  istartm = kwtop
405  istopm = ihi
406  k = kwbot-1
407  DO WHILE ( k .GE. kwtop )
408 
409 * Move bulge down and remove it
410  DO k2 = k, kwbot-1
411  CALL zlaqz1( .true., .true., k2, kwtop, kwtop+jw-1,
412  $ kwbot, a, lda, b, ldb, jw, kwtop, qc, ldqc,
413  $ jw, kwtop, zc, ldzc )
414  END DO
415 
416  k = k-1
417  END DO
418 
419  END IF
420 
421 * Apply Qc and Zc to rest of the matrix
422  IF ( ilschur ) THEN
423  istartm = 1
424  istopm = n
425  ELSE
426  istartm = ilo
427  istopm = ihi
428  END IF
429 
430  IF ( istopm-ihi > 0 ) THEN
431  CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
432  $ a( kwtop, ihi+1 ), lda, czero, work, jw )
433  CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
434  $ ihi+1 ), lda )
435  CALL zgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
436  $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
437  CALL zlacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
438  $ ihi+1 ), ldb )
439  END IF
440  IF ( ilq ) THEN
441  CALL zgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
442  $ ldqc, czero, work, n )
443  CALL zlacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
444  END IF
445 
446  IF ( kwtop-1-istartm+1 > 0 ) THEN
447  CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
448  $ kwtop ), lda, zc, ldzc, czero, work,
449  $ kwtop-istartm )
450  CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
451  $ a( istartm, kwtop ), lda )
452  CALL zgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
453  $ kwtop ), ldb, zc, ldzc, czero, work,
454  $ kwtop-istartm )
455  CALL zlacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
456  $ b( istartm, kwtop ), ldb )
457  END IF
458  IF ( ilz ) THEN
459  CALL zgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
460  $ ldzc, czero, work, n )
461  CALL zlacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
462  END IF
463 
464  END SUBROUTINE
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine zlartg(f, g, c, s, r)
ZLARTG generates a plane rotation with real cosine and complex sine.
Definition: zlartg.f90:118
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
Definition: zgemm.f:187
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 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 ztgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
ZTGEXC
Definition: ztgexc.f:200
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 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 zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:103
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