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