LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \endverbatim
213*>
214*> \param[out] INFO
215*> \verbatim
216*> INFO is INTEGER
217*> = 0: successful exit
218*> < 0: if INFO = -i, the i-th argument had an illegal value
219*> \endverbatim
220*
221* Authors:
222* ========
223*
224*> \author Thijs Steel, KU Leuven, KU Leuven
225*
226*> \date May 2020
227*
228*> \ingroup laqz2
229*>
230* =====================================================================
231 RECURSIVE SUBROUTINE claqz2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
232 $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
233 $ ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
234 $ WORK, LWORK, RWORK, REC, INFO )
235 IMPLICIT NONE
236
237* Arguments
238 LOGICAL, INTENT( IN ) :: ilschur, ilq, ilz
239 INTEGER, INTENT( IN ) :: n, ilo, ihi, nw, lda, ldb, ldq, ldz,
240 $ ldqc, ldzc, lwork, rec
241
242 COMPLEX, INTENT( INOUT ) :: a( lda, * ), b( ldb, * ), q( ldq, * ),
243 $ z( ldz, * ), alpha( * ), beta( * )
244 INTEGER, INTENT( OUT ) :: ns, nd, info
245 COMPLEX :: qc( ldqc, * ), zc( ldzc, * ), work( * )
246 REAL :: rwork( * )
247
248* Parameters
249 COMPLEX czero, cone
250 PARAMETER ( czero = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
251 REAL :: zero, one, half
252 parameter( zero = 0.0, one = 1.0, half = 0.5 )
253
254* Local Scalars
255 INTEGER :: jw, kwtop, kwbot, istopm, istartm, k, k2, ctgexc_info,
256 $ ifst, ilst, lworkreq, qz_small_info
257 REAL :: smlnum, ulp, safmin, safmax, c1, tempr
258 COMPLEX :: s, s1, temp
259
260* External Functions
261 EXTERNAL :: xerbla, claqz0, claqz1, clacpy, claset, cgemm,
262 $ ctgexc, clartg, crot
263 REAL, EXTERNAL :: slamch
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 claqz0( '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( 'CLAQZ2', -info )
294 RETURN
295 END IF
296
297* Get machine constants
298 safmin = slamch( 'SAFE MINIMUM' )
299 safmax = one/safmin
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
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
cgemm
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
Definition cgemm.f:188
clacpy
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
slamch
real function slamch(cmach)
SLAMCH
Definition slamch.f:68
claqz0
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
claqz1
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
claqz2
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:235
clartg
subroutine clartg(f, g, c, s, r)
CLARTG generates a plane rotation with real cosine and complex sine.
Definition clartg.f90:116
claset
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
crot
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
ctgexc
subroutine ctgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, info)
CTGEXC
Definition ctgexc.f:200