LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dlaqz4.f
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
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
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
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
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
dlaqz1
subroutine dlaqz1(A, LDA, B, LDB, SR1, SR2, SI, BETA1, BETA2, V)
DLAQZ1
Definition: dlaqz1.f:127
drot
subroutine drot(N, DX, INCX, DY, INCY, C, S)
DROT
Definition: drot.f:92
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
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