LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
claqz3.f
1 *> \brief \b CLAQZ3
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download CLAQZ3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/CLAQZ3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/CLAQZ3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/CLAQZ3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLAQZ3( 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, 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 *> CLAQZ3 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] ALPHA
96 *> \verbatim
97 *> ALPHA is COMPLEX array. SR contains
98 *> the alpha parts of the shifts to use.
99 *> \endverbatim
100 *>
101 *> \param[in] BETA
102 *> \verbatim
103 *> BETA is COMPLEX array. SS contains
104 *> the scale of the shifts to use.
105 *> \endverbatim
106 *>
107 *> \param[in,out] A
108 *> \verbatim
109 *> A is COMPLEX array, dimension (LDA, N)
110 *> \endverbatim
111 *>
112 *> \param[in] LDA
113 *> \verbatim
114 *> LDA is INTEGER
115 *> The leading dimension of the array A. LDA >= max( 1, N ).
116 *> \endverbatim
117 *>
118 *> \param[in,out] B
119 *> \verbatim
120 *> B is COMPLEX array, dimension (LDB, N)
121 *> \endverbatim
122 *>
123 *> \param[in] LDB
124 *> \verbatim
125 *> LDB is INTEGER
126 *> The leading dimension of the array B. LDB >= max( 1, N ).
127 *> \endverbatim
128 *>
129 *> \param[in,out] Q
130 *> \verbatim
131 *> Q is COMPLEX array, dimension (LDQ, N)
132 *> \endverbatim
133 *>
134 *> \param[in] LDQ
135 *> \verbatim
136 *> LDQ is INTEGER
137 *> \endverbatim
138 *>
139 *> \param[in,out] Z
140 *> \verbatim
141 *> Z is COMPLEX array, dimension (LDZ, N)
142 *> \endverbatim
143 *>
144 *> \param[in] LDZ
145 *> \verbatim
146 *> LDZ is INTEGER
147 *> \endverbatim
148 *>
149 *> \param[in,out] QC
150 *> \verbatim
151 *> QC is COMPLEX array, dimension (LDQC, NBLOCK_DESIRED)
152 *> \endverbatim
153 *>
154 *> \param[in] LDQC
155 *> \verbatim
156 *> LDQC is INTEGER
157 *> \endverbatim
158 *>
159 *> \param[in,out] ZC
160 *> \verbatim
161 *> ZC is COMPLEX array, dimension (LDZC, NBLOCK_DESIRED)
162 *> \endverbatim
163 *>
164 *> \param[in] LDZC
165 *> \verbatim
166 *> LDZ is INTEGER
167 *> \endverbatim
168 *>
169 *> \param[out] WORK
170 *> \verbatim
171 *> WORK is COMPLEX array, dimension (MAX(1,LWORK))
172 *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
173 *> \endverbatim
174 *>
175 *> \param[in] LWORK
176 *> \verbatim
177 *> LWORK is INTEGER
178 *> The dimension of the array WORK. LWORK >= max(1,N).
179 *>
180 *> If LWORK = -1, then a workspace query is assumed; the routine
181 *> only calculates the optimal size of the WORK array, returns
182 *> this value as the first entry of the WORK array, and no error
183 *> message related to LWORK is issued by XERBLA.
184 *> \endverbatim
185 *>
186 *> \param[out] INFO
187 *> \verbatim
188 *> INFO is INTEGER
189 *> = 0: successful exit
190 *> < 0: if INFO = -i, the i-th argument had an illegal value
191 *> \endverbatim
192 *
193 * Authors:
194 * ========
195 *
196 *> \author Thijs Steel, KU Leuven
197 *
198 *> \date May 2020
199 *
200 *> \ingroup laqz3
201 *>
202 * =====================================================================
203  SUBROUTINE claqz3( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NSHIFTS,
204  $ NBLOCK_DESIRED, ALPHA, BETA, A, LDA, B, LDB,
205  $ Q, LDQ, Z, LDZ, QC, LDQC, ZC, LDZC, WORK,
206  $ LWORK, INFO )
207  IMPLICIT NONE
208 
209 * Function arguments
210  LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
211  INTEGER, INTENT( IN ) :: N, ILO, IHI, LDA, LDB, LDQ, LDZ, LWORK,
212  $ nshifts, nblock_desired, ldqc, ldzc
213 
214  COMPLEX, INTENT( INOUT ) :: A( lda, * ), B( ldb, * ), Q( ldq, * ),
215  $ z( ldz, * ), qc( ldqc, * ), zc( ldzc, * ), work( * ),
216  $ alpha( * ), beta( * )
217 
218  INTEGER, INTENT( OUT ) :: INFO
219 
220 * Parameters
221  COMPLEX CZERO, CONE
222  parameter( czero = ( 0.0, 0.0 ), cone = ( 1.0, 0.0 ) )
223  REAL :: ZERO, ONE, HALF
224  parameter( zero = 0.0, one = 1.0, half = 0.5 )
225 
226 * Local scalars
227  INTEGER :: I, J, NS, ISTARTM, ISTOPM, SHEIGHT, SWIDTH, K, NP,
228  $ istartb, istopb, ishift, nblock, npos
229  REAL :: SAFMIN, SAFMAX, C, SCALE
230  COMPLEX :: TEMP, TEMP2, TEMP3, S
231 
232 * External Functions
233  EXTERNAL :: xerbla, slabad, claset, clartg, crot, claqz1, cgemm,
234  $ clacpy
235  REAL, EXTERNAL :: SLAMCH
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( 'CLAQZ3', -info )
251  RETURN
252  END IF
253 
254 *
255 * Executable statements
256 *
257 
258 * Get machine constants
259  safmin = slamch( 'SAFE MINIMUM' )
260  safmax = one/safmin
261  CALL slabad( safmin, safmax )
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  ns = nshifts
276  npos = max( nblock_desired-ns, 1 )
277 
278 
279 * The following block introduces the shifts and chases
280 * them down one by one just enough to make space for
281 * the other shifts. The near-the-diagonal block is
282 * of size (ns+1) x ns.
283 
284  CALL claset( 'FULL', ns+1, ns+1, czero, cone, qc, ldqc )
285  CALL claset( 'FULL', ns, ns, czero, cone, zc, ldzc )
286 
287  DO i = 1, ns
288 * Introduce the shift
289  scale = sqrt( abs( alpha( i ) ) ) * sqrt( abs( beta( i ) ) )
290  IF( scale .GE. safmin .AND. scale .LE. safmax ) THEN
291  alpha( i ) = alpha( i )/scale
292  beta( i ) = beta( i )/scale
293  END IF
294 
295  temp2 = beta( i )*a( ilo, ilo )-alpha( i )*b( ilo, ilo )
296  temp3 = beta( i )*a( ilo+1, ilo )
297 
298  IF ( abs( temp2 ) .GT. safmax .OR.
299  $ abs( temp3 ) .GT. safmax ) THEN
300  temp2 = cone
301  temp3 = czero
302  END IF
303 
304  CALL clartg( temp2, temp3, c, s, temp )
305  CALL crot( ns, a( ilo, ilo ), lda, a( ilo+1, ilo ), lda, c,
306  $ s )
307  CALL crot( ns, b( ilo, ilo ), ldb, b( ilo+1, ilo ), ldb, c,
308  $ s )
309  CALL crot( ns+1, qc( 1, 1 ), 1, qc( 1, 2 ), 1, c, conjg( s ) )
310 
311 * Chase the shift down
312  DO j = 1, ns-i
313 
314  CALL claqz1( .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 cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
330  $ a( ilo, ilo+ns ), lda, czero, work, sheight )
331  CALL clacpy( 'ALL', sheight, swidth, work, sheight, a( ilo,
332  $ ilo+ns ), lda )
333  CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
334  $ b( ilo, ilo+ns ), ldb, czero, work, sheight )
335  CALL clacpy( 'ALL', sheight, swidth, work, sheight, b( ilo,
336  $ ilo+ns ), ldb )
337  END IF
338  IF ( ilq ) THEN
339  CALL cgemm( 'N', 'N', n, sheight, sheight, cone, q( 1, ilo ),
340  $ ldq, qc, ldqc, czero, work, n )
341  CALL clacpy( 'ALL', n, sheight, work, n, q( 1, ilo ), ldq )
342  END IF
343 
344 * Update A(istartm:ilo-1,ilo:ilo+ns-1) and B(istartm:ilo-1,ilo:ilo+ns-1)
345 * from the right with Zc(1:ns,1:ns)
346  sheight = ilo-1-istartm+1
347  swidth = ns
348  IF ( sheight > 0 ) THEN
349  CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
350  $ a( istartm, ilo ), lda, zc, ldzc, czero, work,
351  $ sheight )
352  CALL clacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
353  $ ilo ), lda )
354  CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
355  $ b( istartm, ilo ), ldb, zc, ldzc, czero, work,
356  $ sheight )
357  CALL clacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
358  $ ilo ), ldb )
359  END IF
360  IF ( ilz ) THEN
361  CALL cgemm( 'N', 'N', n, swidth, swidth, cone, z( 1, ilo ),
362  $ ldz, zc, ldzc, czero, work, n )
363  CALL clacpy( 'ALL', n, swidth, work, n, z( 1, ilo ), ldz )
364  END IF
365 
366 * The following block chases the shifts down to the bottom
367 * right block. If possible, a shift is moved down npos
368 * positions at a time
369 
370  k = ilo
371  DO WHILE ( k < ihi-ns )
372  np = min( ihi-ns-k, npos )
373 * Size of the near-the-diagonal block
374  nblock = ns+np
375 * istartb points to the first row we will be updating
376  istartb = k+1
377 * istopb points to the last column we will be updating
378  istopb = k+nblock-1
379 
380  CALL claset( 'FULL', ns+np, ns+np, czero, cone, qc, ldqc )
381  CALL claset( 'FULL', ns+np, ns+np, czero, cone, zc, ldzc )
382 
383 * Near the diagonal shift chase
384  DO i = ns-1, 0, -1
385  DO j = 0, np-1
386 * Move down the block with index k+i+j, updating
387 * the (ns+np x ns+np) block:
388 * (k:k+ns+np,k:k+ns+np-1)
389  CALL claqz1( .true., .true., k+i+j, istartb, istopb, ihi,
390  $ a, lda, b, ldb, nblock, k+1, qc, ldqc,
391  $ nblock, k, zc, ldzc )
392  END DO
393  END DO
394 
395 * Update rest of the pencil
396 
397 * Update A(k+1:k+ns+np, k+ns+np:istopm) and
398 * B(k+1:k+ns+np, k+ns+np:istopm)
399 * from the left with Qc(1:ns+np,1:ns+np)'
400  sheight = ns+np
401  swidth = istopm-( k+ns+np )+1
402  IF ( swidth > 0 ) THEN
403  CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
404  $ ldqc, a( k+1, k+ns+np ), lda, czero, work,
405  $ sheight )
406  CALL clacpy( 'ALL', sheight, swidth, work, sheight, a( k+1,
407  $ k+ns+np ), lda )
408  CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc,
409  $ ldqc, b( k+1, k+ns+np ), ldb, czero, work,
410  $ sheight )
411  CALL clacpy( 'ALL', sheight, swidth, work, sheight, b( k+1,
412  $ k+ns+np ), ldb )
413  END IF
414  IF ( ilq ) THEN
415  CALL cgemm( 'N', 'N', n, nblock, nblock, cone, q( 1, k+1 ),
416  $ ldq, qc, ldqc, czero, work, n )
417  CALL clacpy( 'ALL', n, nblock, work, n, q( 1, k+1 ), ldq )
418  END IF
419 
420 * Update A(istartm:k,k:k+ns+npos-1) and B(istartm:k,k:k+ns+npos-1)
421 * from the right with Zc(1:ns+np,1:ns+np)
422  sheight = k-istartm+1
423  swidth = nblock
424  IF ( sheight > 0 ) THEN
425  CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
426  $ a( istartm, k ), lda, zc, ldzc, czero, work,
427  $ sheight )
428  CALL clacpy( 'ALL', sheight, swidth, work, sheight,
429  $ a( istartm, k ), lda )
430  CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
431  $ b( istartm, k ), ldb, zc, ldzc, czero, work,
432  $ sheight )
433  CALL clacpy( 'ALL', sheight, swidth, work, sheight,
434  $ b( istartm, k ), ldb )
435  END IF
436  IF ( ilz ) THEN
437  CALL cgemm( 'N', 'N', n, nblock, nblock, cone, z( 1, k ),
438  $ ldz, zc, ldzc, czero, work, n )
439  CALL clacpy( 'ALL', n, nblock, work, n, z( 1, k ), ldz )
440  END IF
441 
442  k = k+np
443 
444  END DO
445 
446 * The following block removes the shifts from the bottom right corner
447 * one by one. Updates are initially applied to A(ihi-ns+1:ihi,ihi-ns:ihi).
448 
449  CALL claset( 'FULL', ns, ns, czero, cone, qc, ldqc )
450  CALL claset( 'FULL', ns+1, ns+1, czero, cone, zc, ldzc )
451 
452 * istartb points to the first row we will be updating
453  istartb = ihi-ns+1
454 * istopb points to the last column we will be updating
455  istopb = ihi
456 
457  DO i = 1, ns
458 * Chase the shift down to the bottom right corner
459  DO ishift = ihi-i, ihi-1
460  CALL claqz1( .true., .true., ishift, istartb, istopb, ihi,
461  $ a, lda, b, ldb, ns, ihi-ns+1, qc, ldqc, ns+1,
462  $ ihi-ns, zc, ldzc )
463  END DO
464 
465  END DO
466 
467 * Update rest of the pencil
468 
469 * Update A(ihi-ns+1:ihi, ihi+1:istopm)
470 * from the left with Qc(1:ns,1:ns)'
471  sheight = ns
472  swidth = istopm-( ihi+1 )+1
473  IF ( swidth > 0 ) THEN
474  CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
475  $ a( ihi-ns+1, ihi+1 ), lda, czero, work, sheight )
476  CALL clacpy( 'ALL', sheight, swidth, work, sheight,
477  $ a( ihi-ns+1, ihi+1 ), lda )
478  CALL cgemm( 'C', 'N', sheight, swidth, sheight, cone, qc, ldqc,
479  $ b( ihi-ns+1, ihi+1 ), ldb, czero, work, sheight )
480  CALL clacpy( 'ALL', sheight, swidth, work, sheight,
481  $ b( ihi-ns+1, ihi+1 ), ldb )
482  END IF
483  IF ( ilq ) THEN
484  CALL cgemm( 'N', 'N', n, ns, ns, cone, q( 1, ihi-ns+1 ), ldq,
485  $ qc, ldqc, czero, work, n )
486  CALL clacpy( 'ALL', n, ns, work, n, q( 1, ihi-ns+1 ), ldq )
487  END IF
488 
489 * Update A(istartm:ihi-ns,ihi-ns:ihi)
490 * from the right with Zc(1:ns+1,1:ns+1)
491  sheight = ihi-ns-istartm+1
492  swidth = ns+1
493  IF ( sheight > 0 ) THEN
494  CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
495  $ a( istartm, ihi-ns ), lda, zc, ldzc, czero, work,
496  $ sheight )
497  CALL clacpy( 'ALL', sheight, swidth, work, sheight, a( istartm,
498  $ ihi-ns ), lda )
499  CALL cgemm( 'N', 'N', sheight, swidth, swidth, cone,
500  $ b( istartm, ihi-ns ), ldb, zc, ldzc, czero, work,
501  $ sheight )
502  CALL clacpy( 'ALL', sheight, swidth, work, sheight, b( istartm,
503  $ ihi-ns ), ldb )
504  END IF
505  IF ( ilz ) THEN
506  CALL cgemm( 'N', 'N', n, ns+1, ns+1, cone, z( 1, ihi-ns ), ldz,
507  $ zc, ldzc, czero, work, n )
508  CALL clacpy( 'ALL', n, ns+1, work, n, z( 1, ihi-ns ), ldz )
509  END IF
510 
511  END SUBROUTINE
clacpy
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
Definition: clacpy.f:103
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
cgemm
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
Definition: cgemm.f:187
claset
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: claset.f:106
crot
subroutine crot(N, CX, INCX, CY, INCY, C, S)
CROT applies a plane rotation with real cosine and complex sine to a pair of complex vectors...
Definition: crot.f:103
slabad
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
claqz3
subroutine claqz3(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)
CLAQZ3
Definition: claqz3.f:207
claqz1
subroutine claqz1(ILQ, ILZ, K, ISTARTM, ISTOPM, IHI, A, LDA, B, LDB, NQ, QSTART, Q, LDQ, NZ, ZSTART, Z, LDZ)
CLAQZ1
Definition: claqz1.f:173