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