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