001:       SUBROUTINE ZTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
002:      $                   LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
003:      $                   INFO )
004: *
005: *  -- LAPACK auxiliary routine (version 3.2) --
006: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
007: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
008: *     November 2006
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          TRANS
012:       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N
013:       DOUBLE PRECISION   RDSCAL, RDSUM, SCALE
014: *     ..
015: *     .. Array Arguments ..
016:       COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * ),
017:      $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  ZTGSY2 solves the generalized Sylvester equation
024: *
025: *              A * R - L * B = scale *   C               (1)
026: *              D * R - L * E = scale * F
027: *
028: *  using Level 1 and 2 BLAS, where R and L are unknown M-by-N matrices,
029: *  (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
030: *  N-by-N and M-by-N, respectively. A, B, D and E are upper triangular
031: *  (i.e., (A,D) and (B,E) in generalized Schur form).
032: *
033: *  The solution (R, L) overwrites (C, F). 0 <= SCALE <= 1 is an output
034: *  scaling factor chosen to avoid overflow.
035: *
036: *  In matrix notation solving equation (1) corresponds to solve
037: *  Zx = scale * b, where Z is defined as
038: *
039: *         Z = [ kron(In, A)  -kron(B', Im) ]             (2)
040: *             [ kron(In, D)  -kron(E', Im) ],
041: *
042: *  Ik is the identity matrix of size k and X' is the transpose of X.
043: *  kron(X, Y) is the Kronecker product between the matrices X and Y.
044: *
045: *  If TRANS = 'C', y in the conjugate transposed system Z'y = scale*b
046: *  is solved for, which is equivalent to solve for R and L in
047: *
048: *              A' * R  + D' * L   = scale *  C           (3)
049: *              R  * B' + L  * E'  = scale * -F
050: *
051: *  This case is used to compute an estimate of Dif[(A, D), (B, E)] =
052: *  = sigma_min(Z) using reverse communicaton with ZLACON.
053: *
054: *  ZTGSY2 also (IJOB >= 1) contributes to the computation in ZTGSYL
055: *  of an upper bound on the separation between to matrix pairs. Then
056: *  the input (A, D), (B, E) are sub-pencils of two matrix pairs in
057: *  ZTGSYL.
058: *
059: *  Arguments
060: *  =========
061: *
062: *  TRANS   (input) CHARACTER*1
063: *          = 'N', solve the generalized Sylvester equation (1).
064: *          = 'T': solve the 'transposed' system (3).
065: *
066: *  IJOB    (input) INTEGER
067: *          Specifies what kind of functionality to be performed.
068: *          =0: solve (1) only.
069: *          =1: A contribution from this subsystem to a Frobenius
070: *              norm-based estimate of the separation between two matrix
071: *              pairs is computed. (look ahead strategy is used).
072: *          =2: A contribution from this subsystem to a Frobenius
073: *              norm-based estimate of the separation between two matrix
074: *              pairs is computed. (DGECON on sub-systems is used.)
075: *          Not referenced if TRANS = 'T'.
076: *
077: *  M       (input) INTEGER
078: *          On entry, M specifies the order of A and D, and the row
079: *          dimension of C, F, R and L.
080: *
081: *  N       (input) INTEGER
082: *          On entry, N specifies the order of B and E, and the column
083: *          dimension of C, F, R and L.
084: *
085: *  A       (input) COMPLEX*16 array, dimension (LDA, M)
086: *          On entry, A contains an upper triangular matrix.
087: *
088: *  LDA     (input) INTEGER
089: *          The leading dimension of the matrix A. LDA >= max(1, M).
090: *
091: *  B       (input) COMPLEX*16 array, dimension (LDB, N)
092: *          On entry, B contains an upper triangular matrix.
093: *
094: *  LDB     (input) INTEGER
095: *          The leading dimension of the matrix B. LDB >= max(1, N).
096: *
097: *  C       (input/output) COMPLEX*16 array, dimension (LDC, N)
098: *          On entry, C contains the right-hand-side of the first matrix
099: *          equation in (1).
100: *          On exit, if IJOB = 0, C has been overwritten by the solution
101: *          R.
102: *
103: *  LDC     (input) INTEGER
104: *          The leading dimension of the matrix C. LDC >= max(1, M).
105: *
106: *  D       (input) COMPLEX*16 array, dimension (LDD, M)
107: *          On entry, D contains an upper triangular matrix.
108: *
109: *  LDD     (input) INTEGER
110: *          The leading dimension of the matrix D. LDD >= max(1, M).
111: *
112: *  E       (input) COMPLEX*16 array, dimension (LDE, N)
113: *          On entry, E contains an upper triangular matrix.
114: *
115: *  LDE     (input) INTEGER
116: *          The leading dimension of the matrix E. LDE >= max(1, N).
117: *
118: *  F       (input/output) COMPLEX*16 array, dimension (LDF, N)
119: *          On entry, F contains the right-hand-side of the second matrix
120: *          equation in (1).
121: *          On exit, if IJOB = 0, F has been overwritten by the solution
122: *          L.
123: *
124: *  LDF     (input) INTEGER
125: *          The leading dimension of the matrix F. LDF >= max(1, M).
126: *
127: *  SCALE   (output) DOUBLE PRECISION
128: *          On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
129: *          R and L (C and F on entry) will hold the solutions to a
130: *          slightly perturbed system but the input matrices A, B, D and
131: *          E have not been changed. If SCALE = 0, R and L will hold the
132: *          solutions to the homogeneous system with C = F = 0.
133: *          Normally, SCALE = 1.
134: *
135: *  RDSUM   (input/output) DOUBLE PRECISION
136: *          On entry, the sum of squares of computed contributions to
137: *          the Dif-estimate under computation by ZTGSYL, where the
138: *          scaling factor RDSCAL (see below) has been factored out.
139: *          On exit, the corresponding sum of squares updated with the
140: *          contributions from the current sub-system.
141: *          If TRANS = 'T' RDSUM is not touched.
142: *          NOTE: RDSUM only makes sense when ZTGSY2 is called by
143: *          ZTGSYL.
144: *
145: *  RDSCAL  (input/output) DOUBLE PRECISION
146: *          On entry, scaling factor used to prevent overflow in RDSUM.
147: *          On exit, RDSCAL is updated w.r.t. the current contributions
148: *          in RDSUM.
149: *          If TRANS = 'T', RDSCAL is not touched.
150: *          NOTE: RDSCAL only makes sense when ZTGSY2 is called by
151: *          ZTGSYL.
152: *
153: *  INFO    (output) INTEGER
154: *          On exit, if INFO is set to
155: *            =0: Successful exit
156: *            <0: If INFO = -i, input argument number i is illegal.
157: *            >0: The matrix pairs (A, D) and (B, E) have common or very
158: *                close eigenvalues.
159: *
160: *  Further Details
161: *  ===============
162: *
163: *  Based on contributions by
164: *     Bo Kagstrom and Peter Poromaa, Department of Computing Science,
165: *     Umea University, S-901 87 Umea, Sweden.
166: *
167: *  =====================================================================
168: *
169: *     .. Parameters ..
170:       DOUBLE PRECISION   ZERO, ONE
171:       INTEGER            LDZ
172:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, LDZ = 2 )
173: *     ..
174: *     .. Local Scalars ..
175:       LOGICAL            NOTRAN
176:       INTEGER            I, IERR, J, K
177:       DOUBLE PRECISION   SCALOC
178:       COMPLEX*16         ALPHA
179: *     ..
180: *     .. Local Arrays ..
181:       INTEGER            IPIV( LDZ ), JPIV( LDZ )
182:       COMPLEX*16         RHS( LDZ ), Z( LDZ, LDZ )
183: *     ..
184: *     .. External Functions ..
185:       LOGICAL            LSAME
186:       EXTERNAL           LSAME
187: *     ..
188: *     .. External Subroutines ..
189:       EXTERNAL           XERBLA, ZAXPY, ZGESC2, ZGETC2, ZLATDF, ZSCAL
190: *     ..
191: *     .. Intrinsic Functions ..
192:       INTRINSIC          DCMPLX, DCONJG, MAX
193: *     ..
194: *     .. Executable Statements ..
195: *
196: *     Decode and test input parameters
197: *
198:       INFO = 0
199:       IERR = 0
200:       NOTRAN = LSAME( TRANS, 'N' )
201:       IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'C' ) ) THEN
202:          INFO = -1
203:       ELSE IF( NOTRAN ) THEN
204:          IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
205:             INFO = -2
206:          END IF
207:       END IF
208:       IF( INFO.EQ.0 ) THEN
209:          IF( M.LE.0 ) THEN
210:             INFO = -3
211:          ELSE IF( N.LE.0 ) THEN
212:             INFO = -4
213:          ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
214:             INFO = -5
215:          ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
216:             INFO = -8
217:          ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
218:             INFO = -10
219:          ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
220:             INFO = -12
221:          ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
222:             INFO = -14
223:          ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
224:             INFO = -16
225:          END IF
226:       END IF
227:       IF( INFO.NE.0 ) THEN
228:          CALL XERBLA( 'ZTGSY2', -INFO )
229:          RETURN
230:       END IF
231: *
232:       IF( NOTRAN ) THEN
233: *
234: *        Solve (I, J) - system
235: *           A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
236: *           D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
237: *        for I = M, M - 1, ..., 1; J = 1, 2, ..., N
238: *
239:          SCALE = ONE
240:          SCALOC = ONE
241:          DO 30 J = 1, N
242:             DO 20 I = M, 1, -1
243: *
244: *              Build 2 by 2 system
245: *
246:                Z( 1, 1 ) = A( I, I )
247:                Z( 2, 1 ) = D( I, I )
248:                Z( 1, 2 ) = -B( J, J )
249:                Z( 2, 2 ) = -E( J, J )
250: *
251: *              Set up right hand side(s)
252: *
253:                RHS( 1 ) = C( I, J )
254:                RHS( 2 ) = F( I, J )
255: *
256: *              Solve Z * x = RHS
257: *
258:                CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
259:                IF( IERR.GT.0 )
260:      $            INFO = IERR
261:                IF( IJOB.EQ.0 ) THEN
262:                   CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
263:                   IF( SCALOC.NE.ONE ) THEN
264:                      DO 10 K = 1, N
265:                         CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
266:      $                              C( 1, K ), 1 )
267:                         CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ),
268:      $                              F( 1, K ), 1 )
269:    10                CONTINUE
270:                      SCALE = SCALE*SCALOC
271:                   END IF
272:                ELSE
273:                   CALL ZLATDF( IJOB, LDZ, Z, LDZ, RHS, RDSUM, RDSCAL,
274:      $                         IPIV, JPIV )
275:                END IF
276: *
277: *              Unpack solution vector(s)
278: *
279:                C( I, J ) = RHS( 1 )
280:                F( I, J ) = RHS( 2 )
281: *
282: *              Substitute R(I, J) and L(I, J) into remaining equation.
283: *
284:                IF( I.GT.1 ) THEN
285:                   ALPHA = -RHS( 1 )
286:                   CALL ZAXPY( I-1, ALPHA, A( 1, I ), 1, C( 1, J ), 1 )
287:                   CALL ZAXPY( I-1, ALPHA, D( 1, I ), 1, F( 1, J ), 1 )
288:                END IF
289:                IF( J.LT.N ) THEN
290:                   CALL ZAXPY( N-J, RHS( 2 ), B( J, J+1 ), LDB,
291:      $                        C( I, J+1 ), LDC )
292:                   CALL ZAXPY( N-J, RHS( 2 ), E( J, J+1 ), LDE,
293:      $                        F( I, J+1 ), LDF )
294:                END IF
295: *
296:    20       CONTINUE
297:    30    CONTINUE
298:       ELSE
299: *
300: *        Solve transposed (I, J) - system:
301: *           A(I, I)' * R(I, J) + D(I, I)' * L(J, J) = C(I, J)
302: *           R(I, I) * B(J, J) + L(I, J) * E(J, J)   = -F(I, J)
303: *        for I = 1, 2, ..., M, J = N, N - 1, ..., 1
304: *
305:          SCALE = ONE
306:          SCALOC = ONE
307:          DO 80 I = 1, M
308:             DO 70 J = N, 1, -1
309: *
310: *              Build 2 by 2 system Z'
311: *
312:                Z( 1, 1 ) = DCONJG( A( I, I ) )
313:                Z( 2, 1 ) = -DCONJG( B( J, J ) )
314:                Z( 1, 2 ) = DCONJG( D( I, I ) )
315:                Z( 2, 2 ) = -DCONJG( E( J, J ) )
316: *
317: *
318: *              Set up right hand side(s)
319: *
320:                RHS( 1 ) = C( I, J )
321:                RHS( 2 ) = F( I, J )
322: *
323: *              Solve Z' * x = RHS
324: *
325:                CALL ZGETC2( LDZ, Z, LDZ, IPIV, JPIV, IERR )
326:                IF( IERR.GT.0 )
327:      $            INFO = IERR
328:                CALL ZGESC2( LDZ, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
329:                IF( SCALOC.NE.ONE ) THEN
330:                   DO 40 K = 1, N
331:                      CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), C( 1, K ),
332:      $                           1 )
333:                      CALL ZSCAL( M, DCMPLX( SCALOC, ZERO ), F( 1, K ),
334:      $                           1 )
335:    40             CONTINUE
336:                   SCALE = SCALE*SCALOC
337:                END IF
338: *
339: *              Unpack solution vector(s)
340: *
341:                C( I, J ) = RHS( 1 )
342:                F( I, J ) = RHS( 2 )
343: *
344: *              Substitute R(I, J) and L(I, J) into remaining equation.
345: *
346:                DO 50 K = 1, J - 1
347:                   F( I, K ) = F( I, K ) + RHS( 1 )*DCONJG( B( K, J ) ) +
348:      $                        RHS( 2 )*DCONJG( E( K, J ) )
349:    50          CONTINUE
350:                DO 60 K = I + 1, M
351:                   C( K, J ) = C( K, J ) - DCONJG( A( I, K ) )*RHS( 1 ) -
352:      $                        DCONJG( D( I, K ) )*RHS( 2 )
353:    60          CONTINUE
354: *
355:    70       CONTINUE
356:    80    CONTINUE
357:       END IF
358:       RETURN
359: *
360: *     End of ZTGSY2
361: *
362:       END
363: