Annotation of rpl/lapack/lapack/ztgsy2.f, revision 1.21
1.12 bertrand 1: *> \brief \b ZTGSY2 solves the generalized Sylvester equation (unblocked algorithm).
1.9 bertrand 2: *
3: * =========== DOCUMENTATION ===========
4: *
1.17 bertrand 5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
1.9 bertrand 7: *
8: *> \htmlonly
1.17 bertrand 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">
1.9 bertrand 15: *> [TXT]</a>
1.17 bertrand 16: *> \endhtmlonly
1.9 bertrand 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 )
1.17 bertrand 24: *
1.9 bertrand 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: * ..
1.17 bertrand 34: *
1.9 bertrand 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)] =
1.20 bertrand 70: *> = sigma_min(Z) using reverse communication with ZLACON.
1.9 bertrand 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
1.20 bertrand 84: *> = 'N': solve the generalized Sylvester equation (1).
1.9 bertrand 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: *
1.17 bertrand 242: *> \author Univ. of Tennessee
243: *> \author Univ. of California Berkeley
244: *> \author Univ. of Colorado Denver
245: *> \author NAG Ltd.
1.9 bertrand 246: *
247: *> \ingroup complex16SYauxiliary
248: *
249: *> \par Contributors:
250: * ==================
251: *>
252: *> Bo Kagstrom and Peter Poromaa, Department of Computing Science,
253: *> Umea University, S-901 87 Umea, Sweden.
254: *
255: * =====================================================================
1.1 bertrand 256: SUBROUTINE ZTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
257: $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
258: $ INFO )
259: *
1.21 ! bertrand 260: * -- LAPACK auxiliary routine --
1.1 bertrand 261: * -- LAPACK is a software package provided by Univ. of Tennessee, --
262: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
263: *
264: * .. Scalar Arguments ..
265: CHARACTER TRANS
266: INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
267: DOUBLE PRECISION RDSCAL, RDSUM, SCALE
268: * ..
269: * .. Array Arguments ..
270: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * ),
271: $ D( LDD, * ), E( LDE, * ), F( LDF, * )
272: * ..
273: *
274: * =====================================================================
275: *
276: * .. Parameters ..
277: DOUBLE PRECISION ZERO, ONE
278: INTEGER LDZ
279: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, LDZ = 2 )
280: * ..
281: * .. Local Scalars ..
282: LOGICAL NOTRAN
283: INTEGER I, IERR, J, K
284: DOUBLE PRECISION SCALOC
285: COMPLEX*16 ALPHA
286: * ..
287: * .. Local Arrays ..
288: INTEGER IPIV( LDZ ), JPIV( LDZ )
289: COMPLEX*16 RHS( LDZ ), Z( LDZ, LDZ )
290: * ..
291: * .. External Functions ..
292: LOGICAL LSAME
293: EXTERNAL LSAME
294: * ..
295: * .. External Subroutines ..
296: EXTERNAL XERBLA, ZAXPY, ZGESC2, ZGETC2, ZLATDF, ZSCAL
297: * ..
298: * .. Intrinsic Functions ..
299: INTRINSIC DCMPLX, DCONJG, MAX
300: * ..
301: * .. Executable Statements ..
302: *
303: * Decode and test input parameters
304: *
305: INFO = 0
306: IERR = 0
307: NOTRAN = LSAME( TRANS, 'N' )
308: IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
309: INFO = -1
310: ELSE IF( NOTRAN ) THEN
311: IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
312: INFO = -2
313: END IF
314: END IF
315: IF( INFO.EQ.0 ) THEN
316: IF( M.LE.0 ) THEN
317: INFO = -3
318: ELSE IF( N.LE.0 ) THEN
319: INFO = -4
320: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
1.15 bertrand 321: INFO = -6
1.1 bertrand 322: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
323: INFO = -8
324: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
325: INFO = -10
326: ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
327: INFO = -12
328: ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
329: INFO = -14
330: ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
331: INFO = -16
332: END IF
333: END IF
334: IF( INFO.NE.0 ) THEN
335: CALL XERBLA( 'ZTGSY2', -INFO )
336: RETURN
337: END IF
338: *
339: IF( NOTRAN ) THEN
340: *
341: * Solve (I, J) - system
342: * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
343: * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
344: * for I = M, M - 1, ..., 1; J = 1, 2, ..., N
345: *
346: SCALE = ONE
347: SCALOC = ONE
348: DO 30 J = 1, N
349: DO 20 I = M, 1, -1
350: *
351: * Build 2 by 2 system
352: *
353: Z( 1, 1 ) = A( I, I )
354: Z( 2, 1 ) = D( I, I )
355: Z( 1, 2 ) = -B( J, J )
356: Z( 2, 2 ) = -E( J, J )
357: *
358: * Set up right hand side(s)
359: *
360: RHS( 1 ) = C( I, J )
361: RHS( 2 ) = F( I, J )
362: *
363: * Solve Z * x = RHS
364: *
365: CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
366: IF( IERR.GT.0 )
367: $ INFO = IERR
368: IF( IJOB.EQ.0 ) THEN
369: CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
370: IF( SCALOC.NE.ONE ) THEN
371: DO 10 K = 1, N
372: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
373: $ C( 1, K ), 1 )
374: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
375: $ F( 1, K ), 1 )
376: 10 CONTINUE
377: SCALE = SCALE*SCALOC
378: END IF
379: ELSE
380: CALL ZLATDF( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL,
381: $ IPIV, JPIV )
382: END IF
383: *
384: * Unpack solution vector(s)
385: *
386: C( I, J ) = RHS( 1 )
387: F( I, J ) = RHS( 2 )
388: *
389: * Substitute R(I, J) and L(I, J) into remaining equation.
390: *
391: IF( I.GT.1 ) THEN
392: ALPHA = -RHS( 1 )
393: CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 )
394: CALL ZAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 )
395: END IF
396: IF( J.LT.N ) THEN
397: CALL ZAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB,
398: $ C( I, J+1 ), LDC )
399: CALL ZAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE,
400: $ F( I, J+1 ), LDF )
401: END IF
402: *
403: 20 CONTINUE
404: 30 CONTINUE
405: ELSE
406: *
407: * Solve transposed (I, J) - system:
1.8 bertrand 408: * A(I, I)**H * R(I, J) + D(I, I)**H * L(J, J) = C(I, J)
1.1 bertrand 409: * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
410: * for I = 1, 2, ..., M, J = N, N - 1, ..., 1
411: *
412: SCALE = ONE
413: SCALOC = ONE
414: DO 80 I = 1, M
415: DO 70 J = N, 1, -1
416: *
1.8 bertrand 417: * Build 2 by 2 system Z**H
1.1 bertrand 418: *
419: Z( 1, 1 ) = DCONJG( A( I, I ) )
420: Z( 2, 1 ) = -DCONJG( B( J, J ) )
421: Z( 1, 2 ) = DCONJG( D( I, I ) )
422: Z( 2, 2 ) = -DCONJG( E( J, J ) )
423: *
424: *
425: * Set up right hand side(s)
426: *
427: RHS( 1 ) = C( I, J )
428: RHS( 2 ) = F( I, J )
429: *
1.8 bertrand 430: * Solve Z**H * x = RHS
1.1 bertrand 431: *
432: CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
433: IF( IERR.GT.0 )
434: $ INFO = IERR
435: CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
436: IF( SCALOC.NE.ONE ) THEN
437: DO 40 K = 1, N
438: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
439: $ 1 )
440: CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
441: $ 1 )
442: 40 CONTINUE
443: SCALE = SCALE*SCALOC
444: END IF
445: *
446: * Unpack solution vector(s)
447: *
448: C( I, J ) = RHS( 1 )
449: F( I, J ) = RHS( 2 )
450: *
451: * Substitute R(I, J) and L(I, J) into remaining equation.
452: *
453: DO 50 K = 1, J - 1
454: F( I, K ) = F( I, K ) + RHS( 1 )*DCONJG( B( K, J ) ) +
455: $ RHS( 2 )*DCONJG( E( K, J ) )
456: 50 CONTINUE
457: DO 60 K = I + 1, M
458: C( K, J ) = C( K, J ) - DCONJG( A( I, K ) )*RHS( 1 ) -
459: $ DCONJG( D( I, K ) )*RHS( 2 )
460: 60 CONTINUE
461: *
462: 70 CONTINUE
463: 80 CONTINUE
464: END IF
465: RETURN
466: *
467: * End of ZTGSY2
468: *
469: 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