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