001:       SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
002:      $                   LDC, SCALE, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          TRANA, TRANB
011:       INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N
012:       DOUBLE PRECISION   SCALE
013: *     ..
014: *     .. Array Arguments ..
015:       DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * )
016: *     ..
017: *
018: *  Purpose
019: *  =======
020: *
021: *  DTRSYL solves the real Sylvester matrix equation:
022: *
023: *     op(A)*X + X*op(B) = scale*C or
024: *     op(A)*X - X*op(B) = scale*C,
025: *
026: *  where op(A) = A or A**T, and  A and B are both upper quasi-
027: *  triangular. A is M-by-M and B is N-by-N; the right hand side C and
028: *  the solution X are M-by-N; and scale is an output scale factor, set
029: *  <= 1 to avoid overflow in X.
030: *
031: *  A and B must be in Schur canonical form (as returned by DHSEQR), that
032: *  is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
033: *  each 2-by-2 diagonal block has its diagonal elements equal and its
034: *  off-diagonal elements of opposite sign.
035: *
036: *  Arguments
037: *  =========
038: *
039: *  TRANA   (input) CHARACTER*1
040: *          Specifies the option op(A):
041: *          = 'N': op(A) = A    (No transpose)
042: *          = 'T': op(A) = A**T (Transpose)
043: *          = 'C': op(A) = A**H (Conjugate transpose = Transpose)
044: *
045: *  TRANB   (input) CHARACTER*1
046: *          Specifies the option op(B):
047: *          = 'N': op(B) = B    (No transpose)
048: *          = 'T': op(B) = B**T (Transpose)
049: *          = 'C': op(B) = B**H (Conjugate transpose = Transpose)
050: *
051: *  ISGN    (input) INTEGER
052: *          Specifies the sign in the equation:
053: *          = +1: solve op(A)*X + X*op(B) = scale*C
054: *          = -1: solve op(A)*X - X*op(B) = scale*C
055: *
056: *  M       (input) INTEGER
057: *          The order of the matrix A, and the number of rows in the
058: *          matrices X and C. M >= 0.
059: *
060: *  N       (input) INTEGER
061: *          The order of the matrix B, and the number of columns in the
062: *          matrices X and C. N >= 0.
063: *
064: *  A       (input) DOUBLE PRECISION array, dimension (LDA,M)
065: *          The upper quasi-triangular matrix A, in Schur canonical form.
066: *
067: *  LDA     (input) INTEGER
068: *          The leading dimension of the array A. LDA >= max(1,M).
069: *
070: *  B       (input) DOUBLE PRECISION array, dimension (LDB,N)
071: *          The upper quasi-triangular matrix B, in Schur canonical form.
072: *
073: *  LDB     (input) INTEGER
074: *          The leading dimension of the array B. LDB >= max(1,N).
075: *
076: *  C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)
077: *          On entry, the M-by-N right hand side matrix C.
078: *          On exit, C is overwritten by the solution matrix X.
079: *
080: *  LDC     (input) INTEGER
081: *          The leading dimension of the array C. LDC >= max(1,M)
082: *
083: *  SCALE   (output) DOUBLE PRECISION
084: *          The scale factor, scale, set <= 1 to avoid overflow in X.
085: *
086: *  INFO    (output) INTEGER
087: *          = 0: successful exit
088: *          < 0: if INFO = -i, the i-th argument had an illegal value
089: *          = 1: A and B have common or very close eigenvalues; perturbed
090: *               values were used to solve the equation (but the matrices
091: *               A and B are unchanged).
092: *
093: *  =====================================================================
094: *
095: *     .. Parameters ..
096:       DOUBLE PRECISION   ZERO, ONE
097:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
098: *     ..
099: *     .. Local Scalars ..
100:       LOGICAL            NOTRNA, NOTRNB
101:       INTEGER            IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
102:       DOUBLE PRECISION   A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
103:      $                   SMLNUM, SUML, SUMR, XNORM
104: *     ..
105: *     .. Local Arrays ..
106:       DOUBLE PRECISION   DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
107: *     ..
108: *     .. External Functions ..
109:       LOGICAL            LSAME
110:       DOUBLE PRECISION   DDOT, DLAMCH, DLANGE
111:       EXTERNAL           LSAME, DDOT, DLAMCH, DLANGE
112: *     ..
113: *     .. External Subroutines ..
114:       EXTERNAL           DLABAD, DLALN2, DLASY2, DSCAL, XERBLA
115: *     ..
116: *     .. Intrinsic Functions ..
117:       INTRINSIC          ABS, DBLE, MAX, MIN
118: *     ..
119: *     .. Executable Statements ..
120: *
121: *     Decode and Test input parameters
122: *
123:       NOTRNA = LSAME( TRANA, 'N' )
124:       NOTRNB = LSAME( TRANB, 'N' )
125: *
126:       INFO = 0
127:       IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT.
128:      $    LSAME( TRANA, 'C' ) ) THEN
129:          INFO = -1
130:       ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT.
131:      $         LSAME( TRANB, 'C' ) ) THEN
132:          INFO = -2
133:       ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
134:          INFO = -3
135:       ELSE IF( M.LT.0 ) THEN
136:          INFO = -4
137:       ELSE IF( N.LT.0 ) THEN
138:          INFO = -5
139:       ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
140:          INFO = -7
141:       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
142:          INFO = -9
143:       ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
144:          INFO = -11
145:       END IF
146:       IF( INFO.NE.0 ) THEN
147:          CALL XERBLA( 'DTRSYL', -INFO )
148:          RETURN
149:       END IF
150: *
151: *     Quick return if possible
152: *
153:       SCALE = ONE
154:       IF( M.EQ.0 .OR. N.EQ.0 )
155:      $   RETURN
156: *
157: *     Set constants to control overflow
158: *
159:       EPS = DLAMCH( 'P' )
160:       SMLNUM = DLAMCH( 'S' )
161:       BIGNUM = ONE / SMLNUM
162:       CALL DLABAD( SMLNUM, BIGNUM )
163:       SMLNUM = SMLNUM*DBLE( M*N ) / EPS
164:       BIGNUM = ONE / SMLNUM
165: *
166:       SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ),
167:      $       EPS*DLANGE( 'M', N, N, B, LDB, DUM ) )
168: *
169:       SGN = ISGN
170: *
171:       IF( NOTRNA .AND. NOTRNB ) THEN
172: *
173: *        Solve    A*X + ISGN*X*B = scale*C.
174: *
175: *        The (K,L)th block of X is determined starting from
176: *        bottom-left corner column by column by
177: *
178: *         A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
179: *
180: *        Where
181: *                  M                         L-1
182: *        R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
183: *                I=K+1                       J=1
184: *
185: *        Start column loop (index = L)
186: *        L1 (L2) : column index of the first (first) row of X(K,L).
187: *
188:          LNEXT = 1
189:          DO 60 L = 1, N
190:             IF( L.LT.LNEXT )
191:      $         GO TO 60
192:             IF( L.EQ.N ) THEN
193:                L1 = L
194:                L2 = L
195:             ELSE
196:                IF( B( L+1, L ).NE.ZERO ) THEN
197:                   L1 = L
198:                   L2 = L + 1
199:                   LNEXT = L + 2
200:                ELSE
201:                   L1 = L
202:                   L2 = L
203:                   LNEXT = L + 1
204:                END IF
205:             END IF
206: *
207: *           Start row loop (index = K)
208: *           K1 (K2): row index of the first (last) row of X(K,L).
209: *
210:             KNEXT = M
211:             DO 50 K = M, 1, -1
212:                IF( K.GT.KNEXT )
213:      $            GO TO 50
214:                IF( K.EQ.1 ) THEN
215:                   K1 = K
216:                   K2 = K
217:                ELSE
218:                   IF( A( K, K-1 ).NE.ZERO ) THEN
219:                      K1 = K - 1
220:                      K2 = K
221:                      KNEXT = K - 2
222:                   ELSE
223:                      K1 = K
224:                      K2 = K
225:                      KNEXT = K - 1
226:                   END IF
227:                END IF
228: *
229:                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
230:                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
231:      $                   C( MIN( K1+1, M ), L1 ), 1 )
232:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
233:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
234:                   SCALOC = ONE
235: *
236:                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
237:                   DA11 = ABS( A11 )
238:                   IF( DA11.LE.SMIN ) THEN
239:                      A11 = SMIN
240:                      DA11 = SMIN
241:                      INFO = 1
242:                   END IF
243:                   DB = ABS( VEC( 1, 1 ) )
244:                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
245:                      IF( DB.GT.BIGNUM*DA11 )
246:      $                  SCALOC = ONE / DB
247:                   END IF
248:                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
249: *
250:                   IF( SCALOC.NE.ONE ) THEN
251:                      DO 10 J = 1, N
252:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
253:    10                CONTINUE
254:                      SCALE = SCALE*SCALOC
255:                   END IF
256:                   C( K1, L1 ) = X( 1, 1 )
257: *
258:                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
259: *
260:                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
261:      $                   C( MIN( K2+1, M ), L1 ), 1 )
262:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
263:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
264: *
265:                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
266:      $                   C( MIN( K2+1, M ), L1 ), 1 )
267:                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
268:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
269: *
270:                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
271:      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
272:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
273:                   IF( IERR.NE.0 )
274:      $               INFO = 1
275: *
276:                   IF( SCALOC.NE.ONE ) THEN
277:                      DO 20 J = 1, N
278:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
279:    20                CONTINUE
280:                      SCALE = SCALE*SCALOC
281:                   END IF
282:                   C( K1, L1 ) = X( 1, 1 )
283:                   C( K2, L1 ) = X( 2, 1 )
284: *
285:                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
286: *
287:                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
288:      $                   C( MIN( K1+1, M ), L1 ), 1 )
289:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
290:                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
291: *
292:                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
293:      $                   C( MIN( K1+1, M ), L2 ), 1 )
294:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
295:                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
296: *
297:                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
298:      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
299:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
300:                   IF( IERR.NE.0 )
301:      $               INFO = 1
302: *
303:                   IF( SCALOC.NE.ONE ) THEN
304:                      DO 30 J = 1, N
305:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
306:    30                CONTINUE
307:                      SCALE = SCALE*SCALOC
308:                   END IF
309:                   C( K1, L1 ) = X( 1, 1 )
310:                   C( K1, L2 ) = X( 2, 1 )
311: *
312:                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
313: *
314:                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
315:      $                   C( MIN( K2+1, M ), L1 ), 1 )
316:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
317:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
318: *
319:                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
320:      $                   C( MIN( K2+1, M ), L2 ), 1 )
321:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
322:                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
323: *
324:                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
325:      $                   C( MIN( K2+1, M ), L1 ), 1 )
326:                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
327:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
328: *
329:                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
330:      $                   C( MIN( K2+1, M ), L2 ), 1 )
331:                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
332:                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
333: *
334:                   CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2,
335:      $                         A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC,
336:      $                         2, SCALOC, X, 2, XNORM, IERR )
337:                   IF( IERR.NE.0 )
338:      $               INFO = 1
339: *
340:                   IF( SCALOC.NE.ONE ) THEN
341:                      DO 40 J = 1, N
342:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
343:    40                CONTINUE
344:                      SCALE = SCALE*SCALOC
345:                   END IF
346:                   C( K1, L1 ) = X( 1, 1 )
347:                   C( K1, L2 ) = X( 1, 2 )
348:                   C( K2, L1 ) = X( 2, 1 )
349:                   C( K2, L2 ) = X( 2, 2 )
350:                END IF
351: *
352:    50       CONTINUE
353: *
354:    60    CONTINUE
355: *
356:       ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
357: *
358: *        Solve    A' *X + ISGN*X*B = scale*C.
359: *
360: *        The (K,L)th block of X is determined starting from
361: *        upper-left corner column by column by
362: *
363: *          A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
364: *
365: *        Where
366: *                   K-1                        L-1
367: *          R(K,L) = SUM [A(I,K)'*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
368: *                   I=1                        J=1
369: *
370: *        Start column loop (index = L)
371: *        L1 (L2): column index of the first (last) row of X(K,L)
372: *
373:          LNEXT = 1
374:          DO 120 L = 1, N
375:             IF( L.LT.LNEXT )
376:      $         GO TO 120
377:             IF( L.EQ.N ) THEN
378:                L1 = L
379:                L2 = L
380:             ELSE
381:                IF( B( L+1, L ).NE.ZERO ) THEN
382:                   L1 = L
383:                   L2 = L + 1
384:                   LNEXT = L + 2
385:                ELSE
386:                   L1 = L
387:                   L2 = L
388:                   LNEXT = L + 1
389:                END IF
390:             END IF
391: *
392: *           Start row loop (index = K)
393: *           K1 (K2): row index of the first (last) row of X(K,L)
394: *
395:             KNEXT = 1
396:             DO 110 K = 1, M
397:                IF( K.LT.KNEXT )
398:      $            GO TO 110
399:                IF( K.EQ.M ) THEN
400:                   K1 = K
401:                   K2 = K
402:                ELSE
403:                   IF( A( K+1, K ).NE.ZERO ) THEN
404:                      K1 = K
405:                      K2 = K + 1
406:                      KNEXT = K + 2
407:                   ELSE
408:                      K1 = K
409:                      K2 = K
410:                      KNEXT = K + 1
411:                   END IF
412:                END IF
413: *
414:                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
415:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
416:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
417:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
418:                   SCALOC = ONE
419: *
420:                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
421:                   DA11 = ABS( A11 )
422:                   IF( DA11.LE.SMIN ) THEN
423:                      A11 = SMIN
424:                      DA11 = SMIN
425:                      INFO = 1
426:                   END IF
427:                   DB = ABS( VEC( 1, 1 ) )
428:                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
429:                      IF( DB.GT.BIGNUM*DA11 )
430:      $                  SCALOC = ONE / DB
431:                   END IF
432:                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
433: *
434:                   IF( SCALOC.NE.ONE ) THEN
435:                      DO 70 J = 1, N
436:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
437:    70                CONTINUE
438:                      SCALE = SCALE*SCALOC
439:                   END IF
440:                   C( K1, L1 ) = X( 1, 1 )
441: *
442:                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
443: *
444:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
445:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
446:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
447: *
448:                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
449:                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
450:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
451: *
452:                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
453:      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
454:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
455:                   IF( IERR.NE.0 )
456:      $               INFO = 1
457: *
458:                   IF( SCALOC.NE.ONE ) THEN
459:                      DO 80 J = 1, N
460:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
461:    80                CONTINUE
462:                      SCALE = SCALE*SCALOC
463:                   END IF
464:                   C( K1, L1 ) = X( 1, 1 )
465:                   C( K2, L1 ) = X( 2, 1 )
466: *
467:                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
468: *
469:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
470:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
471:                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
472: *
473:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
474:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
475:                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
476: *
477:                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
478:      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
479:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
480:                   IF( IERR.NE.0 )
481:      $               INFO = 1
482: *
483:                   IF( SCALOC.NE.ONE ) THEN
484:                      DO 90 J = 1, N
485:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
486:    90                CONTINUE
487:                      SCALE = SCALE*SCALOC
488:                   END IF
489:                   C( K1, L1 ) = X( 1, 1 )
490:                   C( K1, L2 ) = X( 2, 1 )
491: *
492:                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
493: *
494:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
495:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
496:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
497: *
498:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
499:                   SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
500:                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
501: *
502:                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
503:                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
504:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
505: *
506:                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
507:                   SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
508:                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
509: *
510:                   CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ),
511:      $                         LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
512:      $                         2, XNORM, IERR )
513:                   IF( IERR.NE.0 )
514:      $               INFO = 1
515: *
516:                   IF( SCALOC.NE.ONE ) THEN
517:                      DO 100 J = 1, N
518:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
519:   100                CONTINUE
520:                      SCALE = SCALE*SCALOC
521:                   END IF
522:                   C( K1, L1 ) = X( 1, 1 )
523:                   C( K1, L2 ) = X( 1, 2 )
524:                   C( K2, L1 ) = X( 2, 1 )
525:                   C( K2, L2 ) = X( 2, 2 )
526:                END IF
527: *
528:   110       CONTINUE
529:   120    CONTINUE
530: *
531:       ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
532: *
533: *        Solve    A'*X + ISGN*X*B' = scale*C.
534: *
535: *        The (K,L)th block of X is determined starting from
536: *        top-right corner column by column by
537: *
538: *           A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L)
539: *
540: *        Where
541: *                     K-1                          N
542: *            R(K,L) = SUM [A(I,K)'*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)'].
543: *                     I=1                        J=L+1
544: *
545: *        Start column loop (index = L)
546: *        L1 (L2): column index of the first (last) row of X(K,L)
547: *
548:          LNEXT = N
549:          DO 180 L = N, 1, -1
550:             IF( L.GT.LNEXT )
551:      $         GO TO 180
552:             IF( L.EQ.1 ) THEN
553:                L1 = L
554:                L2 = L
555:             ELSE
556:                IF( B( L, L-1 ).NE.ZERO ) THEN
557:                   L1 = L - 1
558:                   L2 = L
559:                   LNEXT = L - 2
560:                ELSE
561:                   L1 = L
562:                   L2 = L
563:                   LNEXT = L - 1
564:                END IF
565:             END IF
566: *
567: *           Start row loop (index = K)
568: *           K1 (K2): row index of the first (last) row of X(K,L)
569: *
570:             KNEXT = 1
571:             DO 170 K = 1, M
572:                IF( K.LT.KNEXT )
573:      $            GO TO 170
574:                IF( K.EQ.M ) THEN
575:                   K1 = K
576:                   K2 = K
577:                ELSE
578:                   IF( A( K+1, K ).NE.ZERO ) THEN
579:                      K1 = K
580:                      K2 = K + 1
581:                      KNEXT = K + 2
582:                   ELSE
583:                      K1 = K
584:                      K2 = K
585:                      KNEXT = K + 1
586:                   END IF
587:                END IF
588: *
589:                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
590:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
591:                   SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
592:      $                   B( L1, MIN( L1+1, N ) ), LDB )
593:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
594:                   SCALOC = ONE
595: *
596:                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
597:                   DA11 = ABS( A11 )
598:                   IF( DA11.LE.SMIN ) THEN
599:                      A11 = SMIN
600:                      DA11 = SMIN
601:                      INFO = 1
602:                   END IF
603:                   DB = ABS( VEC( 1, 1 ) )
604:                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
605:                      IF( DB.GT.BIGNUM*DA11 )
606:      $                  SCALOC = ONE / DB
607:                   END IF
608:                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
609: *
610:                   IF( SCALOC.NE.ONE ) THEN
611:                      DO 130 J = 1, N
612:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
613:   130                CONTINUE
614:                      SCALE = SCALE*SCALOC
615:                   END IF
616:                   C( K1, L1 ) = X( 1, 1 )
617: *
618:                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
619: *
620:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
621:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
622:      $                   B( L1, MIN( L2+1, N ) ), LDB )
623:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
624: *
625:                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
626:                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
627:      $                   B( L1, MIN( L2+1, N ) ), LDB )
628:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
629: *
630:                   CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
631:      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
632:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
633:                   IF( IERR.NE.0 )
634:      $               INFO = 1
635: *
636:                   IF( SCALOC.NE.ONE ) THEN
637:                      DO 140 J = 1, N
638:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
639:   140                CONTINUE
640:                      SCALE = SCALE*SCALOC
641:                   END IF
642:                   C( K1, L1 ) = X( 1, 1 )
643:                   C( K2, L1 ) = X( 2, 1 )
644: *
645:                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
646: *
647:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
648:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
649:      $                   B( L1, MIN( L2+1, N ) ), LDB )
650:                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
651: *
652:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
653:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
654:      $                   B( L2, MIN( L2+1, N ) ), LDB )
655:                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
656: *
657:                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
658:      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
659:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
660:                   IF( IERR.NE.0 )
661:      $               INFO = 1
662: *
663:                   IF( SCALOC.NE.ONE ) THEN
664:                      DO 150 J = 1, N
665:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
666:   150                CONTINUE
667:                      SCALE = SCALE*SCALOC
668:                   END IF
669:                   C( K1, L1 ) = X( 1, 1 )
670:                   C( K1, L2 ) = X( 2, 1 )
671: *
672:                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
673: *
674:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
675:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
676:      $                   B( L1, MIN( L2+1, N ) ), LDB )
677:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
678: *
679:                   SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
680:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
681:      $                   B( L2, MIN( L2+1, N ) ), LDB )
682:                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
683: *
684:                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
685:                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
686:      $                   B( L1, MIN( L2+1, N ) ), LDB )
687:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
688: *
689:                   SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
690:                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
691:      $                   B( L2, MIN( L2+1, N ) ), LDB )
692:                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
693: *
694:                   CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
695:      $                         LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
696:      $                         2, XNORM, IERR )
697:                   IF( IERR.NE.0 )
698:      $               INFO = 1
699: *
700:                   IF( SCALOC.NE.ONE ) THEN
701:                      DO 160 J = 1, N
702:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
703:   160                CONTINUE
704:                      SCALE = SCALE*SCALOC
705:                   END IF
706:                   C( K1, L1 ) = X( 1, 1 )
707:                   C( K1, L2 ) = X( 1, 2 )
708:                   C( K2, L1 ) = X( 2, 1 )
709:                   C( K2, L2 ) = X( 2, 2 )
710:                END IF
711: *
712:   170       CONTINUE
713:   180    CONTINUE
714: *
715:       ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
716: *
717: *        Solve    A*X + ISGN*X*B' = scale*C.
718: *
719: *        The (K,L)th block of X is determined starting from
720: *        bottom-right corner column by column by
721: *
722: *            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L)
723: *
724: *        Where
725: *                      M                          N
726: *            R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)'].
727: *                    I=K+1                      J=L+1
728: *
729: *        Start column loop (index = L)
730: *        L1 (L2): column index of the first (last) row of X(K,L)
731: *
732:          LNEXT = N
733:          DO 240 L = N, 1, -1
734:             IF( L.GT.LNEXT )
735:      $         GO TO 240
736:             IF( L.EQ.1 ) THEN
737:                L1 = L
738:                L2 = L
739:             ELSE
740:                IF( B( L, L-1 ).NE.ZERO ) THEN
741:                   L1 = L - 1
742:                   L2 = L
743:                   LNEXT = L - 2
744:                ELSE
745:                   L1 = L
746:                   L2 = L
747:                   LNEXT = L - 1
748:                END IF
749:             END IF
750: *
751: *           Start row loop (index = K)
752: *           K1 (K2): row index of the first (last) row of X(K,L)
753: *
754:             KNEXT = M
755:             DO 230 K = M, 1, -1
756:                IF( K.GT.KNEXT )
757:      $            GO TO 230
758:                IF( K.EQ.1 ) THEN
759:                   K1 = K
760:                   K2 = K
761:                ELSE
762:                   IF( A( K, K-1 ).NE.ZERO ) THEN
763:                      K1 = K - 1
764:                      K2 = K
765:                      KNEXT = K - 2
766:                   ELSE
767:                      K1 = K
768:                      K2 = K
769:                      KNEXT = K - 1
770:                   END IF
771:                END IF
772: *
773:                IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
774:                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
775:      $                   C( MIN( K1+1, M ), L1 ), 1 )
776:                   SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
777:      $                   B( L1, MIN( L1+1, N ) ), LDB )
778:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
779:                   SCALOC = ONE
780: *
781:                   A11 = A( K1, K1 ) + SGN*B( L1, L1 )
782:                   DA11 = ABS( A11 )
783:                   IF( DA11.LE.SMIN ) THEN
784:                      A11 = SMIN
785:                      DA11 = SMIN
786:                      INFO = 1
787:                   END IF
788:                   DB = ABS( VEC( 1, 1 ) )
789:                   IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
790:                      IF( DB.GT.BIGNUM*DA11 )
791:      $                  SCALOC = ONE / DB
792:                   END IF
793:                   X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
794: *
795:                   IF( SCALOC.NE.ONE ) THEN
796:                      DO 190 J = 1, N
797:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
798:   190                CONTINUE
799:                      SCALE = SCALE*SCALOC
800:                   END IF
801:                   C( K1, L1 ) = X( 1, 1 )
802: *
803:                ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
804: *
805:                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
806:      $                   C( MIN( K2+1, M ), L1 ), 1 )
807:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
808:      $                   B( L1, MIN( L2+1, N ) ), LDB )
809:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
810: *
811:                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
812:      $                   C( MIN( K2+1, M ), L1 ), 1 )
813:                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
814:      $                   B( L1, MIN( L2+1, N ) ), LDB )
815:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
816: *
817:                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
818:      $                         LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
819:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
820:                   IF( IERR.NE.0 )
821:      $               INFO = 1
822: *
823:                   IF( SCALOC.NE.ONE ) THEN
824:                      DO 200 J = 1, N
825:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
826:   200                CONTINUE
827:                      SCALE = SCALE*SCALOC
828:                   END IF
829:                   C( K1, L1 ) = X( 1, 1 )
830:                   C( K2, L1 ) = X( 2, 1 )
831: *
832:                ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
833: *
834:                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
835:      $                   C( MIN( K1+1, M ), L1 ), 1 )
836:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
837:      $                   B( L1, MIN( L2+1, N ) ), LDB )
838:                   VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
839: *
840:                   SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
841:      $                   C( MIN( K1+1, M ), L2 ), 1 )
842:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
843:      $                   B( L2, MIN( L2+1, N ) ), LDB )
844:                   VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
845: *
846:                   CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
847:      $                         LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
848:      $                         ZERO, X, 2, SCALOC, XNORM, IERR )
849:                   IF( IERR.NE.0 )
850:      $               INFO = 1
851: *
852:                   IF( SCALOC.NE.ONE ) THEN
853:                      DO 210 J = 1, N
854:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
855:   210                CONTINUE
856:                      SCALE = SCALE*SCALOC
857:                   END IF
858:                   C( K1, L1 ) = X( 1, 1 )
859:                   C( K1, L2 ) = X( 2, 1 )
860: *
861:                ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
862: *
863:                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
864:      $                   C( MIN( K2+1, M ), L1 ), 1 )
865:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
866:      $                   B( L1, MIN( L2+1, N ) ), LDB )
867:                   VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
868: *
869:                   SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
870:      $                   C( MIN( K2+1, M ), L2 ), 1 )
871:                   SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
872:      $                   B( L2, MIN( L2+1, N ) ), LDB )
873:                   VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
874: *
875:                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
876:      $                   C( MIN( K2+1, M ), L1 ), 1 )
877:                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
878:      $                   B( L1, MIN( L2+1, N ) ), LDB )
879:                   VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
880: *
881:                   SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
882:      $                   C( MIN( K2+1, M ), L2 ), 1 )
883:                   SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
884:      $                   B( L2, MIN( L2+1, N ) ), LDB )
885:                   VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
886: *
887:                   CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
888:      $                         LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
889:      $                         2, XNORM, IERR )
890:                   IF( IERR.NE.0 )
891:      $               INFO = 1
892: *
893:                   IF( SCALOC.NE.ONE ) THEN
894:                      DO 220 J = 1, N
895:                         CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
896:   220                CONTINUE
897:                      SCALE = SCALE*SCALOC
898:                   END IF
899:                   C( K1, L1 ) = X( 1, 1 )
900:                   C( K1, L2 ) = X( 1, 2 )
901:                   C( K2, L1 ) = X( 2, 1 )
902:                   C( K2, L2 ) = X( 2, 2 )
903:                END IF
904: *
905:   230       CONTINUE
906:   240    CONTINUE
907: *
908:       END IF
909: *
910:       RETURN
911: *
912: *     End of DTRSYL
913: *
914:       END
915: