LAPACK 3.12.1
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*> Download CLAQZ2 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/CLAQZ2.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/CLAQZ2.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/CLAQZ2.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLAQZ2( 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, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30* $ Z( LDZ, * ), ALPHA( * ), BETA( * )
31* INTEGER, INTENT( OUT ) :: NS, ND, INFO
32* COMPLEX :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
33* REAL :: RWORK( * )
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> CLAQZ2 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*>
54*> \param[in] ILQ
55*> \verbatim
56*> ILQ is LOGICAL
57*> Determines whether or not to update the matrix Q
58*> \endverbatim
59*>
60*> \param[in] ILZ
61*> \verbatim
62*> ILZ is LOGICAL
63*> Determines whether or not to update the matrix Z
64*> \endverbatim
65*>
66*> \param[in] N
67*> \verbatim
68*> N is INTEGER
69*> The order of the matrices A, B, Q, and Z. N >= 0.
70*> \endverbatim
71*>
72*> \param[in] ILO
73*> \verbatim
74*> ILO is INTEGER
75*> \endverbatim
76*>
77*> \param[in] IHI
78*> \verbatim
79*> IHI is INTEGER
80*> ILO and IHI mark the rows and columns of (A,B) which
81*> are to be normalized
82*> \endverbatim
83*>
84*> \param[in] NW
85*> \verbatim
86*> NW is INTEGER
87*> The desired size of the deflation window.
88*> \endverbatim
89*>
90*> \param[in,out] A
91*> \verbatim
92*> A is COMPLEX array, dimension (LDA, N)
93*> \endverbatim
94*>
95*> \param[in] LDA
96*> \verbatim
97*> LDA is INTEGER
98*> The leading dimension of the array A. LDA >= max( 1, N ).
99*> \endverbatim
100*>
101*> \param[in,out] B
102*> \verbatim
103*> B is COMPLEX array, dimension (LDB, N)
104*> \endverbatim
105*>
106*> \param[in] LDB
107*> \verbatim
108*> LDB is INTEGER
109*> The leading dimension of the array B. LDB >= max( 1, N ).
110*> \endverbatim
111*>
112*> \param[in,out] Q
113*> \verbatim
114*> Q is COMPLEX array, dimension (LDQ, N)
115*> \endverbatim
116*>
117*> \param[in] LDQ
118*> \verbatim
119*> LDQ is INTEGER
120*> \endverbatim
121*>
122*> \param[in,out] Z
123*> \verbatim
124*> Z is COMPLEX array, dimension (LDZ, N)
125*> \endverbatim
126*>
127*> \param[in] LDZ
128*> \verbatim
129*> LDZ is INTEGER
130*> \endverbatim
131*>
132*> \param[out] NS
133*> \verbatim
134*> NS is INTEGER
135*> The number of unconverged eigenvalues available to
136*> use as shifts.
137*> \endverbatim
138*>
139*> \param[out] ND
140*> \verbatim
141*> ND is INTEGER
142*> The number of converged eigenvalues found.
143*> \endverbatim
144*>
145*> \param[out] ALPHA
146*> \verbatim
147*> ALPHA is COMPLEX array, dimension (N)
148*> Each scalar alpha defining an eigenvalue
149*> of GNEP.
150*> \endverbatim
151*>
152*> \param[out] BETA
153*> \verbatim
154*> BETA is COMPLEX array, dimension (N)
155*> The scalars beta that define the eigenvalues of GNEP.
156*> Together, the quantities alpha = ALPHA(j) and
157*> beta = BETA(j) represent the j-th eigenvalue of the matrix
158*> pair (A,B), in one of the forms lambda = alpha/beta or
159*> mu = beta/alpha. Since either lambda or mu may overflow,
160*> they should not, in general, be computed.
161*> \endverbatim
162*>
163*> \param[in,out] QC
164*> \verbatim
165*> QC is COMPLEX array, dimension (LDQC, NW)
166*> \endverbatim
167*>
168*> \param[in] LDQC
169*> \verbatim
170*> LDQC is INTEGER
171*> \endverbatim
172*>
173*> \param[in,out] ZC
174*> \verbatim
175*> ZC is COMPLEX array, dimension (LDZC, NW)
176*> \endverbatim
177*>
178*> \param[in] LDZC
179*> \verbatim
180*> LDZ is INTEGER
181*> \endverbatim
182*>
183*> \param[out] WORK
184*> \verbatim
185*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
186*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
187*> \endverbatim
188*>
189*> \param[in] LWORK
190*> \verbatim
191*> LWORK is INTEGER
192*> The dimension of the array WORK. LWORK >= max(1,N).
193*>
194*> If LWORK = -1, then a workspace query is assumed; the routine
195*> only calculates the optimal size of the WORK array, returns
196*> this value as the first entry of the WORK array, and no error
197*> message related to LWORK is issued by XERBLA.
198*> \endverbatim
199*>
200*> \param[out] RWORK
201*> \verbatim
202*> RWORK is REAL array, dimension (N)
203*> \endverbatim
204*>
205*> \param[in] REC
206*> \verbatim
207*> REC is INTEGER
208*> REC indicates the current recursion level. Should be set
209*> to 0 on first call.
210*> \endverbatim
211*>
212*> \param[out] INFO
213*> \verbatim
214*> INFO is INTEGER
215*> = 0: successful exit
216*> < 0: if INFO = -i, the i-th argument had an illegal value
217*> \endverbatim
218*
219* Authors:
220* ========
221*
222*> \author Thijs Steel, KU Leuven, KU Leuven
223*
224*> \date May 2020
225*
226*> \ingroup laqz2
227*>
228* =====================================================================
229 RECURSIVE SUBROUTINE claqz2( ILSCHUR, ILQ, ILZ, N, ILO, IHI,
230 $ 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, 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 ) = cmplx( 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 ulp = slamch( 'PRECISION' )
300 smlnum = safmin*( real( 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 clacpy( 'ALL', jw, jw, a( kwtop, kwtop ), lda, work, jw )
321 CALL clacpy( 'ALL', jw, jw, b( kwtop, kwtop ), ldb,
322 $ 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 ),
338 $ lda )
339 CALL clacpy( '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 ctgexc( .true., .true., jw, a( kwtop, kwtop ),
366 $ lda, b( kwtop, kwtop ), ldb, qc, ldqc,
367 $ zc, ldzc, ifst, ilst, ctgexc_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 ) *conjg( qc( 1,
388 $ 1:jw-nd ) )
389 DO k = kwbot-1, kwtop, -1
390 CALL clartg( 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 crot( ihi-k2+1, a( k, k2 ), lda, a( k+1, k2 ), lda,
396 $ c1,
397 $ s1 )
398 CALL crot( ihi-( k-1 )+1, b( k, k-1 ), ldb, b( k+1,
399 $ k-1 ),
400 $ ldb, c1, s1 )
401 CALL crot( jw, qc( 1, k-kwtop+1 ), 1, qc( 1,
402 $ k+1-kwtop+1 ),
403 $ 1, c1, conjg( 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 claqz1( .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 cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
435 $ a( kwtop, ihi+1 ), lda, czero, work, jw )
436 CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, a( kwtop,
437 $ ihi+1 ), lda )
438 CALL cgemm( 'C', 'N', jw, istopm-ihi, jw, cone, qc, ldqc,
439 $ b( kwtop, ihi+1 ), ldb, czero, work, jw )
440 CALL clacpy( 'ALL', jw, istopm-ihi, work, jw, b( kwtop,
441 $ ihi+1 ), ldb )
442 END IF
443 IF ( ilq ) THEN
444 CALL cgemm( 'N', 'N', n, jw, jw, cone, q( 1, kwtop ), ldq,
445 $ qc,
446 $ ldqc, czero, work, n )
447 CALL clacpy( 'ALL', n, jw, work, n, q( 1, kwtop ), ldq )
448 END IF
449
450 IF ( kwtop-1-istartm+1 > 0 ) THEN
451 CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone,
452 $ a( istartm,
453 $ kwtop ), lda, zc, ldzc, czero, work,
454 $ kwtop-istartm )
455 CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
456 $ a( istartm, kwtop ), lda )
457 CALL cgemm( 'N', 'N', kwtop-istartm, jw, jw, cone,
458 $ b( istartm,
459 $ kwtop ), ldb, zc, ldzc, czero, work,
460 $ kwtop-istartm )
461 CALL clacpy( 'ALL', kwtop-istartm, jw, work, kwtop-istartm,
462 $ b( istartm, kwtop ), ldb )
463 END IF
464 IF ( ilz ) THEN
465 CALL cgemm( 'N', 'N', n, jw, jw, cone, z( 1, kwtop ), ldz,
466 $ zc,
467 $ ldzc, czero, work, n )
468 CALL clacpy( '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
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:101
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:283
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:172
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:234
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:104
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:101
ctgexc
subroutine ctgexc(wantq, wantz, n, a, lda, b, ldb, q, ldq, z, ldz, ifst, ilst, info)
CTGEXC
Definition ctgexc.f:198