1: *> \brief \b ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download ZTGSY2 + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/ztgsy2.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/ztgsy2.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/ztgsy2.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE ZTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
22: * LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
23: * INFO )
24: *
25: * .. Scalar Arguments ..
26: * CHARACTER TRANS
27: * INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
28: * DOUBLE PRECISION RDSCAL, RDSUM, SCALE
29: * ..
30: * .. Array Arguments ..
31: * COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
32: * $ D( LDD, * ), E( LDE, * ), F( LDF, * )
33: * ..
34: *
35: *
36: *> \par Purpose:
37: * =============
38: *>
39: *> \verbatim
40: *>
41: *> ZTGSY2 solves the generalized Sylvester equation
42: *>
43: *> A * R - L * B = scale * C (1)
44: *> D * R - L * E = scale * F
45: *>
46: *> using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
47: *> (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
48: *> N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
49: *> (i.e., (A,D) and (B,E) in generalized Schur form).
50: *>
51: *> The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
52: *> scaling factor chosen to avoid overflow.
53: *>
54: *> In matrix notation solving equation (1) corresponds to solve
55: *> Zx = scale * b, where Z is defined as
56: *>
57: *> Z = [ kron(In, A) -kron(B**H, Im) ] (2)
58: *> [ kron(In, D) -kron(E**H, Im) ],
59: *>
60: *> Ik is the identity matrix of size k and X**H is the conjuguate transpose of X.
61: *> kron(X, Y) is the Kronecker product between the matrices X and Y.
62: *>
63: *> If TRANS = 'C', y in the conjugate transposed system Z**H*y = scale*b
64: *> is solved for, which is equivalent to solve for R and L in
65: *>
66: *> A**H * R + D**H * L = scale * C (3)
67: *> R * B**H + L * E**H = scale * -F
68: *>
69: *> This case is used to compute an estimate of Dif[(A, D), (B, E)] =
70: *> = sigma_min(Z) using reverse communicaton with ZLACON.
71: *>
72: *> ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL
73: *> of an upper bound on the separation between to matrix pairs. Then
74: *> the input (A, D), (B, E) are sub-pencils of two matrix pairs in
75: *> ZTGSYL.
76: *> \endverbatim
77: *
78: * Arguments:
79: * ==========
80: *
81: *> \param[in] TRANS
82: *> \verbatim
83: *> TRANS is CHARACTER*1
84: *> = 'N', solve the generalized Sylvester equation (1).
85: *> = 'T': solve the 'transposed' system (3).
86: *> \endverbatim
87: *>
88: *> \param[in] IJOB
89: *> \verbatim
90: *> IJOB is INTEGER
91: *> Specifies what kind of functionality to be performed.
92: *> =0: solve (1) only.
93: *> =1: A contribution from this subsystem to a Frobenius
94: *> norm-based estimate of the separation between two matrix
95: *> pairs is computed. (look ahead strategy is used).
96: *> =2: A contribution from this subsystem to a Frobenius
97: *> norm-based estimate of the separation between two matrix
98: *> pairs is computed. (DGECON on sub-systems is used.)
99: *> Not referenced if TRANS = 'T'.
100: *> \endverbatim
101: *>
102: *> \param[in] M
103: *> \verbatim
104: *> M is INTEGER
105: *> On entry, M specifies the order of A and D, and the row
106: *> dimension of C, F, R and L.
107: *> \endverbatim
108: *>
109: *> \param[in] N
110: *> \verbatim
111: *> N is INTEGER
112: *> On entry, N specifies the order of B and E, and the column
113: *> dimension of C, F, R and L.
114: *> \endverbatim
115: *>
116: *> \param[in] A
117: *> \verbatim
118: *> A is COMPLEX*16 array, dimension (LDA, M)
119: *> On entry, A contains an upper triangular matrix.
120: *> \endverbatim
121: *>
122: *> \param[in] LDA
123: *> \verbatim
124: *> LDA is INTEGER
125: *> The leading dimension of the matrix A. LDA >= max(1, M).
126: *> \endverbatim
127: *>
128: *> \param[in] B
129: *> \verbatim
130: *> B is COMPLEX*16 array, dimension (LDB, N)
131: *> On entry, B contains an upper triangular matrix.
132: *> \endverbatim
133: *>
134: *> \param[in] LDB
135: *> \verbatim
136: *> LDB is INTEGER
137: *> The leading dimension of the matrix B. LDB >= max(1, N).
138: *> \endverbatim
139: *>
140: *> \param[in,out] C
141: *> \verbatim
142: *> C is COMPLEX*16 array, dimension (LDC, N)
143: *> On entry, C contains the right-hand-side of the first matrix
144: *> equation in (1).
145: *> On exit, if IJOB = 0, C has been overwritten by the solution
146: *> R.
147: *> \endverbatim
148: *>
149: *> \param[in] LDC
150: *> \verbatim
151: *> LDC is INTEGER
152: *> The leading dimension of the matrix C. LDC >= max(1, M).
153: *> \endverbatim
154: *>
155: *> \param[in] D
156: *> \verbatim
157: *> D is COMPLEX*16 array, dimension (LDD, M)
158: *> On entry, D contains an upper triangular matrix.
159: *> \endverbatim
160: *>
161: *> \param[in] LDD
162: *> \verbatim
163: *> LDD is INTEGER
164: *> The leading dimension of the matrix D. LDD >= max(1, M).
165: *> \endverbatim
166: *>
167: *> \param[in] E
168: *> \verbatim
169: *> E is COMPLEX*16 array, dimension (LDE, N)
170: *> On entry, E contains an upper triangular matrix.
171: *> \endverbatim
172: *>
173: *> \param[in] LDE
174: *> \verbatim
175: *> LDE is INTEGER
176: *> The leading dimension of the matrix E. LDE >= max(1, N).
177: *> \endverbatim
178: *>
179: *> \param[in,out] F
180: *> \verbatim
181: *> F is COMPLEX*16 array, dimension (LDF, N)
182: *> On entry, F contains the right-hand-side of the second matrix
183: *> equation in (1).
184: *> On exit, if IJOB = 0, F has been overwritten by the solution
185: *> L.
186: *> \endverbatim
187: *>
188: *> \param[in] LDF
189: *> \verbatim
190: *> LDF is INTEGER
191: *> The leading dimension of the matrix F. LDF >= max(1, M).
192: *> \endverbatim
193: *>
194: *> \param[out] SCALE
195: *> \verbatim
196: *> SCALE is DOUBLE PRECISION
197: *> On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
198: *> R and L (C and F on entry) will hold the solutions to a
199: *> slightly perturbed system but the input matrices A, B, D and
200: *> E have not been changed. If SCALE = 0, R and L will hold the
201: *> solutions to the homogeneous system with C = F = 0.
202: *> Normally, SCALE = 1.
203: *> \endverbatim
204: *>
205: *> \param[in,out] RDSUM
206: *> \verbatim
207: *> RDSUM is DOUBLE PRECISION
208: *> On entry, the sum of squares of computed contributions to
209: *> the Dif-estimate under computation by ZTGSYL, where the
210: *> scaling factor RDSCAL (see below) has been factored out.
211: *> On exit, the corresponding sum of squares updated with the
212: *> contributions from the current sub-system.
213: *> If TRANS = 'T' RDSUM is not touched.
214: *> NOTE: RDSUM only makes sense when ZTGSY2 is called by
215: *> ZTGSYL.
216: *> \endverbatim
217: *>
218: *> \param[in,out] RDSCAL
219: *> \verbatim
220: *> RDSCAL is DOUBLE PRECISION
221: *> On entry, scaling factor used to prevent overflow in RDSUM.
222: *> On exit, RDSCAL is updated w.r.t. the current contributions
223: *> in RDSUM.
224: *> If TRANS = 'T', RDSCAL is not touched.
225: *> NOTE: RDSCAL only makes sense when ZTGSY2 is called by
226: *> ZTGSYL.
227: *> \endverbatim
228: *>
229: *> \param[out] INFO
230: *> \verbatim
231: *> INFO is INTEGER
232: *> On exit, if INFO is set to
233: *> =0: Successful exit
234: *> <0: If INFO = -i, input argument number i is illegal.
235: *> >0: The matrix pairs (A, D) and (B, E) have common or very
236: *> close eigenvalues.
237: *> \endverbatim
238: *
239: * Authors:
240: * ========
241: *
242: *> \author Univ. of Tennessee
243: *> \author Univ. of California Berkeley
244: *> \author Univ. of Colorado Denver
245: *> \author NAG Ltd.
246: *
247: *> \date September 2012
248: *
249: *> \ingroup complex16SYauxiliary
250: *
251: *> \par Contributors:
252: * ==================
253: *>
254: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
255: *> Umea University, S-901 87 Umea, Sweden.
256: *
257: * =====================================================================
258: SUBROUTINE ZTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
259: $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
260: $ INFO )
261: *
262: * -- LAPACK auxiliary routine (version 3.4.2) --
263: * -- LAPACK is a software package provided by Univ. of Tennessee, --
264: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
265: * September 2012
266: *
267: * .. Scalar Arguments ..
268: CHARACTER TRANS
269: INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
270: DOUBLE PRECISION RDSCAL, RDSUM, SCALE
271: * ..
272: * .. Array Arguments ..
273: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
274: $ D( LDD, * ), E( LDE, * ), F( LDF, * )
275: * ..
276: *
277: * =====================================================================
278: *
279: * .. Parameters ..
280: DOUBLE PRECISION ZERO, ONE
281: INTEGER LDZ
282: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, LDZ = 2 )
283: * ..
284: * .. Local Scalars ..
285: LOGICAL NOTRAN
286: INTEGER I, IERR, J, K
287: DOUBLE PRECISION SCALOC
288: COMPLEX*16 ALPHA
289: * ..
290: * .. Local Arrays ..
291: INTEGER IPIV( LDZ ), JPIV( LDZ )
292: COMPLEX*16 RHS( LDZ ), Z( LDZ, LDZ )
293: * ..
294: * .. External Functions ..
295: LOGICAL LSAME
296: EXTERNAL LSAME
297: * ..
298: * .. External Subroutines ..
299: EXTERNAL XERBLA, ZAXPY, ZGESC2, ZGETC2, ZLATDF, ZSCAL
300: * ..
301: * .. Intrinsic Functions ..
302: INTRINSIC DCMPLX, DCONJG, MAX
303: * ..
304: * .. Executable Statements ..
305: *
306: * Decode and test input parameters
307: *
308: INFO = 0
309: IERR = 0
310: NOTRAN = LSAME( TRANS, 'N' )
311: IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
312: INFO = -1
313: ELSE IF( NOTRAN ) THEN
314: IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
315: INFO = -2
316: END IF
317: END IF
318: IF( INFO.EQ.0 ) THEN
319: IF( M.LE.0 ) THEN
320: INFO = -3
321: ELSE IF( N.LE.0 ) THEN
322: INFO = -4
323: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
324: INFO = -5
325: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
326: INFO = -8
327: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
328: INFO = -10
329: ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
330: INFO = -12
331: ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
332: INFO = -14
333: ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
334: INFO = -16
335: END IF
336: END IF
337: IF( INFO.NE.0 ) THEN
338: CALL XERBLA( 'ZTGSY2', -INFO )
339: RETURN
340: END IF
341: *
342: IF( NOTRAN ) THEN
343: *
344: * Solve (I, J) - system
345: * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
346: * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
347: * for I = M, M - 1, ..., 1; J = 1, 2, ..., N
348: *
349: SCALE = ONE
350: SCALOC = ONE
351: DO 30 J = 1, N
352: DO 20 I = M, 1, -1
353: *
354: * Build 2 by 2 system
355: *
356: Z( 1, 1 ) = A( I, I )
357: Z( 2, 1 ) = D( I, I )
358: Z( 1, 2 ) = -B( J, J )
359: Z( 2, 2 ) = -E( J, J )
360: *
361: * Set up right hand side(s)
362: *
363: RHS( 1 ) = C( I, J )
364: RHS( 2 ) = F( I, J )
365: *
366: * Solve Z * x = RHS
367: *
368: CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
369: IF( IERR.GT.0 )
370: $ INFO = IERR
371: IF( IJOB.EQ.0 ) THEN
372: CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
373: IF( SCALOC.NE.ONE ) THEN
374: DO 10 K = 1, N
375: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
376: $ C( 1, K ), 1 )
377: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
378: $ F( 1, K ), 1 )
379: 10 CONTINUE
380: SCALE = SCALE*SCALOC
381: END IF
382: ELSE
383: CALL ZLATDF( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL,
384: $ IPIV, JPIV )
385: END IF
386: *
387: * Unpack solution vector(s)
388: *
389: C( I, J ) = RHS( 1 )
390: F( I, J ) = RHS( 2 )
391: *
392: * Substitute R(I, J) and L(I, J) into remaining equation.
393: *
394: IF( I.GT.1 ) THEN
395: ALPHA = -RHS( 1 )
396: CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 )
397: CALL ZAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 )
398: END IF
399: IF( J.LT.N ) THEN
400: CALL ZAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB,
401: $ C( I, J+1 ), LDC )
402: CALL ZAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE,
403: $ F( I, J+1 ), LDF )
404: END IF
405: *
406: 20 CONTINUE
407: 30 CONTINUE
408: ELSE
409: *
410: * Solve transposed (I, J) - system:
411: * A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
412: * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
413: * for I = 1, 2, ..., M, J = N, N - 1, ..., 1
414: *
415: SCALE = ONE
416: SCALOC = ONE
417: DO 80 I = 1, M
418: DO 70 J = N, 1, -1
419: *
420: * Build 2 by 2 system Z**H
421: *
422: Z( 1, 1 ) = DCONJG( A( I, I ) )
423: Z( 2, 1 ) = -DCONJG( B( J, J ) )
424: Z( 1, 2 ) = DCONJG( D( I, I ) )
425: Z( 2, 2 ) = -DCONJG( E( J, J ) )
426: *
427: *
428: * Set up right hand side(s)
429: *
430: RHS( 1 ) = C( I, J )
431: RHS( 2 ) = F( I, J )
432: *
433: * Solve Z**H * x = RHS
434: *
435: CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
436: IF( IERR.GT.0 )
437: $ INFO = IERR
438: CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
439: IF( SCALOC.NE.ONE ) THEN
440: DO 40 K = 1, N
441: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
442: $ 1 )
443: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
444: $ 1 )
445: 40 CONTINUE
446: SCALE = SCALE*SCALOC
447: END IF
448: *
449: * Unpack solution vector(s)
450: *
451: C( I, J ) = RHS( 1 )
452: F( I, J ) = RHS( 2 )
453: *
454: * Substitute R(I, J) and L(I, J) into remaining equation.
455: *
456: DO 50 K = 1, J - 1
457: F( I, K ) = F( I, K ) + RHS( 1 )*DCONJG( B( K, J ) ) +
458: $ RHS( 2 )*DCONJG( E( K, J ) )
459: 50 CONTINUE
460: DO 60 K = I + 1, M
461: C( K, J ) = C( K, J ) - DCONJG( A( I, K ) )*RHS( 1 ) -
462: $ DCONJG( D( I, K ) )*RHS( 2 )
463: 60 CONTINUE
464: *
465: 70 CONTINUE
466: 80 CONTINUE
467: END IF
468: RETURN
469: *
470: * End of ZTGSY2
471: *
472: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztgsy2.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack