LAPACK  3.10.1
LAPACK: Linear Algebra PACKage
claqz2.f
Go to the documentation of this file.
1 *> \brief \b CLAQZ2
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLAQZ2 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/CLAQZ2.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/CLAQZ2.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/CLAQZ2.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLAQZ2( 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, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
32 * $ Z( LDZ, * ), ALPHA( * ), BETA( * )
33 * INTEGER, INTENT( OUT ) :: NS, ND, INFO
34 * COMPLEX :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
35 * REAL :: RWORK( * )
36 * ..
37 *
38 *
39 *> \par Purpose:
40 * =============
41 *>
42 *> \verbatim
43 *>
44 *> CLAQZ2 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 *>
56 *> \param[in] ILQ
57 *> \verbatim
58 *> ILQ is LOGICAL
59 *> Determines whether or not to update the matrix Q
60 *> \endverbatim
61 *>
62 *> \param[in] ILZ
63 *> \verbatim
64 *> ILZ is LOGICAL
65 *> Determines whether or not to update the matrix Z
66 *> \endverbatim
67 *>
68 *> \param[in] N
69 *> \verbatim
70 *> N is INTEGER
71 *> The order of the matrices A, B, Q, and Z. N >= 0.
72 *> \endverbatim
73 *>
74 *> \param[in] ILO
75 *> \verbatim
76 *> ILO is INTEGER
77 *> \endverbatim
78 *>
79 *> \param[in] IHI
80 *> \verbatim
81 *> IHI is INTEGER
82 *> ILO and IHI mark the rows and columns of (A,B) which
83 *> are to be normalized
84 *> \endverbatim
85 *>
86 *> \param[in] NW
87 *> \verbatim
88 *> NW is INTEGER
89 *> The desired size of the deflation window.
90 *> \endverbatim
91 *>
92 *> \param[in,out] A
93 *> \verbatim
94 *> A is COMPLEX array, dimension (LDA, N)
95 *> \endverbatim
96 *>
97 *> \param[in] LDA
98 *> \verbatim
99 *> LDA is INTEGER
100 *> The leading dimension of the array A. LDA >= max( 1, N ).
101 *> \endverbatim
102 *>
103 *> \param[in,out] B
104 *> \verbatim
105 *> B is COMPLEX array, dimension (LDB, N)
106 *> \endverbatim
107 *>
108 *> \param[in] LDB
109 *> \verbatim
110 *> LDB is INTEGER
111 *> The leading dimension of the array B. LDB >= max( 1, N ).
112 *> \endverbatim
113 *>
114 *> \param[in,out] Q
115 *> \verbatim
116 *> Q is COMPLEX array, dimension (LDQ, N)
117 *> \endverbatim
118 *>
119 *> \param[in] LDQ
120 *> \verbatim
121 *> LDQ is INTEGER
122 *> \endverbatim
123 *>
124 *> \param[in,out] Z
125 *> \verbatim
126 *> Z is COMPLEX array, dimension (LDZ, N)
127 *> \endverbatim
128 *>
129 *> \param[in] LDZ
130 *> \verbatim
131 *> LDZ is INTEGER
132 *> \endverbatim
133 *>
134 *> \param[out] NS
135 *> \verbatim
136 *> NS is INTEGER
137 *> The number of unconverged eigenvalues available to
138 *> use as shifts.
139 *> \endverbatim
140 *>
141 *> \param[out] ND
142 *> \verbatim
143 *> ND is INTEGER
144 *> The number of converged eigenvalues found.
145 *> \endverbatim
146 *>
147 *> \param[out] ALPHA
148 *> \verbatim
149 *> ALPHA is COMPLEX array, dimension (N)
150 *> Each scalar alpha defining an eigenvalue
151 *> of GNEP.
152 *> \endverbatim
153 *>
154 *> \param[out] BETA
155 *> \verbatim
156 *> BETA is COMPLEX array, dimension (N)
157 *> The scalars beta that define the eigenvalues of GNEP.
158 *> Together, the quantities alpha = ALPHA(j) and
159 *> beta = BETA(j) represent the j-th eigenvalue of the matrix
160 *> pair (A,B), in one of the forms lambda = alpha/beta or
161 *> mu = beta/alpha. Since either lambda or mu may overflow,
162 *> they should not, in general, be computed.
163 *> \endverbatim
164 *>
165 *> \param[in,out] QC
166 *> \verbatim
167 *> QC is COMPLEX array, dimension (LDQC, NW)
168 *> \endverbatim
169 *>
170 *> \param[in] LDQC
171 *> \verbatim
172 *> LDQC is INTEGER
173 *> \endverbatim
174 *>
175 *> \param[in,out] ZC
176 *> \verbatim
177 *> ZC is COMPLEX array, dimension (LDZC, NW)
178 *> \endverbatim
179 *>
180 *> \param[in] LDZC
181 *> \verbatim
182 *> LDZ is INTEGER
183 *> \endverbatim
184 *>
185 *> \param[out] WORK
186 *> \verbatim
187 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
188 *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
189 *> \endverbatim
190 *>
191 *> \param[in] LWORK
192 *> \verbatim
193 *> LWORK is INTEGER
194 *> The dimension of the array WORK. LWORK >= max(1,N).
195 *>
196 *> If LWORK = -1, then a workspace query is assumed; the routine
197 *> only calculates the optimal size of the WORK array, returns
198 *> this value as the first entry of the WORK array, and no error
199 *> message related to LWORK is issued by XERBLA.
200 *> \endverbatim
201 *>
202 *> \param[out] RWORK
203 *> \verbatim
204 *> RWORK is REAL array, dimension (N)
205 *> \endverbatim
206 *>
207 *> \param[in] REC
208 *> \verbatim
209 *> REC is INTEGER
210 *> REC indicates the current recursion level. Should be set
211 *> to 0 on first call.
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, KU Leuven
224 *
225 *> \date May 2020
226 *
227 *> \ingroup complexGEcomputational
228 *>
229 * =====================================================================
230  RECURSIVE SUBROUTINE claqz2( 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, INTENT( INOUT ) :: a( lda, * ), b( ldb, * ), q( ldq, * ),
242  $ z( ldz, * ), alpha( * ), beta( * )
243  INTEGER, INTENT( OUT ) :: ns, nd, info
244  COMPLEX :: qc( ldqc, * ), zc( ldzc, * ), work( * )
245  REAL :: rwork( * )
246 
247 * Parameters
248  COMPLEX czero, cone
249  PARAMETER ( czero = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
250  REAL :: zero, one, half
251  parameter( zero = 0.0, one = 1.0, half = 0.5 )
252 
253 * Local Scalars
254  INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, ctgexc_info,
255  $ ifst, ilst, lworkreq, qz_small_info
256  REAL :: smlnum, ulp, safmin, safmax, c1, tempr
257  COMPLEX :: s, s1, temp
258 
259 * External Functions
260  EXTERNAL :: xerbla, claqz0, claqz1, slabad, clacpy, claset, cgemm,
261  $ ctgexc, clartg, crot
262  REAL, EXTERNAL :: slamch
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 claqz0( '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( 'CLAQZ2', -info )
293  RETURN
294  END IF
295 
296 * Get machine constants
297  safmin = slamch( 'SAFE MINIMUM' )
298  safmax = one/safmin
299  CALL slabad( safmin, safmax )
300  ulp = slamch( 'PRECISION' )
301  smlnum = safmin*( real( 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 clacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
322  CALL clacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb, work( jw**2+
323  $ 1 ), jw )
324 
325 * Transform window to real schur form
326  CALL claset( 'FULL', jw, jw, czero, cone, qc, ldqc )
327  CALL claset( 'FULL', jw, jw, czero, cone, zc, ldzc )
328  CALL claqz0( '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 clacpy( 'ALL', jw, jw, work, jw, a( kwtop, kwtop ), lda )
338  CALL clacpy( '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 ctgexc( .true., .true., jw, a( kwtop, kwtop ),
365  $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
366  $ zc, ldzc, ifst, ilst, ctgexc_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 ) *conjg( qc( 1,
387  $ 1:jw-nd ) )
388  DO k = kwbot-1, kwtop, -1
389  CALL clartg( 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 crot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda, c1,
395  $ s1 )
396  CALL crot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1, k-1 ),
397  $ ldb, c1, s1 )
398  CALL crot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1, k+1-kwtop+1 ),
399  $ 1, c1, conjg( 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 claqz1( .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 cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
431  $ a( kwtop, ihi+1 ), lda, czero, work, jw )
432  CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
433  $ ihi+1 ), lda )
434  CALL cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
435  $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
436  CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
437  $ ihi+1 ), ldb )
438  END IF
439  IF ( ilq ) THEN
440  CALL cgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq, qc,
441  $ ldqc, czero, work, n )
442  CALL clacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
443  END IF
444 
445  IF ( kwtop-1-istartm+1 > 0 ) THEN
446  CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, a( istartm,
447  $ kwtop ), lda, zc, ldzc, czero, work,
448  $ kwtop-istartm )
449  CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
450  $ a( istartm, kwtop ), lda )
451  CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone, b( istartm,
452  $ kwtop ), ldb, zc, ldzc, czero, work,
453  $ kwtop-istartm )
454  CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
455  $ b( istartm, kwtop ), ldb )
456  END IF
457  IF ( ilz ) THEN
458  CALL cgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz, zc,
459  $ ldzc, czero, work, n )
460  CALL clacpy( 'ALL', n, jw, work, n, z( 1, kwtop ), ldz )
461  END IF
462 
463  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
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
subroutine claqz1(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
CLAQZ1
Definition: claqz1.f:173
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
subroutine ctgexc(WANTQ, WANTZ, N, A, LDA, B, LDB, Q, LDQ, Z, LDZ, IFST, ILST, INFO)
CTGEXC
Definition: ctgexc.f:200
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
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68