1 SUBROUTINE DTGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
2 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
3 $ IWORK, PQ, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER TRANS
12 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
13 $ PQ
14 DOUBLE PRECISION RDSCAL, RDSUM, SCALE
15 * ..
16 * .. Array Arguments ..
17 INTEGER IWORK( * )
18 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
19 $ D( LDD, * ), E( LDE, * ), F( LDF, * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DTGSY2 solves the generalized Sylvester equation:
26 *
27 * A * R - L * B = scale * C (1)
28 * D * R - L * E = scale * F,
29 *
30 * using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
31 * (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
32 * N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
33 * must be in generalized Schur canonical form, i.e. A, B are upper
34 * quasi triangular and D, E are upper triangular. The solution (R, L)
35 * overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
36 * chosen to avoid overflow.
37 *
38 * In matrix notation solving equation (1) corresponds to solve
39 * Z*x = scale*b, where Z is defined as
40 *
41 * Z = [ kron(In, A) -kron(B**T, Im) ] (2)
42 * [ kron(In, D) -kron(E**T, Im) ],
43 *
44 * Ik is the identity matrix of size k and X**T is the transpose of X.
45 * kron(X, Y) is the Kronecker product between the matrices X and Y.
46 * In the process of solving (1), we solve a number of such systems
47 * where Dim(In), Dim(In) = 1 or 2.
48 *
49 * If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
50 * which is equivalent to solve for R and L in
51 *
52 * A**T * R + D**T * L = scale * C (3)
53 * R * B**T + L * E**T = scale * -F
54 *
55 * This case is used to compute an estimate of Dif[(A, D), (B, E)] =
56 * sigma_min(Z) using reverse communicaton with DLACON.
57 *
58 * DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
59 * of an upper bound on the separation between to matrix pairs. Then
60 * the input (A, D), (B, E) are sub-pencils of the matrix pair in
61 * DTGSYL. See DTGSYL for details.
62 *
63 * Arguments
64 * =========
65 *
66 * TRANS (input) CHARACTER*1
67 * = 'N', solve the generalized Sylvester equation (1).
68 * = 'T': solve the 'transposed' system (3).
69 *
70 * IJOB (input) INTEGER
71 * Specifies what kind of functionality to be performed.
72 * = 0: solve (1) only.
73 * = 1: A contribution from this subsystem to a Frobenius
74 * norm-based estimate of the separation between two matrix
75 * pairs is computed. (look ahead strategy is used).
76 * = 2: A contribution from this subsystem to a Frobenius
77 * norm-based estimate of the separation between two matrix
78 * pairs is computed. (DGECON on sub-systems is used.)
79 * Not referenced if TRANS = 'T'.
80 *
81 * M (input) INTEGER
82 * On entry, M specifies the order of A and D, and the row
83 * dimension of C, F, R and L.
84 *
85 * N (input) INTEGER
86 * On entry, N specifies the order of B and E, and the column
87 * dimension of C, F, R and L.
88 *
89 * A (input) DOUBLE PRECISION array, dimension (LDA, M)
90 * On entry, A contains an upper quasi triangular matrix.
91 *
92 * LDA (input) INTEGER
93 * The leading dimension of the matrix A. LDA >= max(1, M).
94 *
95 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
96 * On entry, B contains an upper quasi triangular matrix.
97 *
98 * LDB (input) INTEGER
99 * The leading dimension of the matrix B. LDB >= max(1, N).
100 *
101 * C (input/output) DOUBLE PRECISION array, dimension (LDC, N)
102 * On entry, C contains the right-hand-side of the first matrix
103 * equation in (1).
104 * On exit, if IJOB = 0, C has been overwritten by the
105 * solution R.
106 *
107 * LDC (input) INTEGER
108 * The leading dimension of the matrix C. LDC >= max(1, M).
109 *
110 * D (input) DOUBLE PRECISION array, dimension (LDD, M)
111 * On entry, D contains an upper triangular matrix.
112 *
113 * LDD (input) INTEGER
114 * The leading dimension of the matrix D. LDD >= max(1, M).
115 *
116 * E (input) DOUBLE PRECISION array, dimension (LDE, N)
117 * On entry, E contains an upper triangular matrix.
118 *
119 * LDE (input) INTEGER
120 * The leading dimension of the matrix E. LDE >= max(1, N).
121 *
122 * F (input/output) DOUBLE PRECISION array, dimension (LDF, N)
123 * On entry, F contains the right-hand-side of the second matrix
124 * equation in (1).
125 * On exit, if IJOB = 0, F has been overwritten by the
126 * solution L.
127 *
128 * LDF (input) INTEGER
129 * The leading dimension of the matrix F. LDF >= max(1, M).
130 *
131 * SCALE (output) DOUBLE PRECISION
132 * On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
133 * R and L (C and F on entry) will hold the solutions to a
134 * slightly perturbed system but the input matrices A, B, D and
135 * E have not been changed. If SCALE = 0, R and L will hold the
136 * solutions to the homogeneous system with C = F = 0. Normally,
137 * SCALE = 1.
138 *
139 * RDSUM (input/output) DOUBLE PRECISION
140 * On entry, the sum of squares of computed contributions to
141 * the Dif-estimate under computation by DTGSYL, where the
142 * scaling factor RDSCAL (see below) has been factored out.
143 * On exit, the corresponding sum of squares updated with the
144 * contributions from the current sub-system.
145 * If TRANS = 'T' RDSUM is not touched.
146 * NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL.
147 *
148 * RDSCAL (input/output) DOUBLE PRECISION
149 * On entry, scaling factor used to prevent overflow in RDSUM.
150 * On exit, RDSCAL is updated w.r.t. the current contributions
151 * in RDSUM.
152 * If TRANS = 'T', RDSCAL is not touched.
153 * NOTE: RDSCAL only makes sense when DTGSY2 is called by
154 * DTGSYL.
155 *
156 * IWORK (workspace) INTEGER array, dimension (M+N+2)
157 *
158 * PQ (output) INTEGER
159 * On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
160 * 8-by-8) solved by this routine.
161 *
162 * INFO (output) INTEGER
163 * On exit, if INFO is set to
164 * =0: Successful exit
165 * <0: If INFO = -i, the i-th argument had an illegal value.
166 * >0: The matrix pairs (A, D) and (B, E) have common or very
167 * close eigenvalues.
168 *
169 * Further Details
170 * ===============
171 *
172 * Based on contributions by
173 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
174 * Umea University, S-901 87 Umea, Sweden.
175 *
176 * =====================================================================
177 * Replaced various illegal calls to DCOPY by calls to DLASET.
178 * Sven Hammarling, 27/5/02.
179 *
180 * .. Parameters ..
181 INTEGER LDZ
182 PARAMETER ( LDZ = 8 )
183 DOUBLE PRECISION ZERO, ONE
184 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
185 * ..
186 * .. Local Scalars ..
187 LOGICAL NOTRAN
188 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
189 $ K, MB, NB, P, Q, ZDIM
190 DOUBLE PRECISION ALPHA, SCALOC
191 * ..
192 * .. Local Arrays ..
193 INTEGER IPIV( LDZ ), JPIV( LDZ )
194 DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ )
195 * ..
196 * .. External Functions ..
197 LOGICAL LSAME
198 EXTERNAL LSAME
199 * ..
200 * .. External Subroutines ..
201 EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DGER, DGESC2,
202 $ DGETC2, DLASET, DLATDF, DSCAL, XERBLA
203 * ..
204 * .. Intrinsic Functions ..
205 INTRINSIC MAX
206 * ..
207 * .. Executable Statements ..
208 *
209 * Decode and test input parameters
210 *
211 INFO = 0
212 IERR = 0
213 NOTRAN = LSAME( TRANS, 'N' )
214 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
215 INFO = -1
216 ELSE IF( NOTRAN ) THEN
217 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
218 INFO = -2
219 END IF
220 END IF
221 IF( INFO.EQ.0 ) THEN
222 IF( M.LE.0 ) THEN
223 INFO = -3
224 ELSE IF( N.LE.0 ) THEN
225 INFO = -4
226 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
227 INFO = -5
228 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
229 INFO = -8
230 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
231 INFO = -10
232 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
233 INFO = -12
234 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
235 INFO = -14
236 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
237 INFO = -16
238 END IF
239 END IF
240 IF( INFO.NE.0 ) THEN
241 CALL XERBLA( 'DTGSY2', -INFO )
242 RETURN
243 END IF
244 *
245 * Determine block structure of A
246 *
247 PQ = 0
248 P = 0
249 I = 1
250 10 CONTINUE
251 IF( I.GT.M )
252 $ GO TO 20
253 P = P + 1
254 IWORK( P ) = I
255 IF( I.EQ.M )
256 $ GO TO 20
257 IF( A( I+1, I ).NE.ZERO ) THEN
258 I = I + 2
259 ELSE
260 I = I + 1
261 END IF
262 GO TO 10
263 20 CONTINUE
264 IWORK( P+1 ) = M + 1
265 *
266 * Determine block structure of B
267 *
268 Q = P + 1
269 J = 1
270 30 CONTINUE
271 IF( J.GT.N )
272 $ GO TO 40
273 Q = Q + 1
274 IWORK( Q ) = J
275 IF( J.EQ.N )
276 $ GO TO 40
277 IF( B( J+1, J ).NE.ZERO ) THEN
278 J = J + 2
279 ELSE
280 J = J + 1
281 END IF
282 GO TO 30
283 40 CONTINUE
284 IWORK( Q+1 ) = N + 1
285 PQ = P*( Q-P-1 )
286 *
287 IF( NOTRAN ) THEN
288 *
289 * Solve (I, J) - subsystem
290 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
291 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
292 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
293 *
294 SCALE = ONE
295 SCALOC = ONE
296 DO 120 J = P + 2, Q
297 JS = IWORK( J )
298 JSP1 = JS + 1
299 JE = IWORK( J+1 ) - 1
300 NB = JE - JS + 1
301 DO 110 I = P, 1, -1
302 *
303 IS = IWORK( I )
304 ISP1 = IS + 1
305 IE = IWORK( I+1 ) - 1
306 MB = IE - IS + 1
307 ZDIM = MB*NB*2
308 *
309 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
310 *
311 * Build a 2-by-2 system Z * x = RHS
312 *
313 Z( 1, 1 ) = A( IS, IS )
314 Z( 2, 1 ) = D( IS, IS )
315 Z( 1, 2 ) = -B( JS, JS )
316 Z( 2, 2 ) = -E( JS, JS )
317 *
318 * Set up right hand side(s)
319 *
320 RHS( 1 ) = C( IS, JS )
321 RHS( 2 ) = F( IS, JS )
322 *
323 * Solve Z * x = RHS
324 *
325 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
326 IF( IERR.GT.0 )
327 $ INFO = IERR
328 *
329 IF( IJOB.EQ.0 ) THEN
330 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
331 $ SCALOC )
332 IF( SCALOC.NE.ONE ) THEN
333 DO 50 K = 1, N
334 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
335 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
336 50 CONTINUE
337 SCALE = SCALE*SCALOC
338 END IF
339 ELSE
340 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
341 $ RDSCAL, IPIV, JPIV )
342 END IF
343 *
344 * Unpack solution vector(s)
345 *
346 C( IS, JS ) = RHS( 1 )
347 F( IS, JS ) = RHS( 2 )
348 *
349 * Substitute R(I, J) and L(I, J) into remaining
350 * equation.
351 *
352 IF( I.GT.1 ) THEN
353 ALPHA = -RHS( 1 )
354 CALL DAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ),
355 $ 1 )
356 CALL DAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ),
357 $ 1 )
358 END IF
359 IF( J.LT.Q ) THEN
360 CALL DAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB,
361 $ C( IS, JE+1 ), LDC )
362 CALL DAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE,
363 $ F( IS, JE+1 ), LDF )
364 END IF
365 *
366 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
367 *
368 * Build a 4-by-4 system Z * x = RHS
369 *
370 Z( 1, 1 ) = A( IS, IS )
371 Z( 2, 1 ) = ZERO
372 Z( 3, 1 ) = D( IS, IS )
373 Z( 4, 1 ) = ZERO
374 *
375 Z( 1, 2 ) = ZERO
376 Z( 2, 2 ) = A( IS, IS )
377 Z( 3, 2 ) = ZERO
378 Z( 4, 2 ) = D( IS, IS )
379 *
380 Z( 1, 3 ) = -B( JS, JS )
381 Z( 2, 3 ) = -B( JS, JSP1 )
382 Z( 3, 3 ) = -E( JS, JS )
383 Z( 4, 3 ) = -E( JS, JSP1 )
384 *
385 Z( 1, 4 ) = -B( JSP1, JS )
386 Z( 2, 4 ) = -B( JSP1, JSP1 )
387 Z( 3, 4 ) = ZERO
388 Z( 4, 4 ) = -E( JSP1, JSP1 )
389 *
390 * Set up right hand side(s)
391 *
392 RHS( 1 ) = C( IS, JS )
393 RHS( 2 ) = C( IS, JSP1 )
394 RHS( 3 ) = F( IS, JS )
395 RHS( 4 ) = F( IS, JSP1 )
396 *
397 * Solve Z * x = RHS
398 *
399 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
400 IF( IERR.GT.0 )
401 $ INFO = IERR
402 *
403 IF( IJOB.EQ.0 ) THEN
404 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
405 $ SCALOC )
406 IF( SCALOC.NE.ONE ) THEN
407 DO 60 K = 1, N
408 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
409 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
410 60 CONTINUE
411 SCALE = SCALE*SCALOC
412 END IF
413 ELSE
414 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
415 $ RDSCAL, IPIV, JPIV )
416 END IF
417 *
418 * Unpack solution vector(s)
419 *
420 C( IS, JS ) = RHS( 1 )
421 C( IS, JSP1 ) = RHS( 2 )
422 F( IS, JS ) = RHS( 3 )
423 F( IS, JSP1 ) = RHS( 4 )
424 *
425 * Substitute R(I, J) and L(I, J) into remaining
426 * equation.
427 *
428 IF( I.GT.1 ) THEN
429 CALL DGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ),
430 $ 1, C( 1, JS ), LDC )
431 CALL DGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ),
432 $ 1, F( 1, JS ), LDF )
433 END IF
434 IF( J.LT.Q ) THEN
435 CALL DAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB,
436 $ C( IS, JE+1 ), LDC )
437 CALL DAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE,
438 $ F( IS, JE+1 ), LDF )
439 CALL DAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB,
440 $ C( IS, JE+1 ), LDC )
441 CALL DAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE,
442 $ F( IS, JE+1 ), LDF )
443 END IF
444 *
445 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
446 *
447 * Build a 4-by-4 system Z * x = RHS
448 *
449 Z( 1, 1 ) = A( IS, IS )
450 Z( 2, 1 ) = A( ISP1, IS )
451 Z( 3, 1 ) = D( IS, IS )
452 Z( 4, 1 ) = ZERO
453 *
454 Z( 1, 2 ) = A( IS, ISP1 )
455 Z( 2, 2 ) = A( ISP1, ISP1 )
456 Z( 3, 2 ) = D( IS, ISP1 )
457 Z( 4, 2 ) = D( ISP1, ISP1 )
458 *
459 Z( 1, 3 ) = -B( JS, JS )
460 Z( 2, 3 ) = ZERO
461 Z( 3, 3 ) = -E( JS, JS )
462 Z( 4, 3 ) = ZERO
463 *
464 Z( 1, 4 ) = ZERO
465 Z( 2, 4 ) = -B( JS, JS )
466 Z( 3, 4 ) = ZERO
467 Z( 4, 4 ) = -E( JS, JS )
468 *
469 * Set up right hand side(s)
470 *
471 RHS( 1 ) = C( IS, JS )
472 RHS( 2 ) = C( ISP1, JS )
473 RHS( 3 ) = F( IS, JS )
474 RHS( 4 ) = F( ISP1, JS )
475 *
476 * Solve Z * x = RHS
477 *
478 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
479 IF( IERR.GT.0 )
480 $ INFO = IERR
481 IF( IJOB.EQ.0 ) THEN
482 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
483 $ SCALOC )
484 IF( SCALOC.NE.ONE ) THEN
485 DO 70 K = 1, N
486 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
487 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
488 70 CONTINUE
489 SCALE = SCALE*SCALOC
490 END IF
491 ELSE
492 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
493 $ RDSCAL, IPIV, JPIV )
494 END IF
495 *
496 * Unpack solution vector(s)
497 *
498 C( IS, JS ) = RHS( 1 )
499 C( ISP1, JS ) = RHS( 2 )
500 F( IS, JS ) = RHS( 3 )
501 F( ISP1, JS ) = RHS( 4 )
502 *
503 * Substitute R(I, J) and L(I, J) into remaining
504 * equation.
505 *
506 IF( I.GT.1 ) THEN
507 CALL DGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA,
508 $ RHS( 1 ), 1, ONE, C( 1, JS ), 1 )
509 CALL DGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD,
510 $ RHS( 1 ), 1, ONE, F( 1, JS ), 1 )
511 END IF
512 IF( J.LT.Q ) THEN
513 CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
514 $ B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC )
515 CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
516 $ E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF )
517 END IF
518 *
519 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
520 *
521 * Build an 8-by-8 system Z * x = RHS
522 *
523 CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
524 *
525 Z( 1, 1 ) = A( IS, IS )
526 Z( 2, 1 ) = A( ISP1, IS )
527 Z( 5, 1 ) = D( IS, IS )
528 *
529 Z( 1, 2 ) = A( IS, ISP1 )
530 Z( 2, 2 ) = A( ISP1, ISP1 )
531 Z( 5, 2 ) = D( IS, ISP1 )
532 Z( 6, 2 ) = D( ISP1, ISP1 )
533 *
534 Z( 3, 3 ) = A( IS, IS )
535 Z( 4, 3 ) = A( ISP1, IS )
536 Z( 7, 3 ) = D( IS, IS )
537 *
538 Z( 3, 4 ) = A( IS, ISP1 )
539 Z( 4, 4 ) = A( ISP1, ISP1 )
540 Z( 7, 4 ) = D( IS, ISP1 )
541 Z( 8, 4 ) = D( ISP1, ISP1 )
542 *
543 Z( 1, 5 ) = -B( JS, JS )
544 Z( 3, 5 ) = -B( JS, JSP1 )
545 Z( 5, 5 ) = -E( JS, JS )
546 Z( 7, 5 ) = -E( JS, JSP1 )
547 *
548 Z( 2, 6 ) = -B( JS, JS )
549 Z( 4, 6 ) = -B( JS, JSP1 )
550 Z( 6, 6 ) = -E( JS, JS )
551 Z( 8, 6 ) = -E( JS, JSP1 )
552 *
553 Z( 1, 7 ) = -B( JSP1, JS )
554 Z( 3, 7 ) = -B( JSP1, JSP1 )
555 Z( 7, 7 ) = -E( JSP1, JSP1 )
556 *
557 Z( 2, 8 ) = -B( JSP1, JS )
558 Z( 4, 8 ) = -B( JSP1, JSP1 )
559 Z( 8, 8 ) = -E( JSP1, JSP1 )
560 *
561 * Set up right hand side(s)
562 *
563 K = 1
564 II = MB*NB + 1
565 DO 80 JJ = 0, NB - 1
566 CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
567 CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
568 K = K + MB
569 II = II + MB
570 80 CONTINUE
571 *
572 * Solve Z * x = RHS
573 *
574 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
575 IF( IERR.GT.0 )
576 $ INFO = IERR
577 IF( IJOB.EQ.0 ) THEN
578 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
579 $ SCALOC )
580 IF( SCALOC.NE.ONE ) THEN
581 DO 90 K = 1, N
582 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
583 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
584 90 CONTINUE
585 SCALE = SCALE*SCALOC
586 END IF
587 ELSE
588 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
589 $ RDSCAL, IPIV, JPIV )
590 END IF
591 *
592 * Unpack solution vector(s)
593 *
594 K = 1
595 II = MB*NB + 1
596 DO 100 JJ = 0, NB - 1
597 CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
598 CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
599 K = K + MB
600 II = II + MB
601 100 CONTINUE
602 *
603 * Substitute R(I, J) and L(I, J) into remaining
604 * equation.
605 *
606 IF( I.GT.1 ) THEN
607 CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
608 $ A( 1, IS ), LDA, RHS( 1 ), MB, ONE,
609 $ C( 1, JS ), LDC )
610 CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
611 $ D( 1, IS ), LDD, RHS( 1 ), MB, ONE,
612 $ F( 1, JS ), LDF )
613 END IF
614 IF( J.LT.Q ) THEN
615 K = MB*NB + 1
616 CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
617 $ MB, B( JS, JE+1 ), LDB, ONE,
618 $ C( IS, JE+1 ), LDC )
619 CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
620 $ MB, E( JS, JE+1 ), LDE, ONE,
621 $ F( IS, JE+1 ), LDF )
622 END IF
623 *
624 END IF
625 *
626 110 CONTINUE
627 120 CONTINUE
628 ELSE
629 *
630 * Solve (I, J) - subsystem
631 * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
632 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
633 * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
634 *
635 SCALE = ONE
636 SCALOC = ONE
637 DO 200 I = 1, P
638 *
639 IS = IWORK( I )
640 ISP1 = IS + 1
641 IE = IWORK ( I+1 ) - 1
642 MB = IE - IS + 1
643 DO 190 J = Q, P + 2, -1
644 *
645 JS = IWORK( J )
646 JSP1 = JS + 1
647 JE = IWORK( J+1 ) - 1
648 NB = JE - JS + 1
649 ZDIM = MB*NB*2
650 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
651 *
652 * Build a 2-by-2 system Z**T * x = RHS
653 *
654 Z( 1, 1 ) = A( IS, IS )
655 Z( 2, 1 ) = -B( JS, JS )
656 Z( 1, 2 ) = D( IS, IS )
657 Z( 2, 2 ) = -E( JS, JS )
658 *
659 * Set up right hand side(s)
660 *
661 RHS( 1 ) = C( IS, JS )
662 RHS( 2 ) = F( IS, JS )
663 *
664 * Solve Z**T * x = RHS
665 *
666 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
667 IF( IERR.GT.0 )
668 $ INFO = IERR
669 *
670 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
671 IF( SCALOC.NE.ONE ) THEN
672 DO 130 K = 1, N
673 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
674 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
675 130 CONTINUE
676 SCALE = SCALE*SCALOC
677 END IF
678 *
679 * Unpack solution vector(s)
680 *
681 C( IS, JS ) = RHS( 1 )
682 F( IS, JS ) = RHS( 2 )
683 *
684 * Substitute R(I, J) and L(I, J) into remaining
685 * equation.
686 *
687 IF( J.GT.P+2 ) THEN
688 ALPHA = RHS( 1 )
689 CALL DAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
690 $ LDF )
691 ALPHA = RHS( 2 )
692 CALL DAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ),
693 $ LDF )
694 END IF
695 IF( I.LT.P ) THEN
696 ALPHA = -RHS( 1 )
697 CALL DAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
698 $ C( IE+1, JS ), 1 )
699 ALPHA = -RHS( 2 )
700 CALL DAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD,
701 $ C( IE+1, JS ), 1 )
702 END IF
703 *
704 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
705 *
706 * Build a 4-by-4 system Z**T * x = RHS
707 *
708 Z( 1, 1 ) = A( IS, IS )
709 Z( 2, 1 ) = ZERO
710 Z( 3, 1 ) = -B( JS, JS )
711 Z( 4, 1 ) = -B( JSP1, JS )
712 *
713 Z( 1, 2 ) = ZERO
714 Z( 2, 2 ) = A( IS, IS )
715 Z( 3, 2 ) = -B( JS, JSP1 )
716 Z( 4, 2 ) = -B( JSP1, JSP1 )
717 *
718 Z( 1, 3 ) = D( IS, IS )
719 Z( 2, 3 ) = ZERO
720 Z( 3, 3 ) = -E( JS, JS )
721 Z( 4, 3 ) = ZERO
722 *
723 Z( 1, 4 ) = ZERO
724 Z( 2, 4 ) = D( IS, IS )
725 Z( 3, 4 ) = -E( JS, JSP1 )
726 Z( 4, 4 ) = -E( JSP1, JSP1 )
727 *
728 * Set up right hand side(s)
729 *
730 RHS( 1 ) = C( IS, JS )
731 RHS( 2 ) = C( IS, JSP1 )
732 RHS( 3 ) = F( IS, JS )
733 RHS( 4 ) = F( IS, JSP1 )
734 *
735 * Solve Z**T * x = RHS
736 *
737 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
738 IF( IERR.GT.0 )
739 $ INFO = IERR
740 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
741 IF( SCALOC.NE.ONE ) THEN
742 DO 140 K = 1, N
743 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
744 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
745 140 CONTINUE
746 SCALE = SCALE*SCALOC
747 END IF
748 *
749 * Unpack solution vector(s)
750 *
751 C( IS, JS ) = RHS( 1 )
752 C( IS, JSP1 ) = RHS( 2 )
753 F( IS, JS ) = RHS( 3 )
754 F( IS, JSP1 ) = RHS( 4 )
755 *
756 * Substitute R(I, J) and L(I, J) into remaining
757 * equation.
758 *
759 IF( J.GT.P+2 ) THEN
760 CALL DAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
761 $ F( IS, 1 ), LDF )
762 CALL DAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
763 $ F( IS, 1 ), LDF )
764 CALL DAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
765 $ F( IS, 1 ), LDF )
766 CALL DAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
767 $ F( IS, 1 ), LDF )
768 END IF
769 IF( I.LT.P ) THEN
770 CALL DGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,
771 $ RHS( 1 ), 1, C( IE+1, JS ), LDC )
772 CALL DGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD,
773 $ RHS( 3 ), 1, C( IE+1, JS ), LDC )
774 END IF
775 *
776 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
777 *
778 * Build a 4-by-4 system Z**T * x = RHS
779 *
780 Z( 1, 1 ) = A( IS, IS )
781 Z( 2, 1 ) = A( IS, ISP1 )
782 Z( 3, 1 ) = -B( JS, JS )
783 Z( 4, 1 ) = ZERO
784 *
785 Z( 1, 2 ) = A( ISP1, IS )
786 Z( 2, 2 ) = A( ISP1, ISP1 )
787 Z( 3, 2 ) = ZERO
788 Z( 4, 2 ) = -B( JS, JS )
789 *
790 Z( 1, 3 ) = D( IS, IS )
791 Z( 2, 3 ) = D( IS, ISP1 )
792 Z( 3, 3 ) = -E( JS, JS )
793 Z( 4, 3 ) = ZERO
794 *
795 Z( 1, 4 ) = ZERO
796 Z( 2, 4 ) = D( ISP1, ISP1 )
797 Z( 3, 4 ) = ZERO
798 Z( 4, 4 ) = -E( JS, JS )
799 *
800 * Set up right hand side(s)
801 *
802 RHS( 1 ) = C( IS, JS )
803 RHS( 2 ) = C( ISP1, JS )
804 RHS( 3 ) = F( IS, JS )
805 RHS( 4 ) = F( ISP1, JS )
806 *
807 * Solve Z**T * x = RHS
808 *
809 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
810 IF( IERR.GT.0 )
811 $ INFO = IERR
812 *
813 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
814 IF( SCALOC.NE.ONE ) THEN
815 DO 150 K = 1, N
816 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
817 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
818 150 CONTINUE
819 SCALE = SCALE*SCALOC
820 END IF
821 *
822 * Unpack solution vector(s)
823 *
824 C( IS, JS ) = RHS( 1 )
825 C( ISP1, JS ) = RHS( 2 )
826 F( IS, JS ) = RHS( 3 )
827 F( ISP1, JS ) = RHS( 4 )
828 *
829 * Substitute R(I, J) and L(I, J) into remaining
830 * equation.
831 *
832 IF( J.GT.P+2 ) THEN
833 CALL DGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
834 $ 1, F( IS, 1 ), LDF )
835 CALL DGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ),
836 $ 1, F( IS, 1 ), LDF )
837 END IF
838 IF( I.LT.P ) THEN
839 CALL DGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ),
840 $ LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
841 $ 1 )
842 CALL DGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ),
843 $ LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ),
844 $ 1 )
845 END IF
846 *
847 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
848 *
849 * Build an 8-by-8 system Z**T * x = RHS
850 *
851 CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
852 *
853 Z( 1, 1 ) = A( IS, IS )
854 Z( 2, 1 ) = A( IS, ISP1 )
855 Z( 5, 1 ) = -B( JS, JS )
856 Z( 7, 1 ) = -B( JSP1, JS )
857 *
858 Z( 1, 2 ) = A( ISP1, IS )
859 Z( 2, 2 ) = A( ISP1, ISP1 )
860 Z( 6, 2 ) = -B( JS, JS )
861 Z( 8, 2 ) = -B( JSP1, JS )
862 *
863 Z( 3, 3 ) = A( IS, IS )
864 Z( 4, 3 ) = A( IS, ISP1 )
865 Z( 5, 3 ) = -B( JS, JSP1 )
866 Z( 7, 3 ) = -B( JSP1, JSP1 )
867 *
868 Z( 3, 4 ) = A( ISP1, IS )
869 Z( 4, 4 ) = A( ISP1, ISP1 )
870 Z( 6, 4 ) = -B( JS, JSP1 )
871 Z( 8, 4 ) = -B( JSP1, JSP1 )
872 *
873 Z( 1, 5 ) = D( IS, IS )
874 Z( 2, 5 ) = D( IS, ISP1 )
875 Z( 5, 5 ) = -E( JS, JS )
876 *
877 Z( 2, 6 ) = D( ISP1, ISP1 )
878 Z( 6, 6 ) = -E( JS, JS )
879 *
880 Z( 3, 7 ) = D( IS, IS )
881 Z( 4, 7 ) = D( IS, ISP1 )
882 Z( 5, 7 ) = -E( JS, JSP1 )
883 Z( 7, 7 ) = -E( JSP1, JSP1 )
884 *
885 Z( 4, 8 ) = D( ISP1, ISP1 )
886 Z( 6, 8 ) = -E( JS, JSP1 )
887 Z( 8, 8 ) = -E( JSP1, JSP1 )
888 *
889 * Set up right hand side(s)
890 *
891 K = 1
892 II = MB*NB + 1
893 DO 160 JJ = 0, NB - 1
894 CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
895 CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
896 K = K + MB
897 II = II + MB
898 160 CONTINUE
899 *
900 *
901 * Solve Z**T * x = RHS
902 *
903 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
904 IF( IERR.GT.0 )
905 $ INFO = IERR
906 *
907 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
908 IF( SCALOC.NE.ONE ) THEN
909 DO 170 K = 1, N
910 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
911 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
912 170 CONTINUE
913 SCALE = SCALE*SCALOC
914 END IF
915 *
916 * Unpack solution vector(s)
917 *
918 K = 1
919 II = MB*NB + 1
920 DO 180 JJ = 0, NB - 1
921 CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
922 CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
923 K = K + MB
924 II = II + MB
925 180 CONTINUE
926 *
927 * Substitute R(I, J) and L(I, J) into remaining
928 * equation.
929 *
930 IF( J.GT.P+2 ) THEN
931 CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
932 $ C( IS, JS ), LDC, B( 1, JS ), LDB, ONE,
933 $ F( IS, 1 ), LDF )
934 CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
935 $ F( IS, JS ), LDF, E( 1, JS ), LDE, ONE,
936 $ F( IS, 1 ), LDF )
937 END IF
938 IF( I.LT.P ) THEN
939 CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
940 $ A( IS, IE+1 ), LDA, C( IS, JS ), LDC,
941 $ ONE, C( IE+1, JS ), LDC )
942 CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
943 $ D( IS, IE+1 ), LDD, F( IS, JS ), LDF,
944 $ ONE, C( IE+1, JS ), LDC )
945 END IF
946 *
947 END IF
948 *
949 190 CONTINUE
950 200 CONTINUE
951 *
952 END IF
953 RETURN
954 *
955 * End of DTGSY2
956 *
957 END
2 $ LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
3 $ IWORK, PQ, INFO )
4 *
5 * -- LAPACK auxiliary routine (version 3.3.1) --
6 * -- LAPACK is a software package provided by Univ. of Tennessee, --
7 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
8 * -- April 2011 --
9 *
10 * .. Scalar Arguments ..
11 CHARACTER TRANS
12 INTEGER IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
13 $ PQ
14 DOUBLE PRECISION RDSCAL, RDSUM, SCALE
15 * ..
16 * .. Array Arguments ..
17 INTEGER IWORK( * )
18 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ),
19 $ D( LDD, * ), E( LDE, * ), F( LDF, * )
20 * ..
21 *
22 * Purpose
23 * =======
24 *
25 * DTGSY2 solves the generalized Sylvester equation:
26 *
27 * A * R - L * B = scale * C (1)
28 * D * R - L * E = scale * F,
29 *
30 * using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
31 * (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
32 * N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
33 * must be in generalized Schur canonical form, i.e. A, B are upper
34 * quasi triangular and D, E are upper triangular. The solution (R, L)
35 * overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
36 * chosen to avoid overflow.
37 *
38 * In matrix notation solving equation (1) corresponds to solve
39 * Z*x = scale*b, where Z is defined as
40 *
41 * Z = [ kron(In, A) -kron(B**T, Im) ] (2)
42 * [ kron(In, D) -kron(E**T, Im) ],
43 *
44 * Ik is the identity matrix of size k and X**T is the transpose of X.
45 * kron(X, Y) is the Kronecker product between the matrices X and Y.
46 * In the process of solving (1), we solve a number of such systems
47 * where Dim(In), Dim(In) = 1 or 2.
48 *
49 * If TRANS = 'T', solve the transposed system Z**T*y = scale*b for y,
50 * which is equivalent to solve for R and L in
51 *
52 * A**T * R + D**T * L = scale * C (3)
53 * R * B**T + L * E**T = scale * -F
54 *
55 * This case is used to compute an estimate of Dif[(A, D), (B, E)] =
56 * sigma_min(Z) using reverse communicaton with DLACON.
57 *
58 * DTGSY2 also (IJOB >= 1) contributes to the computation in DTGSYL
59 * of an upper bound on the separation between to matrix pairs. Then
60 * the input (A, D), (B, E) are sub-pencils of the matrix pair in
61 * DTGSYL. See DTGSYL for details.
62 *
63 * Arguments
64 * =========
65 *
66 * TRANS (input) CHARACTER*1
67 * = 'N', solve the generalized Sylvester equation (1).
68 * = 'T': solve the 'transposed' system (3).
69 *
70 * IJOB (input) INTEGER
71 * Specifies what kind of functionality to be performed.
72 * = 0: solve (1) only.
73 * = 1: A contribution from this subsystem to a Frobenius
74 * norm-based estimate of the separation between two matrix
75 * pairs is computed. (look ahead strategy is used).
76 * = 2: A contribution from this subsystem to a Frobenius
77 * norm-based estimate of the separation between two matrix
78 * pairs is computed. (DGECON on sub-systems is used.)
79 * Not referenced if TRANS = 'T'.
80 *
81 * M (input) INTEGER
82 * On entry, M specifies the order of A and D, and the row
83 * dimension of C, F, R and L.
84 *
85 * N (input) INTEGER
86 * On entry, N specifies the order of B and E, and the column
87 * dimension of C, F, R and L.
88 *
89 * A (input) DOUBLE PRECISION array, dimension (LDA, M)
90 * On entry, A contains an upper quasi triangular matrix.
91 *
92 * LDA (input) INTEGER
93 * The leading dimension of the matrix A. LDA >= max(1, M).
94 *
95 * B (input) DOUBLE PRECISION array, dimension (LDB, N)
96 * On entry, B contains an upper quasi triangular matrix.
97 *
98 * LDB (input) INTEGER
99 * The leading dimension of the matrix B. LDB >= max(1, N).
100 *
101 * C (input/output) DOUBLE PRECISION array, dimension (LDC, N)
102 * On entry, C contains the right-hand-side of the first matrix
103 * equation in (1).
104 * On exit, if IJOB = 0, C has been overwritten by the
105 * solution R.
106 *
107 * LDC (input) INTEGER
108 * The leading dimension of the matrix C. LDC >= max(1, M).
109 *
110 * D (input) DOUBLE PRECISION array, dimension (LDD, M)
111 * On entry, D contains an upper triangular matrix.
112 *
113 * LDD (input) INTEGER
114 * The leading dimension of the matrix D. LDD >= max(1, M).
115 *
116 * E (input) DOUBLE PRECISION array, dimension (LDE, N)
117 * On entry, E contains an upper triangular matrix.
118 *
119 * LDE (input) INTEGER
120 * The leading dimension of the matrix E. LDE >= max(1, N).
121 *
122 * F (input/output) DOUBLE PRECISION array, dimension (LDF, N)
123 * On entry, F contains the right-hand-side of the second matrix
124 * equation in (1).
125 * On exit, if IJOB = 0, F has been overwritten by the
126 * solution L.
127 *
128 * LDF (input) INTEGER
129 * The leading dimension of the matrix F. LDF >= max(1, M).
130 *
131 * SCALE (output) DOUBLE PRECISION
132 * On exit, 0 <= SCALE <= 1. If 0 < SCALE < 1, the solutions
133 * R and L (C and F on entry) will hold the solutions to a
134 * slightly perturbed system but the input matrices A, B, D and
135 * E have not been changed. If SCALE = 0, R and L will hold the
136 * solutions to the homogeneous system with C = F = 0. Normally,
137 * SCALE = 1.
138 *
139 * RDSUM (input/output) DOUBLE PRECISION
140 * On entry, the sum of squares of computed contributions to
141 * the Dif-estimate under computation by DTGSYL, where the
142 * scaling factor RDSCAL (see below) has been factored out.
143 * On exit, the corresponding sum of squares updated with the
144 * contributions from the current sub-system.
145 * If TRANS = 'T' RDSUM is not touched.
146 * NOTE: RDSUM only makes sense when DTGSY2 is called by DTGSYL.
147 *
148 * RDSCAL (input/output) DOUBLE PRECISION
149 * On entry, scaling factor used to prevent overflow in RDSUM.
150 * On exit, RDSCAL is updated w.r.t. the current contributions
151 * in RDSUM.
152 * If TRANS = 'T', RDSCAL is not touched.
153 * NOTE: RDSCAL only makes sense when DTGSY2 is called by
154 * DTGSYL.
155 *
156 * IWORK (workspace) INTEGER array, dimension (M+N+2)
157 *
158 * PQ (output) INTEGER
159 * On exit, the number of subsystems (of size 2-by-2, 4-by-4 and
160 * 8-by-8) solved by this routine.
161 *
162 * INFO (output) INTEGER
163 * On exit, if INFO is set to
164 * =0: Successful exit
165 * <0: If INFO = -i, the i-th argument had an illegal value.
166 * >0: The matrix pairs (A, D) and (B, E) have common or very
167 * close eigenvalues.
168 *
169 * Further Details
170 * ===============
171 *
172 * Based on contributions by
173 * Bo Kagstrom and Peter Poromaa, Department of Computing Science,
174 * Umea University, S-901 87 Umea, Sweden.
175 *
176 * =====================================================================
177 * Replaced various illegal calls to DCOPY by calls to DLASET.
178 * Sven Hammarling, 27/5/02.
179 *
180 * .. Parameters ..
181 INTEGER LDZ
182 PARAMETER ( LDZ = 8 )
183 DOUBLE PRECISION ZERO, ONE
184 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
185 * ..
186 * .. Local Scalars ..
187 LOGICAL NOTRAN
188 INTEGER I, IE, IERR, II, IS, ISP1, J, JE, JJ, JS, JSP1,
189 $ K, MB, NB, P, Q, ZDIM
190 DOUBLE PRECISION ALPHA, SCALOC
191 * ..
192 * .. Local Arrays ..
193 INTEGER IPIV( LDZ ), JPIV( LDZ )
194 DOUBLE PRECISION RHS( LDZ ), Z( LDZ, LDZ )
195 * ..
196 * .. External Functions ..
197 LOGICAL LSAME
198 EXTERNAL LSAME
199 * ..
200 * .. External Subroutines ..
201 EXTERNAL DAXPY, DCOPY, DGEMM, DGEMV, DGER, DGESC2,
202 $ DGETC2, DLASET, DLATDF, DSCAL, XERBLA
203 * ..
204 * .. Intrinsic Functions ..
205 INTRINSIC MAX
206 * ..
207 * .. Executable Statements ..
208 *
209 * Decode and test input parameters
210 *
211 INFO = 0
212 IERR = 0
213 NOTRAN = LSAME( TRANS, 'N' )
214 IF( .NOT.NOTRAN .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
215 INFO = -1
216 ELSE IF( NOTRAN ) THEN
217 IF( ( IJOB.LT.0 ) .OR. ( IJOB.GT.2 ) ) THEN
218 INFO = -2
219 END IF
220 END IF
221 IF( INFO.EQ.0 ) THEN
222 IF( M.LE.0 ) THEN
223 INFO = -3
224 ELSE IF( N.LE.0 ) THEN
225 INFO = -4
226 ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
227 INFO = -5
228 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
229 INFO = -8
230 ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
231 INFO = -10
232 ELSE IF( LDD.LT.MAX( 1, M ) ) THEN
233 INFO = -12
234 ELSE IF( LDE.LT.MAX( 1, N ) ) THEN
235 INFO = -14
236 ELSE IF( LDF.LT.MAX( 1, M ) ) THEN
237 INFO = -16
238 END IF
239 END IF
240 IF( INFO.NE.0 ) THEN
241 CALL XERBLA( 'DTGSY2', -INFO )
242 RETURN
243 END IF
244 *
245 * Determine block structure of A
246 *
247 PQ = 0
248 P = 0
249 I = 1
250 10 CONTINUE
251 IF( I.GT.M )
252 $ GO TO 20
253 P = P + 1
254 IWORK( P ) = I
255 IF( I.EQ.M )
256 $ GO TO 20
257 IF( A( I+1, I ).NE.ZERO ) THEN
258 I = I + 2
259 ELSE
260 I = I + 1
261 END IF
262 GO TO 10
263 20 CONTINUE
264 IWORK( P+1 ) = M + 1
265 *
266 * Determine block structure of B
267 *
268 Q = P + 1
269 J = 1
270 30 CONTINUE
271 IF( J.GT.N )
272 $ GO TO 40
273 Q = Q + 1
274 IWORK( Q ) = J
275 IF( J.EQ.N )
276 $ GO TO 40
277 IF( B( J+1, J ).NE.ZERO ) THEN
278 J = J + 2
279 ELSE
280 J = J + 1
281 END IF
282 GO TO 30
283 40 CONTINUE
284 IWORK( Q+1 ) = N + 1
285 PQ = P*( Q-P-1 )
286 *
287 IF( NOTRAN ) THEN
288 *
289 * Solve (I, J) - subsystem
290 * A(I, I) * R(I, J) - L(I, J) * B(J, J) = C(I, J)
291 * D(I, I) * R(I, J) - L(I, J) * E(J, J) = F(I, J)
292 * for I = P, P - 1, ..., 1; J = 1, 2, ..., Q
293 *
294 SCALE = ONE
295 SCALOC = ONE
296 DO 120 J = P + 2, Q
297 JS = IWORK( J )
298 JSP1 = JS + 1
299 JE = IWORK( J+1 ) - 1
300 NB = JE - JS + 1
301 DO 110 I = P, 1, -1
302 *
303 IS = IWORK( I )
304 ISP1 = IS + 1
305 IE = IWORK( I+1 ) - 1
306 MB = IE - IS + 1
307 ZDIM = MB*NB*2
308 *
309 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
310 *
311 * Build a 2-by-2 system Z * x = RHS
312 *
313 Z( 1, 1 ) = A( IS, IS )
314 Z( 2, 1 ) = D( IS, IS )
315 Z( 1, 2 ) = -B( JS, JS )
316 Z( 2, 2 ) = -E( JS, JS )
317 *
318 * Set up right hand side(s)
319 *
320 RHS( 1 ) = C( IS, JS )
321 RHS( 2 ) = F( IS, JS )
322 *
323 * Solve Z * x = RHS
324 *
325 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
326 IF( IERR.GT.0 )
327 $ INFO = IERR
328 *
329 IF( IJOB.EQ.0 ) THEN
330 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
331 $ SCALOC )
332 IF( SCALOC.NE.ONE ) THEN
333 DO 50 K = 1, N
334 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
335 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
336 50 CONTINUE
337 SCALE = SCALE*SCALOC
338 END IF
339 ELSE
340 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
341 $ RDSCAL, IPIV, JPIV )
342 END IF
343 *
344 * Unpack solution vector(s)
345 *
346 C( IS, JS ) = RHS( 1 )
347 F( IS, JS ) = RHS( 2 )
348 *
349 * Substitute R(I, J) and L(I, J) into remaining
350 * equation.
351 *
352 IF( I.GT.1 ) THEN
353 ALPHA = -RHS( 1 )
354 CALL DAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ),
355 $ 1 )
356 CALL DAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ),
357 $ 1 )
358 END IF
359 IF( J.LT.Q ) THEN
360 CALL DAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB,
361 $ C( IS, JE+1 ), LDC )
362 CALL DAXPY( N-JE, RHS( 2 ), E( JS, JE+1 ), LDE,
363 $ F( IS, JE+1 ), LDF )
364 END IF
365 *
366 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
367 *
368 * Build a 4-by-4 system Z * x = RHS
369 *
370 Z( 1, 1 ) = A( IS, IS )
371 Z( 2, 1 ) = ZERO
372 Z( 3, 1 ) = D( IS, IS )
373 Z( 4, 1 ) = ZERO
374 *
375 Z( 1, 2 ) = ZERO
376 Z( 2, 2 ) = A( IS, IS )
377 Z( 3, 2 ) = ZERO
378 Z( 4, 2 ) = D( IS, IS )
379 *
380 Z( 1, 3 ) = -B( JS, JS )
381 Z( 2, 3 ) = -B( JS, JSP1 )
382 Z( 3, 3 ) = -E( JS, JS )
383 Z( 4, 3 ) = -E( JS, JSP1 )
384 *
385 Z( 1, 4 ) = -B( JSP1, JS )
386 Z( 2, 4 ) = -B( JSP1, JSP1 )
387 Z( 3, 4 ) = ZERO
388 Z( 4, 4 ) = -E( JSP1, JSP1 )
389 *
390 * Set up right hand side(s)
391 *
392 RHS( 1 ) = C( IS, JS )
393 RHS( 2 ) = C( IS, JSP1 )
394 RHS( 3 ) = F( IS, JS )
395 RHS( 4 ) = F( IS, JSP1 )
396 *
397 * Solve Z * x = RHS
398 *
399 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
400 IF( IERR.GT.0 )
401 $ INFO = IERR
402 *
403 IF( IJOB.EQ.0 ) THEN
404 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
405 $ SCALOC )
406 IF( SCALOC.NE.ONE ) THEN
407 DO 60 K = 1, N
408 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
409 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
410 60 CONTINUE
411 SCALE = SCALE*SCALOC
412 END IF
413 ELSE
414 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
415 $ RDSCAL, IPIV, JPIV )
416 END IF
417 *
418 * Unpack solution vector(s)
419 *
420 C( IS, JS ) = RHS( 1 )
421 C( IS, JSP1 ) = RHS( 2 )
422 F( IS, JS ) = RHS( 3 )
423 F( IS, JSP1 ) = RHS( 4 )
424 *
425 * Substitute R(I, J) and L(I, J) into remaining
426 * equation.
427 *
428 IF( I.GT.1 ) THEN
429 CALL DGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ),
430 $ 1, C( 1, JS ), LDC )
431 CALL DGER( IS-1, NB, -ONE, D( 1, IS ), 1, RHS( 1 ),
432 $ 1, F( 1, JS ), LDF )
433 END IF
434 IF( J.LT.Q ) THEN
435 CALL DAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB,
436 $ C( IS, JE+1 ), LDC )
437 CALL DAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE,
438 $ F( IS, JE+1 ), LDF )
439 CALL DAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB,
440 $ C( IS, JE+1 ), LDC )
441 CALL DAXPY( N-JE, RHS( 4 ), E( JSP1, JE+1 ), LDE,
442 $ F( IS, JE+1 ), LDF )
443 END IF
444 *
445 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
446 *
447 * Build a 4-by-4 system Z * x = RHS
448 *
449 Z( 1, 1 ) = A( IS, IS )
450 Z( 2, 1 ) = A( ISP1, IS )
451 Z( 3, 1 ) = D( IS, IS )
452 Z( 4, 1 ) = ZERO
453 *
454 Z( 1, 2 ) = A( IS, ISP1 )
455 Z( 2, 2 ) = A( ISP1, ISP1 )
456 Z( 3, 2 ) = D( IS, ISP1 )
457 Z( 4, 2 ) = D( ISP1, ISP1 )
458 *
459 Z( 1, 3 ) = -B( JS, JS )
460 Z( 2, 3 ) = ZERO
461 Z( 3, 3 ) = -E( JS, JS )
462 Z( 4, 3 ) = ZERO
463 *
464 Z( 1, 4 ) = ZERO
465 Z( 2, 4 ) = -B( JS, JS )
466 Z( 3, 4 ) = ZERO
467 Z( 4, 4 ) = -E( JS, JS )
468 *
469 * Set up right hand side(s)
470 *
471 RHS( 1 ) = C( IS, JS )
472 RHS( 2 ) = C( ISP1, JS )
473 RHS( 3 ) = F( IS, JS )
474 RHS( 4 ) = F( ISP1, JS )
475 *
476 * Solve Z * x = RHS
477 *
478 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
479 IF( IERR.GT.0 )
480 $ INFO = IERR
481 IF( IJOB.EQ.0 ) THEN
482 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
483 $ SCALOC )
484 IF( SCALOC.NE.ONE ) THEN
485 DO 70 K = 1, N
486 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
487 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
488 70 CONTINUE
489 SCALE = SCALE*SCALOC
490 END IF
491 ELSE
492 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
493 $ RDSCAL, IPIV, JPIV )
494 END IF
495 *
496 * Unpack solution vector(s)
497 *
498 C( IS, JS ) = RHS( 1 )
499 C( ISP1, JS ) = RHS( 2 )
500 F( IS, JS ) = RHS( 3 )
501 F( ISP1, JS ) = RHS( 4 )
502 *
503 * Substitute R(I, J) and L(I, J) into remaining
504 * equation.
505 *
506 IF( I.GT.1 ) THEN
507 CALL DGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA,
508 $ RHS( 1 ), 1, ONE, C( 1, JS ), 1 )
509 CALL DGEMV( 'N', IS-1, MB, -ONE, D( 1, IS ), LDD,
510 $ RHS( 1 ), 1, ONE, F( 1, JS ), 1 )
511 END IF
512 IF( J.LT.Q ) THEN
513 CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
514 $ B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC )
515 CALL DGER( MB, N-JE, ONE, RHS( 3 ), 1,
516 $ E( JS, JE+1 ), LDE, F( IS, JE+1 ), LDF )
517 END IF
518 *
519 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
520 *
521 * Build an 8-by-8 system Z * x = RHS
522 *
523 CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
524 *
525 Z( 1, 1 ) = A( IS, IS )
526 Z( 2, 1 ) = A( ISP1, IS )
527 Z( 5, 1 ) = D( IS, IS )
528 *
529 Z( 1, 2 ) = A( IS, ISP1 )
530 Z( 2, 2 ) = A( ISP1, ISP1 )
531 Z( 5, 2 ) = D( IS, ISP1 )
532 Z( 6, 2 ) = D( ISP1, ISP1 )
533 *
534 Z( 3, 3 ) = A( IS, IS )
535 Z( 4, 3 ) = A( ISP1, IS )
536 Z( 7, 3 ) = D( IS, IS )
537 *
538 Z( 3, 4 ) = A( IS, ISP1 )
539 Z( 4, 4 ) = A( ISP1, ISP1 )
540 Z( 7, 4 ) = D( IS, ISP1 )
541 Z( 8, 4 ) = D( ISP1, ISP1 )
542 *
543 Z( 1, 5 ) = -B( JS, JS )
544 Z( 3, 5 ) = -B( JS, JSP1 )
545 Z( 5, 5 ) = -E( JS, JS )
546 Z( 7, 5 ) = -E( JS, JSP1 )
547 *
548 Z( 2, 6 ) = -B( JS, JS )
549 Z( 4, 6 ) = -B( JS, JSP1 )
550 Z( 6, 6 ) = -E( JS, JS )
551 Z( 8, 6 ) = -E( JS, JSP1 )
552 *
553 Z( 1, 7 ) = -B( JSP1, JS )
554 Z( 3, 7 ) = -B( JSP1, JSP1 )
555 Z( 7, 7 ) = -E( JSP1, JSP1 )
556 *
557 Z( 2, 8 ) = -B( JSP1, JS )
558 Z( 4, 8 ) = -B( JSP1, JSP1 )
559 Z( 8, 8 ) = -E( JSP1, JSP1 )
560 *
561 * Set up right hand side(s)
562 *
563 K = 1
564 II = MB*NB + 1
565 DO 80 JJ = 0, NB - 1
566 CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
567 CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
568 K = K + MB
569 II = II + MB
570 80 CONTINUE
571 *
572 * Solve Z * x = RHS
573 *
574 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
575 IF( IERR.GT.0 )
576 $ INFO = IERR
577 IF( IJOB.EQ.0 ) THEN
578 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
579 $ SCALOC )
580 IF( SCALOC.NE.ONE ) THEN
581 DO 90 K = 1, N
582 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
583 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
584 90 CONTINUE
585 SCALE = SCALE*SCALOC
586 END IF
587 ELSE
588 CALL DLATDF( IJOB, ZDIM, Z, LDZ, RHS, RDSUM,
589 $ RDSCAL, IPIV, JPIV )
590 END IF
591 *
592 * Unpack solution vector(s)
593 *
594 K = 1
595 II = MB*NB + 1
596 DO 100 JJ = 0, NB - 1
597 CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
598 CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
599 K = K + MB
600 II = II + MB
601 100 CONTINUE
602 *
603 * Substitute R(I, J) and L(I, J) into remaining
604 * equation.
605 *
606 IF( I.GT.1 ) THEN
607 CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
608 $ A( 1, IS ), LDA, RHS( 1 ), MB, ONE,
609 $ C( 1, JS ), LDC )
610 CALL DGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
611 $ D( 1, IS ), LDD, RHS( 1 ), MB, ONE,
612 $ F( 1, JS ), LDF )
613 END IF
614 IF( J.LT.Q ) THEN
615 K = MB*NB + 1
616 CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
617 $ MB, B( JS, JE+1 ), LDB, ONE,
618 $ C( IS, JE+1 ), LDC )
619 CALL DGEMM( 'N', 'N', MB, N-JE, NB, ONE, RHS( K ),
620 $ MB, E( JS, JE+1 ), LDE, ONE,
621 $ F( IS, JE+1 ), LDF )
622 END IF
623 *
624 END IF
625 *
626 110 CONTINUE
627 120 CONTINUE
628 ELSE
629 *
630 * Solve (I, J) - subsystem
631 * A(I, I)**T * R(I, J) + D(I, I)**T * L(J, J) = C(I, J)
632 * R(I, I) * B(J, J) + L(I, J) * E(J, J) = -F(I, J)
633 * for I = 1, 2, ..., P, J = Q, Q - 1, ..., 1
634 *
635 SCALE = ONE
636 SCALOC = ONE
637 DO 200 I = 1, P
638 *
639 IS = IWORK( I )
640 ISP1 = IS + 1
641 IE = IWORK ( I+1 ) - 1
642 MB = IE - IS + 1
643 DO 190 J = Q, P + 2, -1
644 *
645 JS = IWORK( J )
646 JSP1 = JS + 1
647 JE = IWORK( J+1 ) - 1
648 NB = JE - JS + 1
649 ZDIM = MB*NB*2
650 IF( ( MB.EQ.1 ) .AND. ( NB.EQ.1 ) ) THEN
651 *
652 * Build a 2-by-2 system Z**T * x = RHS
653 *
654 Z( 1, 1 ) = A( IS, IS )
655 Z( 2, 1 ) = -B( JS, JS )
656 Z( 1, 2 ) = D( IS, IS )
657 Z( 2, 2 ) = -E( JS, JS )
658 *
659 * Set up right hand side(s)
660 *
661 RHS( 1 ) = C( IS, JS )
662 RHS( 2 ) = F( IS, JS )
663 *
664 * Solve Z**T * x = RHS
665 *
666 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
667 IF( IERR.GT.0 )
668 $ INFO = IERR
669 *
670 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
671 IF( SCALOC.NE.ONE ) THEN
672 DO 130 K = 1, N
673 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
674 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
675 130 CONTINUE
676 SCALE = SCALE*SCALOC
677 END IF
678 *
679 * Unpack solution vector(s)
680 *
681 C( IS, JS ) = RHS( 1 )
682 F( IS, JS ) = RHS( 2 )
683 *
684 * Substitute R(I, J) and L(I, J) into remaining
685 * equation.
686 *
687 IF( J.GT.P+2 ) THEN
688 ALPHA = RHS( 1 )
689 CALL DAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
690 $ LDF )
691 ALPHA = RHS( 2 )
692 CALL DAXPY( JS-1, ALPHA, E( 1, JS ), 1, F( IS, 1 ),
693 $ LDF )
694 END IF
695 IF( I.LT.P ) THEN
696 ALPHA = -RHS( 1 )
697 CALL DAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
698 $ C( IE+1, JS ), 1 )
699 ALPHA = -RHS( 2 )
700 CALL DAXPY( M-IE, ALPHA, D( IS, IE+1 ), LDD,
701 $ C( IE+1, JS ), 1 )
702 END IF
703 *
704 ELSE IF( ( MB.EQ.1 ) .AND. ( NB.EQ.2 ) ) THEN
705 *
706 * Build a 4-by-4 system Z**T * x = RHS
707 *
708 Z( 1, 1 ) = A( IS, IS )
709 Z( 2, 1 ) = ZERO
710 Z( 3, 1 ) = -B( JS, JS )
711 Z( 4, 1 ) = -B( JSP1, JS )
712 *
713 Z( 1, 2 ) = ZERO
714 Z( 2, 2 ) = A( IS, IS )
715 Z( 3, 2 ) = -B( JS, JSP1 )
716 Z( 4, 2 ) = -B( JSP1, JSP1 )
717 *
718 Z( 1, 3 ) = D( IS, IS )
719 Z( 2, 3 ) = ZERO
720 Z( 3, 3 ) = -E( JS, JS )
721 Z( 4, 3 ) = ZERO
722 *
723 Z( 1, 4 ) = ZERO
724 Z( 2, 4 ) = D( IS, IS )
725 Z( 3, 4 ) = -E( JS, JSP1 )
726 Z( 4, 4 ) = -E( JSP1, JSP1 )
727 *
728 * Set up right hand side(s)
729 *
730 RHS( 1 ) = C( IS, JS )
731 RHS( 2 ) = C( IS, JSP1 )
732 RHS( 3 ) = F( IS, JS )
733 RHS( 4 ) = F( IS, JSP1 )
734 *
735 * Solve Z**T * x = RHS
736 *
737 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
738 IF( IERR.GT.0 )
739 $ INFO = IERR
740 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
741 IF( SCALOC.NE.ONE ) THEN
742 DO 140 K = 1, N
743 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
744 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
745 140 CONTINUE
746 SCALE = SCALE*SCALOC
747 END IF
748 *
749 * Unpack solution vector(s)
750 *
751 C( IS, JS ) = RHS( 1 )
752 C( IS, JSP1 ) = RHS( 2 )
753 F( IS, JS ) = RHS( 3 )
754 F( IS, JSP1 ) = RHS( 4 )
755 *
756 * Substitute R(I, J) and L(I, J) into remaining
757 * equation.
758 *
759 IF( J.GT.P+2 ) THEN
760 CALL DAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
761 $ F( IS, 1 ), LDF )
762 CALL DAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
763 $ F( IS, 1 ), LDF )
764 CALL DAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
765 $ F( IS, 1 ), LDF )
766 CALL DAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
767 $ F( IS, 1 ), LDF )
768 END IF
769 IF( I.LT.P ) THEN
770 CALL DGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,
771 $ RHS( 1 ), 1, C( IE+1, JS ), LDC )
772 CALL DGER( M-IE, NB, -ONE, D( IS, IE+1 ), LDD,
773 $ RHS( 3 ), 1, C( IE+1, JS ), LDC )
774 END IF
775 *
776 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.1 ) ) THEN
777 *
778 * Build a 4-by-4 system Z**T * x = RHS
779 *
780 Z( 1, 1 ) = A( IS, IS )
781 Z( 2, 1 ) = A( IS, ISP1 )
782 Z( 3, 1 ) = -B( JS, JS )
783 Z( 4, 1 ) = ZERO
784 *
785 Z( 1, 2 ) = A( ISP1, IS )
786 Z( 2, 2 ) = A( ISP1, ISP1 )
787 Z( 3, 2 ) = ZERO
788 Z( 4, 2 ) = -B( JS, JS )
789 *
790 Z( 1, 3 ) = D( IS, IS )
791 Z( 2, 3 ) = D( IS, ISP1 )
792 Z( 3, 3 ) = -E( JS, JS )
793 Z( 4, 3 ) = ZERO
794 *
795 Z( 1, 4 ) = ZERO
796 Z( 2, 4 ) = D( ISP1, ISP1 )
797 Z( 3, 4 ) = ZERO
798 Z( 4, 4 ) = -E( JS, JS )
799 *
800 * Set up right hand side(s)
801 *
802 RHS( 1 ) = C( IS, JS )
803 RHS( 2 ) = C( ISP1, JS )
804 RHS( 3 ) = F( IS, JS )
805 RHS( 4 ) = F( ISP1, JS )
806 *
807 * Solve Z**T * x = RHS
808 *
809 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
810 IF( IERR.GT.0 )
811 $ INFO = IERR
812 *
813 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
814 IF( SCALOC.NE.ONE ) THEN
815 DO 150 K = 1, N
816 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
817 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
818 150 CONTINUE
819 SCALE = SCALE*SCALOC
820 END IF
821 *
822 * Unpack solution vector(s)
823 *
824 C( IS, JS ) = RHS( 1 )
825 C( ISP1, JS ) = RHS( 2 )
826 F( IS, JS ) = RHS( 3 )
827 F( ISP1, JS ) = RHS( 4 )
828 *
829 * Substitute R(I, J) and L(I, J) into remaining
830 * equation.
831 *
832 IF( J.GT.P+2 ) THEN
833 CALL DGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
834 $ 1, F( IS, 1 ), LDF )
835 CALL DGER( MB, JS-1, ONE, RHS( 3 ), 1, E( 1, JS ),
836 $ 1, F( IS, 1 ), LDF )
837 END IF
838 IF( I.LT.P ) THEN
839 CALL DGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ),
840 $ LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
841 $ 1 )
842 CALL DGEMV( 'T', MB, M-IE, -ONE, D( IS, IE+1 ),
843 $ LDD, RHS( 3 ), 1, ONE, C( IE+1, JS ),
844 $ 1 )
845 END IF
846 *
847 ELSE IF( ( MB.EQ.2 ) .AND. ( NB.EQ.2 ) ) THEN
848 *
849 * Build an 8-by-8 system Z**T * x = RHS
850 *
851 CALL DLASET( 'F', LDZ, LDZ, ZERO, ZERO, Z, LDZ )
852 *
853 Z( 1, 1 ) = A( IS, IS )
854 Z( 2, 1 ) = A( IS, ISP1 )
855 Z( 5, 1 ) = -B( JS, JS )
856 Z( 7, 1 ) = -B( JSP1, JS )
857 *
858 Z( 1, 2 ) = A( ISP1, IS )
859 Z( 2, 2 ) = A( ISP1, ISP1 )
860 Z( 6, 2 ) = -B( JS, JS )
861 Z( 8, 2 ) = -B( JSP1, JS )
862 *
863 Z( 3, 3 ) = A( IS, IS )
864 Z( 4, 3 ) = A( IS, ISP1 )
865 Z( 5, 3 ) = -B( JS, JSP1 )
866 Z( 7, 3 ) = -B( JSP1, JSP1 )
867 *
868 Z( 3, 4 ) = A( ISP1, IS )
869 Z( 4, 4 ) = A( ISP1, ISP1 )
870 Z( 6, 4 ) = -B( JS, JSP1 )
871 Z( 8, 4 ) = -B( JSP1, JSP1 )
872 *
873 Z( 1, 5 ) = D( IS, IS )
874 Z( 2, 5 ) = D( IS, ISP1 )
875 Z( 5, 5 ) = -E( JS, JS )
876 *
877 Z( 2, 6 ) = D( ISP1, ISP1 )
878 Z( 6, 6 ) = -E( JS, JS )
879 *
880 Z( 3, 7 ) = D( IS, IS )
881 Z( 4, 7 ) = D( IS, ISP1 )
882 Z( 5, 7 ) = -E( JS, JSP1 )
883 Z( 7, 7 ) = -E( JSP1, JSP1 )
884 *
885 Z( 4, 8 ) = D( ISP1, ISP1 )
886 Z( 6, 8 ) = -E( JS, JSP1 )
887 Z( 8, 8 ) = -E( JSP1, JSP1 )
888 *
889 * Set up right hand side(s)
890 *
891 K = 1
892 II = MB*NB + 1
893 DO 160 JJ = 0, NB - 1
894 CALL DCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
895 CALL DCOPY( MB, F( IS, JS+JJ ), 1, RHS( II ), 1 )
896 K = K + MB
897 II = II + MB
898 160 CONTINUE
899 *
900 *
901 * Solve Z**T * x = RHS
902 *
903 CALL DGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
904 IF( IERR.GT.0 )
905 $ INFO = IERR
906 *
907 CALL DGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
908 IF( SCALOC.NE.ONE ) THEN
909 DO 170 K = 1, N
910 CALL DSCAL( M, SCALOC, C( 1, K ), 1 )
911 CALL DSCAL( M, SCALOC, F( 1, K ), 1 )
912 170 CONTINUE
913 SCALE = SCALE*SCALOC
914 END IF
915 *
916 * Unpack solution vector(s)
917 *
918 K = 1
919 II = MB*NB + 1
920 DO 180 JJ = 0, NB - 1
921 CALL DCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
922 CALL DCOPY( MB, RHS( II ), 1, F( IS, JS+JJ ), 1 )
923 K = K + MB
924 II = II + MB
925 180 CONTINUE
926 *
927 * Substitute R(I, J) and L(I, J) into remaining
928 * equation.
929 *
930 IF( J.GT.P+2 ) THEN
931 CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
932 $ C( IS, JS ), LDC, B( 1, JS ), LDB, ONE,
933 $ F( IS, 1 ), LDF )
934 CALL DGEMM( 'N', 'T', MB, JS-1, NB, ONE,
935 $ F( IS, JS ), LDF, E( 1, JS ), LDE, ONE,
936 $ F( IS, 1 ), LDF )
937 END IF
938 IF( I.LT.P ) THEN
939 CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
940 $ A( IS, IE+1 ), LDA, C( IS, JS ), LDC,
941 $ ONE, C( IE+1, JS ), LDC )
942 CALL DGEMM( 'T', 'N', M-IE, NB, MB, -ONE,
943 $ D( IS, IE+1 ), LDD, F( IS, JS ), LDF,
944 $ ONE, C( IE+1, JS ), LDC )
945 END IF
946 *
947 END IF
948 *
949 190 CONTINUE
950 200 CONTINUE
951 *
952 END IF
953 RETURN
954 *
955 * End of DTGSY2
956 *
957 END