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