1: SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
2: $ LDC, SCALE, INFO )
3: *
4: * -- LAPACK routine (version 3.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
8: *
9: * .. Scalar Arguments ..
10: CHARACTER TRANA, TRANB
11: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
12: DOUBLE PRECISION SCALE
13: * ..
14: * .. Array Arguments ..
15: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
16: * ..
17: *
18: * Purpose
19: * =======
20: *
21: * DTRSYL solves the real Sylvester matrix equation:
22: *
23: * op(A)*X + X*op(B) = scale*C or
24: * op(A)*X - X*op(B) = scale*C,
25: *
26: * where op(A) = A or A**T, and A and B are both upper quasi-
27: * triangular. A is M-by-M and B is N-by-N; the right hand side C and
28: * the solution X are M-by-N; and scale is an output scale factor, set
29: * <= 1 to avoid overflow in X.
30: *
31: * A and B must be in Schur canonical form (as returned by DHSEQR), that
32: * is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
33: * each 2-by-2 diagonal block has its diagonal elements equal and its
34: * off-diagonal elements of opposite sign.
35: *
36: * Arguments
37: * =========
38: *
39: * TRANA (input) CHARACTER*1
40: * Specifies the option op(A):
41: * = 'N': op(A) = A (No transpose)
42: * = 'T': op(A) = A**T (Transpose)
43: * = 'C': op(A) = A**H (Conjugate transpose = Transpose)
44: *
45: * TRANB (input) CHARACTER*1
46: * Specifies the option op(B):
47: * = 'N': op(B) = B (No transpose)
48: * = 'T': op(B) = B**T (Transpose)
49: * = 'C': op(B) = B**H (Conjugate transpose = Transpose)
50: *
51: * ISGN (input) INTEGER
52: * Specifies the sign in the equation:
53: * = +1: solve op(A)*X + X*op(B) = scale*C
54: * = -1: solve op(A)*X - X*op(B) = scale*C
55: *
56: * M (input) INTEGER
57: * The order of the matrix A, and the number of rows in the
58: * matrices X and C. M >= 0.
59: *
60: * N (input) INTEGER
61: * The order of the matrix B, and the number of columns in the
62: * matrices X and C. N >= 0.
63: *
64: * A (input) DOUBLE PRECISION array, dimension (LDA,M)
65: * The upper quasi-triangular matrix A, in Schur canonical form.
66: *
67: * LDA (input) INTEGER
68: * The leading dimension of the array A. LDA >= max(1,M).
69: *
70: * B (input) DOUBLE PRECISION array, dimension (LDB,N)
71: * The upper quasi-triangular matrix B, in Schur canonical form.
72: *
73: * LDB (input) INTEGER
74: * The leading dimension of the array B. LDB >= max(1,N).
75: *
76: * C (input/output) DOUBLE PRECISION array, dimension (LDC,N)
77: * On entry, the M-by-N right hand side matrix C.
78: * On exit, C is overwritten by the solution matrix X.
79: *
80: * LDC (input) INTEGER
81: * The leading dimension of the array C. LDC >= max(1,M)
82: *
83: * SCALE (output) DOUBLE PRECISION
84: * The scale factor, scale, set <= 1 to avoid overflow in X.
85: *
86: * INFO (output) INTEGER
87: * = 0: successful exit
88: * < 0: if INFO = -i, the i-th argument had an illegal value
89: * = 1: A and B have common or very close eigenvalues; perturbed
90: * values were used to solve the equation (but the matrices
91: * A and B are unchanged).
92: *
93: * =====================================================================
94: *
95: * .. Parameters ..
96: DOUBLE PRECISION ZERO, ONE
97: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
98: * ..
99: * .. 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrsyl.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack