001:       SUBROUTINE STGSY2( TRANS, IJOB, M, N, A, LDA, B, LDB, C, LDC, D,
002:      $                   LDD, E, LDE, F, LDF, SCALE, RDSUM, RDSCAL,
003:      $                   IWORK, PQ, 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: *     January 2007
009: *
010: *     .. Scalar Arguments ..
011:       CHARACTER          TRANS
012:       INTEGER            IJOB, INFO, LDA, LDB, LDC, LDD, LDE, LDF, M, N,
013:      $                   PQ
014:       REAL               RDSCAL, RDSUM, SCALE
015: *     ..
016: *     .. Array Arguments ..
017:       INTEGER            IWORK( * )
018:       REAL               A( LDA, * ), B( LDB, * ), C( LDC, * ),
019:      $                   D( LDD, * ), E( LDE, * ), F( LDF, * )
020: *     ..
021: *
022: *  Purpose
023: *  =======
024: *
025: *  STGSY2 solves the generalized Sylvester equation:
026: *
027: *              A * R - L * B = scale * C                (1)
028: *              D * R - L * E = scale * F,
029: *
030: *  using Level 1 and 2 BLAS. where R and L are unknown M-by-N matrices,
031: *  (A, D), (B, E) and (C, F) are given matrix pairs of size M-by-M,
032: *  N-by-N and M-by-N, respectively, with real entries. (A, D) and (B, E)
033: *  must be in generalized Schur canonical form, i.e. A, B are upper
034: *  quasi triangular and D, E are upper triangular. The solution (R, L)
035: *  overwrites (C, F). 0 <= SCALE <= 1 is an output scaling factor
036: *  chosen to avoid overflow.
037: *
038: *  In matrix notation solving equation (1) corresponds to solve
039: *  Z*x = scale*b, where Z is defined as
040: *
041: *         Z = [ kron(In, A)  -kron(B', Im) ]             (2)
042: *             [ kron(In, D)  -kron(E', Im) ],
043: *
044: *  Ik is the identity matrix of size k and X' is the transpose of X.
045: *  kron(X, Y) is the Kronecker product between the matrices X and Y.
046: *  In the process of solving (1), we solve a number of such systems
047: *  where Dim(In), Dim(In) = 1 or 2.
048: *
049: *  If TRANS = 'T', solve the transposed system Z'*y = scale*b for y,
050: *  which is equivalent to solve for R and L in
051: *
052: *              A' * R  + D' * L   = scale *  C           (3)
053: *              R  * B' + L  * E'  = scale * -F
054: *
055: *  This case is used to compute an estimate of Dif[(A, D), (B, E)] =
056: *  sigma_min(Z) using reverse communicaton with SLACON.
057: *
058: *  STGSY2 also (IJOB >= 1) contributes to the computation in STGSYL
059: *  of an upper bound on the separation between to matrix pairs. Then
060: *  the input (A, D), (B, E) are sub-pencils of the matrix pair in
061: *  STGSYL. See STGSYL for details.
062: *
063: *  Arguments
064: *  =========
065: *
066: *  TRANS   (input) CHARACTER*1
067: *          = 'N', solve the generalized Sylvester equation (1).
068: *          = 'T': solve the 'transposed' system (3).
069: *
070: *  IJOB    (input) INTEGER
071: *          Specifies what kind of functionality to be performed.
072: *          = 0: solve (1) only.
073: *          = 1: A contribution from this subsystem to a Frobenius
074: *               norm-based estimate of the separation between two matrix
075: *               pairs is computed. (look ahead strategy is used).
076: *          = 2: A contribution from this subsystem to a Frobenius
077: *               norm-based estimate of the separation between two matrix
078: *               pairs is computed. (SGECON on sub-systems is used.)
079: *          Not referenced if TRANS = 'T'.
080: *
081: *  M       (input) INTEGER
082: *          On entry, M specifies the order of A and D, and the row
083: *          dimension of C, F, R and L.
084: *
085: *  N       (input) INTEGER
086: *          On entry, N specifies the order of B and E, and the column
087: *          dimension of C, F, R and L.
088: *
089: *  A       (input) REAL array, dimension (LDA, M)
090: *          On entry, A contains an upper quasi triangular matrix.
091: *
092: *  LDA     (input) INTEGER
093: *          The leading dimension of the matrix A. LDA >= max(1, M).
094: *
095: *  B       (input) REAL array, dimension (LDB, N)
096: *          On entry, B contains an upper quasi triangular matrix.
097: *
098: *  LDB     (input) INTEGER
099: *          The leading dimension of the matrix B. LDB >= max(1, N).
100: *
101: *  C       (input/output) REAL 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) REAL 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) REAL 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) REAL 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) REAL
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) REAL
140: *          On entry, the sum of squares of computed contributions to
141: *          the Dif-estimate under computation by STGSYL, 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 STGSY2 is called by STGSYL.
147: *
148: *  RDSCAL  (input/output) REAL
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 STGSY2 is called by
154: *                STGSYL.
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 SCOPY by calls to SLASET.
178: *  Sven Hammarling, 27/5/02.
179: *
180: *     .. Parameters ..
181:       INTEGER            LDZ
182:       PARAMETER          ( LDZ = 8 )
183:       REAL               ZERO, ONE
184:       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+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:       REAL               ALPHA, SCALOC
191: *     ..
192: *     .. Local Arrays ..
193:       INTEGER            IPIV( LDZ ), JPIV( LDZ )
194:       REAL               RHS( LDZ ), Z( LDZ, LDZ )
195: *     ..
196: *     .. External Functions ..
197:       LOGICAL            LSAME
198:       EXTERNAL           LSAME
199: *     ..
200: *     .. External Subroutines ..
201:       EXTERNAL           SAXPY, SCOPY, SGEMM, SGEMV, SGER, SGESC2,
202:      $                   SGETC2, SSCAL, SLASET, SLATDF, 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( 'STGSY2', -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 SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
326:                   IF( IERR.GT.0 )
327:      $               INFO = IERR
328: *
329:                   IF( IJOB.EQ.0 ) THEN
330:                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
331:      $                            SCALOC )
332:                      IF( SCALOC.NE.ONE ) THEN
333:                         DO 50 K = 1, N
334:                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
335:                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
336:    50                   CONTINUE
337:                         SCALE = SCALE*SCALOC
338:                      END IF
339:                   ELSE
340:                      CALL SLATDF( 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 SAXPY( IS-1, ALPHA, A( 1, IS ), 1, C( 1, JS ),
355:      $                           1 )
356:                      CALL SAXPY( IS-1, ALPHA, D( 1, IS ), 1, F( 1, JS ),
357:      $                           1 )
358:                   END IF
359:                   IF( J.LT.Q ) THEN
360:                      CALL SAXPY( N-JE, RHS( 2 ), B( JS, JE+1 ), LDB,
361:      $                           C( IS, JE+1 ), LDC )
362:                      CALL SAXPY( 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 SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
400:                   IF( IERR.GT.0 )
401:      $               INFO = IERR
402: *
403:                   IF( IJOB.EQ.0 ) THEN
404:                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
405:      $                            SCALOC )
406:                      IF( SCALOC.NE.ONE ) THEN
407:                         DO 60 K = 1, N
408:                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
409:                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
410:    60                   CONTINUE
411:                         SCALE = SCALE*SCALOC
412:                      END IF
413:                   ELSE
414:                      CALL SLATDF( 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 SGER( IS-1, NB, -ONE, A( 1, IS ), 1, RHS( 1 ),
430:      $                          1, C( 1, JS ), LDC )
431:                      CALL SGER( 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 SAXPY( N-JE, RHS( 3 ), B( JS, JE+1 ), LDB,
436:      $                           C( IS, JE+1 ), LDC )
437:                      CALL SAXPY( N-JE, RHS( 3 ), E( JS, JE+1 ), LDE,
438:      $                           F( IS, JE+1 ), LDF )
439:                      CALL SAXPY( N-JE, RHS( 4 ), B( JSP1, JE+1 ), LDB,
440:      $                           C( IS, JE+1 ), LDC )
441:                      CALL SAXPY( 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 SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
479:                   IF( IERR.GT.0 )
480:      $               INFO = IERR
481:                   IF( IJOB.EQ.0 ) THEN
482:                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
483:      $                            SCALOC )
484:                      IF( SCALOC.NE.ONE ) THEN
485:                         DO 70 K = 1, N
486:                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
487:                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
488:    70                   CONTINUE
489:                         SCALE = SCALE*SCALOC
490:                      END IF
491:                   ELSE
492:                      CALL SLATDF( 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 SGEMV( 'N', IS-1, MB, -ONE, A( 1, IS ), LDA,
508:      $                           RHS( 1 ), 1, ONE, C( 1, JS ), 1 )
509:                      CALL SGEMV( '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 SGER( MB, N-JE, ONE, RHS( 3 ), 1,
514:      $                          B( JS, JE+1 ), LDB, C( IS, JE+1 ), LDC )
515:                      CALL SGER( 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 SLASET( '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 SCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
567:                      CALL SCOPY( 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 SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
575:                   IF( IERR.GT.0 )
576:      $               INFO = IERR
577:                   IF( IJOB.EQ.0 ) THEN
578:                      CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV,
579:      $                            SCALOC )
580:                      IF( SCALOC.NE.ONE ) THEN
581:                         DO 90 K = 1, N
582:                            CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
583:                            CALL SSCAL( M, SCALOC, F( 1, K ), 1 )
584:    90                   CONTINUE
585:                         SCALE = SCALE*SCALOC
586:                      END IF
587:                   ELSE
588:                      CALL SLATDF( 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 SCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
598:                      CALL SCOPY( 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 SGEMM( 'N', 'N', IS-1, NB, MB, -ONE,
608:      $                           A( 1, IS ), LDA, RHS( 1 ), MB, ONE,
609:      $                           C( 1, JS ), LDC )
610:                      CALL SGEMM( '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 SGEMM( '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 SGEMM( '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)' * R(I, J) + D(I, I)' * 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' * 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' * x = RHS
665: *
666:                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
667:                   IF( IERR.GT.0 )
668:      $               INFO = IERR
669: *
670:                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
671:                   IF( SCALOC.NE.ONE ) THEN
672:                      DO 130 K = 1, N
673:                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
674:                         CALL SSCAL( 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 SAXPY( JS-1, ALPHA, B( 1, JS ), 1, F( IS, 1 ),
690:      $                           LDF )
691:                      ALPHA = RHS( 2 )
692:                      CALL SAXPY( 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 SAXPY( M-IE, ALPHA, A( IS, IE+1 ), LDA,
698:      $                           C( IE+1, JS ), 1 )
699:                      ALPHA = -RHS( 2 )
700:                      CALL SAXPY( 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' * 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' * x = RHS
736: *
737:                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
738:                   IF( IERR.GT.0 )
739:      $               INFO = IERR
740:                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
741:                   IF( SCALOC.NE.ONE ) THEN
742:                      DO 140 K = 1, N
743:                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
744:                         CALL SSCAL( 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 SAXPY( JS-1, RHS( 1 ), B( 1, JS ), 1,
761:      $                           F( IS, 1 ), LDF )
762:                      CALL SAXPY( JS-1, RHS( 2 ), B( 1, JSP1 ), 1,
763:      $                           F( IS, 1 ), LDF )
764:                      CALL SAXPY( JS-1, RHS( 3 ), E( 1, JS ), 1,
765:      $                           F( IS, 1 ), LDF )
766:                      CALL SAXPY( JS-1, RHS( 4 ), E( 1, JSP1 ), 1,
767:      $                           F( IS, 1 ), LDF )
768:                   END IF
769:                   IF( I.LT.P ) THEN
770:                      CALL SGER( M-IE, NB, -ONE, A( IS, IE+1 ), LDA,
771:      $                          RHS( 1 ), 1, C( IE+1, JS ), LDC )
772:                      CALL SGER( 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' * 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' * x = RHS
808: *
809:                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
810:                   IF( IERR.GT.0 )
811:      $               INFO = IERR
812: *
813:                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
814:                   IF( SCALOC.NE.ONE ) THEN
815:                      DO 150 K = 1, N
816:                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
817:                         CALL SSCAL( 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 SGER( MB, JS-1, ONE, RHS( 1 ), 1, B( 1, JS ),
834:      $                          1, F( IS, 1 ), LDF )
835:                      CALL SGER( 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 SGEMV( 'T', MB, M-IE, -ONE, A( IS, IE+1 ),
840:      $                           LDA, RHS( 1 ), 1, ONE, C( IE+1, JS ),
841:      $                           1 )
842:                      CALL SGEMV( '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' * x = RHS
850: *
851:                   CALL SLASET( '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 SCOPY( MB, C( IS, JS+JJ ), 1, RHS( K ), 1 )
895:                      CALL SCOPY( 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' * x = RHS
902: *
903:                   CALL SGETC2( ZDIM, Z, LDZ, IPIV, JPIV, IERR )
904:                   IF( IERR.GT.0 )
905:      $               INFO = IERR
906: *
907:                   CALL SGESC2( ZDIM, Z, LDZ, RHS, IPIV, JPIV, SCALOC )
908:                   IF( SCALOC.NE.ONE ) THEN
909:                      DO 170 K = 1, N
910:                         CALL SSCAL( M, SCALOC, C( 1, K ), 1 )
911:                         CALL SSCAL( 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 SCOPY( MB, RHS( K ), 1, C( IS, JS+JJ ), 1 )
922:                      CALL SCOPY( 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 SGEMM( 'N', 'T', MB, JS-1, NB, ONE,
932:      $                           C( IS, JS ), LDC, B( 1, JS ), LDB, ONE,
933:      $                           F( IS, 1 ), LDF )
934:                      CALL SGEMM( '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 SGEMM( '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 SGEMM( '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 STGSY2
956: *
957:       END
958: