1: *> \brief \b ZLAQZ2
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZLAQZ2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ZLAQZ2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ZLAQZ2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ZLAQZ2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW, A, LDA, B,
22: * $ LDB, Q, LDQ, Z, LDZ, NS, ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
23: * $ WORK, LWORK, RWORK, REC, INFO )
24: * IMPLICIT NONE
25: *
26: * Arguments
27: * LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
28: * INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
29: * $ LDQC, LDZC, LWORK, REC
30: *
31: * COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
32: * $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
33: * INTEGER, INTENT( OUT ) :: NS, ND, INFO
34: * COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
35: * DOUBLE PRECISION :: RWORK( * )
36: * ..
37: *
38: *
39: *> \par Purpose:
40: * =============
41: *>
42: *> \verbatim
43: *>
44: *> ZLAQZ2 performs AED
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: *> ILO and IHI mark the rows and columns of (A,B) which
82: *> are to be normalized
83: *> \endverbatim
84: *>
85: *> \param[in] NW
86: *> \verbatim
87: *> NW is INTEGER
88: *> The desired size of the deflation window.
89: *> \endverbatim
90: *>
91: *> \param[in,out] A
92: *> \verbatim
93: *> A is COMPLEX*16 array, dimension (LDA, N)
94: *> \endverbatim
95: *>
96: *> \param[in] LDA
97: *> \verbatim
98: *> LDA is INTEGER
99: *> The leading dimension of the array A. LDA >= max( 1, N ).
100: *> \endverbatim
101: *>
102: *> \param[in,out] B
103: *> \verbatim
104: *> B is COMPLEX*16 array, dimension (LDB, N)
105: *> \endverbatim
106: *>
107: *> \param[in] LDB
108: *> \verbatim
109: *> LDB is INTEGER
110: *> The leading dimension of the array B. LDB >= max( 1, N ).
111: *> \endverbatim
112: *>
113: *> \param[in,out] Q
114: *> \verbatim
115: *> Q is COMPLEX*16 array, dimension (LDQ, N)
116: *> \endverbatim
117: *>
118: *> \param[in] LDQ
119: *> \verbatim
120: *> LDQ is INTEGER
121: *> \endverbatim
122: *>
123: *> \param[in,out] Z
124: *> \verbatim
125: *> Z is COMPLEX*16 array, dimension (LDZ, N)
126: *> \endverbatim
127: *>
128: *> \param[in] LDZ
129: *> \verbatim
130: *> LDZ is INTEGER
131: *> \endverbatim
132: *>
133: *> \param[out] NS
134: *> \verbatim
135: *> NS is INTEGER
136: *> The number of unconverged eigenvalues available to
137: *> use as shifts.
138: *> \endverbatim
139: *>
140: *> \param[out] ND
141: *> \verbatim
142: *> ND is INTEGER
143: *> The number of converged eigenvalues found.
144: *> \endverbatim
145: *>
146: *> \param[out] ALPHA
147: *> \verbatim
148: *> ALPHA is COMPLEX*16 array, dimension (N)
149: *> Each scalar alpha defining an eigenvalue
150: *> of GNEP.
151: *> \endverbatim
152: *>
153: *> \param[out] BETA
154: *> \verbatim
155: *> BETA is COMPLEX*16 array, dimension (N)
156: *> The scalars beta that define the eigenvalues of GNEP.
157: *> Together, the quantities alpha = ALPHA(j) and
158: *> beta = BETA(j) represent the j-th eigenvalue of the matrix
159: *> pair (A,B), in one of the forms lambda = alpha/beta or
160: *> mu = beta/alpha. Since either lambda or mu may overflow,
161: *> they should not, in general, be computed.
162: *> \endverbatim
163: *>
164: *> \param[in,out] QC
165: *> \verbatim
166: *> QC is COMPLEX*16 array, dimension (LDQC, NW)
167: *> \endverbatim
168: *>
169: *> \param[in] LDQC
170: *> \verbatim
171: *> LDQC is INTEGER
172: *> \endverbatim
173: *>
174: *> \param[in,out] ZC
175: *> \verbatim
176: *> ZC is COMPLEX*16 array, dimension (LDZC, NW)
177: *> \endverbatim
178: *>
179: *> \param[in] LDZC
180: *> \verbatim
181: *> LDZ is INTEGER
182: *> \endverbatim
183: *>
184: *> \param[out] WORK
185: *> \verbatim
186: *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
187: *> On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
188: *> \endverbatim
189: *>
190: *> \param[in] LWORK
191: *> \verbatim
192: *> LWORK is INTEGER
193: *> The dimension of the array WORK. LWORK >= max(1,N).
194: *>
195: *> If LWORK = -1, then a workspace query is assumed; the routine
196: *> only calculates the optimal size of the WORK array, returns
197: *> this value as the first entry of the WORK array, and no error
198: *> message related to LWORK is issued by XERBLA.
199: *> \endverbatim
200: *>
201: *> \param[out] RWORK
202: *> \verbatim
203: *> RWORK is DOUBLE PRECISION array, dimension (N)
204: *> \endverbatim
205: *>
206: *> \param[in] REC
207: *> \verbatim
208: *> REC is INTEGER
209: *> REC indicates the current recursion level. Should be set
210: *> to 0 on first call.
211: *> \endverbatim
212: *>
213: *> \param[out] INFO
214: *> \verbatim
215: *> INFO is INTEGER
216: *> = 0: successful exit
217: *> < 0: if INFO = -i, the i-th argument had an illegal value
218: *> \endverbatim
219: *
220: * Authors:
221: * ========
222: *
223: *> \author Thijs Steel, KU Leuven
224: *
225: *> \date May 2020
226: *
227: *> \ingroup complex16GEcomputational
228: *>
229: * =====================================================================
230: RECURSIVE SUBROUTINE ZLAQZ2( ILSCHUR, ILQ, ILZ, N, ILO, IHI, NW,
231: $ A, LDA, B, LDB, Q, LDQ, Z, LDZ, NS,
232: $ ND, ALPHA, BETA, QC, LDQC, ZC, LDZC,
233: $ WORK, LWORK, RWORK, REC, INFO )
234: IMPLICIT NONE
235:
236: * Arguments
237: LOGICAL, INTENT( IN ) :: ILSCHUR, ILQ, ILZ
238: INTEGER, INTENT( IN ) :: N, ILO, IHI, NW, LDA, LDB, LDQ, LDZ,
239: $ LDQC, LDZC, LWORK, REC
240:
241: COMPLEX*16, INTENT( INOUT ) :: A( LDA, * ), B( LDB, * ), Q( LDQ,
242: $ * ), Z( LDZ, * ), ALPHA( * ), BETA( * )
243: INTEGER, INTENT( OUT ) :: NS, ND, INFO
244: COMPLEX*16 :: QC( LDQC, * ), ZC( LDZC, * ), WORK( * )
245: DOUBLE PRECISION :: RWORK( * )
246:
247: * Parameters
248: COMPLEX*16 CZERO, CONE
249: PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ), CONE = ( 1.0D+0,
250: $ 0.0D+0 ) )
251: DOUBLE PRECISION :: ZERO, ONE, HALF
252: PARAMETER( ZERO = 0.0D0, ONE = 1.0D0, HALF = 0.5D0 )
253:
254: * Local Scalars
255: INTEGER :: JW, KWTOP, KWBOT, ISTOPM, ISTARTM, K, K2, ZTGEXC_INFO,
256: $ IFST, ILST, LWORKREQ, QZ_SMALL_INFO
257: DOUBLE PRECISION ::SMLNUM, ULP, SAFMIN, SAFMAX, C1, TEMPR
258: COMPLEX*16 :: S, S1, TEMP
259:
260: * External Functions
261: EXTERNAL :: XERBLA, ZLAQZ0, ZLAQZ1, DLABAD, ZLACPY, ZLASET, ZGEMM,
262: $ ZTGEXC, ZLARTG, ZROT
263: DOUBLE PRECISION, EXTERNAL :: DLAMCH
264:
265: INFO = 0
266:
267: * Set up deflation window
268: JW = MIN( NW, IHI-ILO+1 )
269: KWTOP = IHI-JW+1
270: IF ( KWTOP .EQ. ILO ) THEN
271: S = CZERO
272: ELSE
273: S = A( KWTOP, KWTOP-1 )
274: END IF
275:
276: * Determine required workspace
277: IFST = 1
278: ILST = JW
279: CALL ZLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
280: $ B( KWTOP, KWTOP ), LDB, ALPHA, BETA, QC, LDQC, ZC,
281: $ LDZC, WORK, -1, RWORK, REC+1, QZ_SMALL_INFO )
282: LWORKREQ = INT( WORK( 1 ) )+2*JW**2
283: LWORKREQ = MAX( LWORKREQ, N*NW, 2*NW**2+N )
284: IF ( LWORK .EQ.-1 ) THEN
285: * workspace query, quick return
286: WORK( 1 ) = LWORKREQ
287: RETURN
288: ELSE IF ( LWORK .LT. LWORKREQ ) THEN
289: INFO = -26
290: END IF
291:
292: IF( INFO.NE.0 ) THEN
293: CALL XERBLA( 'ZLAQZ2', -INFO )
294: RETURN
295: END IF
296:
297: * Get machine constants
298: SAFMIN = DLAMCH( 'SAFE MINIMUM' )
299: SAFMAX = ONE/SAFMIN
300: CALL DLABAD( SAFMIN, SAFMAX )
301: ULP = DLAMCH( 'PRECISION' )
302: SMLNUM = SAFMIN*( DBLE( N )/ULP )
303:
304: IF ( IHI .EQ. KWTOP ) THEN
305: * 1 by 1 deflation window, just try a regular deflation
306: ALPHA( KWTOP ) = A( KWTOP, KWTOP )
307: BETA( KWTOP ) = B( KWTOP, KWTOP )
308: NS = 1
309: ND = 0
310: IF ( ABS( S ) .LE. MAX( SMLNUM, ULP*ABS( A( KWTOP,
311: $ KWTOP ) ) ) ) THEN
312: NS = 0
313: ND = 1
314: IF ( KWTOP .GT. ILO ) THEN
315: A( KWTOP, KWTOP-1 ) = CZERO
316: END IF
317: END IF
318: END IF
319:
320:
321: * Store window in case of convergence failure
322: CALL ZLACPY( 'ALL', JW, JW, A( KWTOP, KWTOP ), LDA, WORK, JW )
323: CALL ZLACPY( 'ALL', JW, JW, B( KWTOP, KWTOP ), LDB, WORK( JW**2+
324: $ 1 ), JW )
325:
326: * Transform window to real schur form
327: CALL ZLASET( 'FULL', JW, JW, CZERO, CONE, QC, LDQC )
328: CALL ZLASET( 'FULL', JW, JW, CZERO, CONE, ZC, LDZC )
329: CALL ZLAQZ0( 'S', 'V', 'V', JW, 1, JW, A( KWTOP, KWTOP ), LDA,
330: $ B( KWTOP, KWTOP ), LDB, ALPHA, BETA, QC, LDQC, ZC,
331: $ LDZC, WORK( 2*JW**2+1 ), LWORK-2*JW**2, RWORK,
332: $ REC+1, QZ_SMALL_INFO )
333:
334: IF( QZ_SMALL_INFO .NE. 0 ) THEN
335: * Convergence failure, restore the window and exit
336: ND = 0
337: NS = JW-QZ_SMALL_INFO
338: CALL ZLACPY( 'ALL', JW, JW, WORK, JW, A( KWTOP, KWTOP ), LDA )
339: CALL ZLACPY( 'ALL', JW, JW, WORK( JW**2+1 ), JW, B( KWTOP,
340: $ KWTOP ), LDB )
341: RETURN
342: END IF
343:
344: * Deflation detection loop
345: IF ( KWTOP .EQ. ILO .OR. S .EQ. CZERO ) THEN
346: KWBOT = KWTOP-1
347: ELSE
348: KWBOT = IHI
349: K = 1
350: K2 = 1
351: DO WHILE ( K .LE. JW )
352: * Try to deflate eigenvalue
353: TEMPR = ABS( A( KWBOT, KWBOT ) )
354: IF( TEMPR .EQ. ZERO ) THEN
355: TEMPR = ABS( S )
356: END IF
357: IF ( ( ABS( S*QC( 1, KWBOT-KWTOP+1 ) ) ) .LE. MAX( ULP*
358: $ TEMPR, SMLNUM ) ) THEN
359: * Deflatable
360: KWBOT = KWBOT-1
361: ELSE
362: * Not deflatable, move out of the way
363: IFST = KWBOT-KWTOP+1
364: ILST = K2
365: CALL ZTGEXC( .TRUE., .TRUE., JW, A( KWTOP, KWTOP ),
366: $ LDA, B( KWTOP, KWTOP ), LDB, QC, LDQC,
367: $ ZC, LDZC, IFST, ILST, ZTGEXC_INFO )
368: K2 = K2+1
369: END IF
370:
371: K = K+1
372: END DO
373: END IF
374:
375: * Store eigenvalues
376: ND = IHI-KWBOT
377: NS = JW-ND
378: K = KWTOP
379: DO WHILE ( K .LE. IHI )
380: ALPHA( K ) = A( K, K )
381: BETA( K ) = B( K, K )
382: K = K+1
383: END DO
384:
385: IF ( KWTOP .NE. ILO .AND. S .NE. CZERO ) THEN
386: * Reflect spike back, this will create optimally packed bulges
387: A( KWTOP:KWBOT, KWTOP-1 ) = A( KWTOP, KWTOP-1 ) *DCONJG( QC( 1,
388: $ 1:JW-ND ) )
389: DO K = KWBOT-1, KWTOP, -1
390: CALL ZLARTG( A( K, KWTOP-1 ), A( K+1, KWTOP-1 ), C1, S1,
391: $ TEMP )
392: A( K, KWTOP-1 ) = TEMP
393: A( K+1, KWTOP-1 ) = CZERO
394: K2 = MAX( KWTOP, K-1 )
395: CALL ZROT( IHI-K2+1, A( K, K2 ), LDA, A( K+1, K2 ), LDA, C1,
396: $ S1 )
397: CALL ZROT( IHI-( K-1 )+1, B( K, K-1 ), LDB, B( K+1, K-1 ),
398: $ LDB, C1, S1 )
399: CALL ZROT( JW, QC( 1, K-KWTOP+1 ), 1, QC( 1, K+1-KWTOP+1 ),
400: $ 1, C1, DCONJG( S1 ) )
401: END DO
402:
403: * Chase bulges down
404: ISTARTM = KWTOP
405: ISTOPM = IHI
406: K = KWBOT-1
407: DO WHILE ( K .GE. KWTOP )
408:
409: * Move bulge down and remove it
410: DO K2 = K, KWBOT-1
411: CALL ZLAQZ1( .TRUE., .TRUE., K2, KWTOP, KWTOP+JW-1,
412: $ KWBOT, A, LDA, B, LDB, JW, KWTOP, QC, LDQC,
413: $ JW, KWTOP, ZC, LDZC )
414: END DO
415:
416: K = K-1
417: END DO
418:
419: END IF
420:
421: * Apply Qc and Zc to rest of the matrix
422: IF ( ILSCHUR ) THEN
423: ISTARTM = 1
424: ISTOPM = N
425: ELSE
426: ISTARTM = ILO
427: ISTOPM = IHI
428: END IF
429:
430: IF ( ISTOPM-IHI > 0 ) THEN
431: CALL ZGEMM( 'C', 'N', JW, ISTOPM-IHI, JW, CONE, QC, LDQC,
432: $ A( KWTOP, IHI+1 ), LDA, CZERO, WORK, JW )
433: CALL ZLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, A( KWTOP,
434: $ IHI+1 ), LDA )
435: CALL ZGEMM( 'C', 'N', JW, ISTOPM-IHI, JW, CONE, QC, LDQC,
436: $ B( KWTOP, IHI+1 ), LDB, CZERO, WORK, JW )
437: CALL ZLACPY( 'ALL', JW, ISTOPM-IHI, WORK, JW, B( KWTOP,
438: $ IHI+1 ), LDB )
439: END IF
440: IF ( ILQ ) THEN
441: CALL ZGEMM( 'N', 'N', N, JW, JW, CONE, Q( 1, KWTOP ), LDQ, QC,
442: $ LDQC, CZERO, WORK, N )
443: CALL ZLACPY( 'ALL', N, JW, WORK, N, Q( 1, KWTOP ), LDQ )
444: END IF
445:
446: IF ( KWTOP-1-ISTARTM+1 > 0 ) THEN
447: CALL ZGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, CONE, A( ISTARTM,
448: $ KWTOP ), LDA, ZC, LDZC, CZERO, WORK,
449: $ KWTOP-ISTARTM )
450: CALL ZLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
451: $ A( ISTARTM, KWTOP ), LDA )
452: CALL ZGEMM( 'N', 'N', KWTOP-ISTARTM, JW, JW, CONE, B( ISTARTM,
453: $ KWTOP ), LDB, ZC, LDZC, CZERO, WORK,
454: $ KWTOP-ISTARTM )
455: CALL ZLACPY( 'ALL', KWTOP-ISTARTM, JW, WORK, KWTOP-ISTARTM,
456: $ B( ISTARTM, KWTOP ), LDB )
457: END IF
458: IF ( ILZ ) THEN
459: CALL ZGEMM( 'N', 'N', N, JW, JW, CONE, Z( 1, KWTOP ), LDZ, ZC,
460: $ LDZC, CZERO, WORK, N )
461: CALL ZLACPY( 'ALL', N, JW, WORK, N, Z( 1, KWTOP ), LDZ )
462: END IF
463:
464: END SUBROUTINE
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to zlaqz2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack