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
CVSweb interface <joel.bertrand@systella.fr>