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
515