LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
claqz3.f
Go to the documentation of this file.
1*> \brief \b CLAQZ3
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> Download CLAQZ3 + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/CLAQZ3.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/CLAQZ3.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/CLAQZ3.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE CLAQZ3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
20* $ NBLOCK_DESIRED, ALPHA, BETA, A, LDA, B, LDB, Q, LDQ, Z, LDZ,
21* $ QC, LDQC, ZC, LDZC, WORK, LWORK, INFO )
22* IMPLICIT NONE
23*
24* Function arguments
25* LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
26* INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
27* $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
28*
29* COMPLEX, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
30* $ Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
31* $ ALPHA( * ), BETA( * )
32*
33* INTEGER, INTENT( OUT ) :: INFO
34* ..
35*
36*
37*> \par Purpose:
38* =============
39*>
40*> \verbatim
41*>
42*> CLAQZ3 Executes a single multishift QZ sweep
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*> \endverbatim
80*>
81*> \param[in] NSHIFTS
82*> \verbatim
83*> NSHIFTS is INTEGER
84*> The desired number of shifts to use
85*> \endverbatim
86*>
87*> \param[in] NBLOCK_DESIRED
88*> \verbatim
89*> NBLOCK_DESIRED is INTEGER
90*> The desired size of the computational windows
91*> \endverbatim
92*>
93*> \param[in] ALPHA
94*> \verbatim
95*> ALPHA is COMPLEX array. SR contains
96*> the alpha parts of the shifts to use.
97*> \endverbatim
98*>
99*> \param[in] BETA
100*> \verbatim
101*> BETA is COMPLEX array. SS contains
102*> the scale of the shifts to use.
103*> \endverbatim
104*>
105*> \param[in,out] A
106*> \verbatim
107*> A is COMPLEX array, dimension (LDA, N)
108*> \endverbatim
109*>
110*> \param[in] LDA
111*> \verbatim
112*> LDA is INTEGER
113*> The leading dimension of the array A. LDA >= max( 1, N ).
114*> \endverbatim
115*>
116*> \param[in,out] B
117*> \verbatim
118*> B is COMPLEX array, dimension (LDB, N)
119*> \endverbatim
120*>
121*> \param[in] LDB
122*> \verbatim
123*> LDB is INTEGER
124*> The leading dimension of the array B. LDB >= max( 1, N ).
125*> \endverbatim
126*>
127*> \param[in,out] Q
128*> \verbatim
129*> Q is COMPLEX array, dimension (LDQ, N)
130*> \endverbatim
131*>
132*> \param[in] LDQ
133*> \verbatim
134*> LDQ is INTEGER
135*> \endverbatim
136*>
137*> \param[in,out] Z
138*> \verbatim
139*> Z is COMPLEX array, dimension (LDZ, N)
140*> \endverbatim
141*>
142*> \param[in] LDZ
143*> \verbatim
144*> LDZ is INTEGER
145*> \endverbatim
146*>
147*> \param[in,out] QC
148*> \verbatim
149*> QC is COMPLEX array, dimension (LDQC, NBLOCK_DESIRED)
150*> \endverbatim
151*>
152*> \param[in] LDQC
153*> \verbatim
154*> LDQC is INTEGER
155*> \endverbatim
156*>
157*> \param[in,out] ZC
158*> \verbatim
159*> ZC is COMPLEX array, dimension (LDZC, NBLOCK_DESIRED)
160*> \endverbatim
161*>
162*> \param[in] LDZC
163*> \verbatim
164*> LDZ is INTEGER
165*> \endverbatim
166*>
167*> \param[out] WORK
168*> \verbatim
169*> WORK is COMPLEX array, dimension (MAX(1,LWORK))
170*> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
171*> \endverbatim
172*>
173*> \param[in] LWORK
174*> \verbatim
175*> LWORK is INTEGER
176*> The dimension of the array WORK. LWORK >= max(1,N).
177*>
178*> If LWORK = -1, then a workspace query is assumed; the routine
179*> only calculates the optimal size of the WORK array, returns
180*> this value as the first entry of the WORK array, and no error
181*> message related to LWORK is issued by XERBLA.
182*> \endverbatim
183*>
184*> \param[out] INFO
185*> \verbatim
186*> INFO is INTEGER
187*> = 0: successful exit
188*> < 0: if INFO = -i, the i-th argument had an illegal value
189*> \endverbatim
190*
191* Authors:
192* ========
193*
194*> \author Thijs Steel, KU Leuven
195*
196*> \date May 2020
197*
198*> \ingroup laqz3
199*>
200* =====================================================================
201 SUBROUTINE claqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
202 $ NBLOCK_DESIRED, ALPHA, BETA, A, LDA, B, LDB,
203 $ Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK,
204 $ LWORK, INFO )
205 IMPLICIT NONE
206
207* Function arguments
208 LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
209 INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
210 $ NSHIFTS, NBLOCK_DESIRED, LDQC, LDZC
211
212 COMPLEX, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ, * ),
213 $ Z( LDZ, * ), QC( LDQC, * ), ZC( LDZC, * ), WORK( * ),
214 $ ALPHA( * ), BETA( * )
215
216 INTEGER, INTENT( OUT ) :: INFO
217
218* Parameters
219 COMPLEX CZERO, CONE
220 PARAMETER ( CZERO = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
221 REAL :: ZERO, ONE, HALF
222 parameter( zero = 0.0, one = 1.0, half = 0.5 )
223
224* Local scalars
225 INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
226 $ ISTARTB, ISTOPB, ISHIFT, NBLOCK, NPOS
227 REAL :: SAFMIN, SAFMAX, C, SCALE
228 COMPLEX :: TEMP, TEMP2, TEMP3, S
229
230* External Functions
231 EXTERNAL :: xerbla, claset, clartg, crot,
232 $ claqz1, cgemm, clacpy
233 REAL, EXTERNAL :: SLAMCH
234
235 info = 0
236 IF ( nblock_desired .LT. nshifts+1 ) THEN
237 info = -8
238 END IF
239 IF ( lwork .EQ.-1 ) THEN
240* workspace query, quick return
241 work( 1 ) = cmplx( n*nblock_desired )
242 RETURN
243 ELSE IF ( lwork .LT. n*nblock_desired ) THEN
244 info = -25
245 END IF
246
247 IF( info.NE.0 ) THEN
248 CALL xerbla( 'CLAQZ3', -info )
249 RETURN
250 END IF
251
252*
253* Executable statements
254*
255
256* Get machine constants
257 safmin = slamch( 'SAFE MINIMUM' )
258 safmax = one/safmin
259
260 IF ( ilo .GE. ihi ) THEN
261 RETURN
262 END IF
263
264 IF ( ilschur ) THEN
265 istartm = 1
266 istopm = n
267 ELSE
268 istartm = ilo
269 istopm = ihi
270 END IF
271
272 ns = nshifts
273 npos = max( nblock_desired-ns, 1 )
274
275
276* The following block introduces the shifts and chases
277* them down one by one just enough to make space for
278* the other shifts. The near-the-diagonal block is
279* of size (ns+1) x ns.
280
281 CALL claset( 'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
282 CALL claset( 'FULL', ns, ns, czero, cone, zc, ldzc )
283
284 DO i = 1, ns
285* Introduce the shift
286 scale = sqrt( abs( alpha( i ) ) ) * sqrt( abs( beta( i ) ) )
287 IF( scale .GE. safmin .AND. scale .LE. safmax ) THEN
288 alpha( i ) = alpha( i )/scale
289 beta( i ) = beta( i )/scale
290 END IF
291
292 temp2 = beta( i )*a( ilo, ilo )-alpha( i )*b( ilo, ilo )
293 temp3 = beta( i )*a( ilo+1, ilo )
294
295 IF ( abs( temp2 ) .GT. safmax .OR.
296 $ abs( temp3 ) .GT. safmax ) THEN
297 temp2 = cone
298 temp3 = czero
299 END IF
300
301 CALL clartg( temp2, temp3, c, s, temp )
302 CALL crot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda,
303 $ c, s )
304 CALL crot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb,
305 $ c, s )
306 CALL crot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c,
307 $ conjg( s ) )
308
309* Chase the shift down
310 DO j = 1, ns-i
311
312 CALL claqz1( .true., .true., j, 1, ns, ihi-ilo+1, a( ilo,
313 $ ilo ), lda, b( ilo, ilo ), ldb, ns+1, 1, qc,
314 $ ldqc, ns, 1, zc, ldzc )
315
316 END DO
317
318 END DO
319
320* Update the rest of the pencil
321
322* Update A(ilo:ilo+ns,ilo+ns:istopm) and B(ilo:ilo+ns,ilo+ns:istopm)
323* from the left with Qc(1:ns+1,1:ns+1)'
324 sheight = ns+1
325 swidth = istopm-( ilo+ns )+1
326 IF ( swidth > 0 ) THEN
327 CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
328 $ ldqc,
329 $ a( ilo, ilo+ns ), lda, czero, work, sheight )
330 CALL clacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
331 $ ilo+ns ), lda )
332 CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
333 $ ldqc,
334 $ b( ilo, ilo+ns ), ldb, czero, work, sheight )
335 CALL clacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
336 $ ilo+ns ), ldb )
337 END IF
338 IF ( ilq ) THEN
339 CALL cgemm( 'N', 'N', n, sheight, sheight, cone, q( 1, ilo ),
340 $ ldq, qc, ldqc, czero, work, n )
341 CALL clacpy( 'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
342 END IF
343
344* Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
345* from the right with Zc(1:ns,1:ns)
346 sheight = ilo-1-istartm+1
347 swidth = ns
348 IF ( sheight > 0 ) THEN
349 CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
350 $ a( istartm, ilo ), lda, zc, ldzc, czero, work,
351 $ sheight )
352 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
353 $ a( istartm,
354 $ ilo ), lda )
355 CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
356 $ b( istartm, ilo ), ldb, zc, ldzc, czero, work,
357 $ sheight )
358 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
359 $ b( istartm,
360 $ ilo ), ldb )
361 END IF
362 IF ( ilz ) THEN
363 CALL cgemm( 'N', 'N', n, swidth, swidth, cone, z( 1, ilo ),
364 $ ldz, zc, ldzc, czero, work, n )
365 CALL clacpy( 'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
366 END IF
367
368* The following block chases the shifts down to the bottom
369* right block. If possible, a shift is moved down npos
370* positions at a time
371
372 k = ilo
373 DO WHILE ( k < ihi-ns )
374 np = min( ihi-ns-k, npos )
375* Size of the near-the-diagonal block
376 nblock = ns+np
377* istartb points to the first row we will be updating
378 istartb = k+1
379* istopb points to the last column we will be updating
380 istopb = k+nblock-1
381
382 CALL claset( 'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
383 CALL claset( 'FULL', ns+np, ns+np, czero, cone, zc, ldzc )
384
385* Near the diagonal shift chase
386 DO i = ns-1, 0, -1
387 DO j = 0, np-1
388* Move down the block with index k+i+j, updating
389* the (ns+np x ns+np) block:
390* (k:k+ns+np,k:k+ns+np-1)
391 CALL claqz1( .true., .true., k+i+j, istartb, istopb,
392 $ ihi,
393 $ a, lda, b, ldb, nblock, k+1, qc, ldqc,
394 $ nblock, k, zc, ldzc )
395 END DO
396 END DO
397
398* Update rest of the pencil
399
400* Update A(k+1:k+ns+np, k+ns+np:istopm) and
401* B(k+1:k+ns+np, k+ns+np:istopm)
402* from the left with Qc(1:ns+np,1:ns+np)'
403 sheight = ns+np
404 swidth = istopm-( k+ns+np )+1
405 IF ( swidth > 0 ) THEN
406 CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
407 $ ldqc, a( k+1, k+ns+np ), lda, czero, work,
408 $ sheight )
409 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
410 $ a( k+1,
411 $ k+ns+np ), lda )
412 CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
413 $ ldqc, b( k+1, k+ns+np ), ldb, czero, work,
414 $ sheight )
415 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
416 $ b( k+1,
417 $ k+ns+np ), ldb )
418 END IF
419 IF ( ilq ) THEN
420 CALL cgemm( 'N', 'N', n, nblock, nblock, cone, q( 1,
421 $ k+1 ),
422 $ ldq, qc, ldqc, czero, work, n )
423 CALL clacpy( 'ALL', n, nblock, work, n, q( 1, k+1 ),
424 $ ldq )
425 END IF
426
427* Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
428* from the right with Zc(1:ns+np,1:ns+np)
429 sheight = k-istartm+1
430 swidth = nblock
431 IF ( sheight > 0 ) THEN
432 CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
433 $ a( istartm, k ), lda, zc, ldzc, czero, work,
434 $ sheight )
435 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
436 $ a( istartm, k ), lda )
437 CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
438 $ b( istartm, k ), ldb, zc, ldzc, czero, work,
439 $ sheight )
440 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
441 $ b( istartm, k ), ldb )
442 END IF
443 IF ( ilz ) THEN
444 CALL cgemm( 'N', 'N', n, nblock, nblock, cone, z( 1, k ),
445 $ ldz, zc, ldzc, czero, work, n )
446 CALL clacpy( 'ALL', n, nblock, work, n, z( 1, k ), ldz )
447 END IF
448
449 k = k+np
450
451 END DO
452
453* The following block removes the shifts from the bottom right corner
454* one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
455
456 CALL claset( 'FULL', ns, ns, czero, cone, qc, ldqc )
457 CALL claset( 'FULL', ns+1, ns+1, czero, cone, zc, ldzc )
458
459* istartb points to the first row we will be updating
460 istartb = ihi-ns+1
461* istopb points to the last column we will be updating
462 istopb = ihi
463
464 DO i = 1, ns
465* Chase the shift down to the bottom right corner
466 DO ishift = ihi-i, ihi-1
467 CALL claqz1( .true., .true., ishift, istartb, istopb,
468 $ ihi,
469 $ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
470 $ ihi-ns, zc, ldzc )
471 END DO
472
473 END DO
474
475* Update rest of the pencil
476
477* Update A(ihi-ns+1:ihi, ihi+1:istopm)
478* from the left with Qc(1:ns,1:ns)'
479 sheight = ns
480 swidth = istopm-( ihi+1 )+1
481 IF ( swidth > 0 ) THEN
482 CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
483 $ ldqc,
484 $ a( ihi-ns+1, ihi+1 ), lda, czero, work, sheight )
485 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
486 $ a( ihi-ns+1, ihi+1 ), lda )
487 CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
488 $ ldqc,
489 $ b( ihi-ns+1, ihi+1 ), ldb, czero, work, sheight )
490 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
491 $ b( ihi-ns+1, ihi+1 ), ldb )
492 END IF
493 IF ( ilq ) THEN
494 CALL cgemm( 'N', 'N', n, ns, ns, cone, q( 1, ihi-ns+1 ),
495 $ ldq,
496 $ qc, ldqc, czero, work, n )
497 CALL clacpy( 'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
498 END IF
499
500* Update A(istartm:ihi-ns,ihi-ns:ihi)
501* from the right with Zc(1:ns+1,1:ns+1)
502 sheight = ihi-ns-istartm+1
503 swidth = ns+1
504 IF ( sheight > 0 ) THEN
505 CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
506 $ a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
507 $ sheight )
508 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
509 $ a( istartm,
510 $ ihi-ns ), lda )
511 CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
512 $ b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
513 $ sheight )
514 CALL clacpy( 'ALL', sheight, swidth, work, sheight,
515 $ b( istartm,
516 $ ihi-ns ), ldb )
517 END IF
518 IF ( ilz ) THEN
519 CALL cgemm( 'N', 'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ),
520 $ ldz,
521 $ zc, ldzc, czero, work, n )
522 CALL clacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
523 END IF
524
525 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
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
claqz3
subroutine claqz3(ilschur, ilq, ilz, n, ilo, ihi, nshifts, nblock_desired, alpha, beta, a, lda, b, ldb, q, ldq, z, ldz, qc, ldqc, zc, ldzc, work, lwork, info)
CLAQZ3
Definition claqz3.f:205
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