1: *> \brief \b DTRSYL
2: *
3: * =========== DOCUMENTATION ===========
4: *
5: * Online html documentation available at
6: * http://www.netlib.org/lapack/explore-html/
7: *
8: *> \htmlonly
9: *> Download DTRSYL + dependencies
10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsyl.f">
11: *> [TGZ]</a>
12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsyl.f">
13: *> [ZIP]</a>
14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsyl.f">
15: *> [TXT]</a>
16: *> \endhtmlonly
17: *
18: * Definition:
19: * ===========
20: *
21: * SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
22: * LDC, SCALE, INFO )
23: *
24: * .. Scalar Arguments ..
25: * CHARACTER TRANA, TRANB
26: * INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
27: * DOUBLE PRECISION SCALE
28: * ..
29: * .. Array Arguments ..
30: * DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
31: * ..
32: *
33: *
34: *> \par Purpose:
35: * =============
36: *>
37: *> \verbatim
38: *>
39: *> DTRSYL solves the real Sylvester matrix equation:
40: *>
41: *> op(A)*X + X*op(B) = scale*C or
42: *> op(A)*X - X*op(B) = scale*C,
43: *>
44: *> where op(A) = A or A**T, and A and B are both upper quasi-
45: *> triangular. A is M-by-M and B is N-by-N; the right hand side C and
46: *> the solution X are M-by-N; and scale is an output scale factor, set
47: *> <= 1 to avoid overflow in X.
48: *>
49: *> A and B must be in Schur canonical form (as returned by DHSEQR), that
50: *> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
51: *> each 2-by-2 diagonal block has its diagonal elements equal and its
52: *> off-diagonal elements of opposite sign.
53: *> \endverbatim
54: *
55: * Arguments:
56: * ==========
57: *
58: *> \param[in] TRANA
59: *> \verbatim
60: *> TRANA is CHARACTER*1
61: *> Specifies the option op(A):
62: *> = 'N': op(A) = A (No transpose)
63: *> = 'T': op(A) = A**T (Transpose)
64: *> = 'C': op(A) = A**H (Conjugate transpose = Transpose)
65: *> \endverbatim
66: *>
67: *> \param[in] TRANB
68: *> \verbatim
69: *> TRANB is CHARACTER*1
70: *> Specifies the option op(B):
71: *> = 'N': op(B) = B (No transpose)
72: *> = 'T': op(B) = B**T (Transpose)
73: *> = 'C': op(B) = B**H (Conjugate transpose = Transpose)
74: *> \endverbatim
75: *>
76: *> \param[in] ISGN
77: *> \verbatim
78: *> ISGN is INTEGER
79: *> Specifies the sign in the equation:
80: *> = +1: solve op(A)*X + X*op(B) = scale*C
81: *> = -1: solve op(A)*X - X*op(B) = scale*C
82: *> \endverbatim
83: *>
84: *> \param[in] M
85: *> \verbatim
86: *> M is INTEGER
87: *> The order of the matrix A, and the number of rows in the
88: *> matrices X and C. M >= 0.
89: *> \endverbatim
90: *>
91: *> \param[in] N
92: *> \verbatim
93: *> N is INTEGER
94: *> The order of the matrix B, and the number of columns in the
95: *> matrices X and C. N >= 0.
96: *> \endverbatim
97: *>
98: *> \param[in] A
99: *> \verbatim
100: *> A is DOUBLE PRECISION array, dimension (LDA,M)
101: *> The upper quasi-triangular matrix A, in Schur canonical form.
102: *> \endverbatim
103: *>
104: *> \param[in] LDA
105: *> \verbatim
106: *> LDA is INTEGER
107: *> The leading dimension of the array A. LDA >= max(1,M).
108: *> \endverbatim
109: *>
110: *> \param[in] B
111: *> \verbatim
112: *> B is DOUBLE PRECISION array, dimension (LDB,N)
113: *> The upper quasi-triangular matrix B, in Schur canonical form.
114: *> \endverbatim
115: *>
116: *> \param[in] LDB
117: *> \verbatim
118: *> LDB is INTEGER
119: *> The leading dimension of the array B. LDB >= max(1,N).
120: *> \endverbatim
121: *>
122: *> \param[in,out] C
123: *> \verbatim
124: *> C is DOUBLE PRECISION array, dimension (LDC,N)
125: *> On entry, the M-by-N right hand side matrix C.
126: *> On exit, C is overwritten by the solution matrix X.
127: *> \endverbatim
128: *>
129: *> \param[in] LDC
130: *> \verbatim
131: *> LDC is INTEGER
132: *> The leading dimension of the array C. LDC >= max(1,M)
133: *> \endverbatim
134: *>
135: *> \param[out] SCALE
136: *> \verbatim
137: *> SCALE is DOUBLE PRECISION
138: *> The scale factor, scale, set <= 1 to avoid overflow in X.
139: *> \endverbatim
140: *>
141: *> \param[out] INFO
142: *> \verbatim
143: *> INFO is INTEGER
144: *> = 0: successful exit
145: *> < 0: if INFO = -i, the i-th argument had an illegal value
146: *> = 1: A and B have common or very close eigenvalues; perturbed
147: *> values were used to solve the equation (but the matrices
148: *> A and B are unchanged).
149: *> \endverbatim
150: *
151: * Authors:
152: * ========
153: *
154: *> \author Univ. of Tennessee
155: *> \author Univ. of California Berkeley
156: *> \author Univ. of Colorado Denver
157: *> \author NAG Ltd.
158: *
159: *> \ingroup doubleSYcomputational
160: *
161: * =====================================================================
162: SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
163: $ LDC, SCALE, INFO )
164: *
165: * -- LAPACK computational routine --
166: * -- LAPACK is a software package provided by Univ. of Tennessee, --
167: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168: *
169: * .. Scalar Arguments ..
170: CHARACTER TRANA, TRANB
171: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
172: DOUBLE PRECISION SCALE
173: * ..
174: * .. Array Arguments ..
175: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
176: * ..
177: *
178: * =====================================================================
179: *
180: * .. Parameters ..
181: DOUBLE PRECISION ZERO, ONE
182: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
183: * ..
184: * .. Local Scalars ..
185: LOGICAL NOTRNA, NOTRNB
186: INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
187: DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
188: $ SMLNUM, SUML, SUMR, XNORM
189: * ..
190: * .. Local Arrays ..
191: DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
192: * ..
193: * .. External Functions ..
194: LOGICAL LSAME
195: DOUBLE PRECISION DDOT, DLAMCH, DLANGE
196: EXTERNAL LSAME, DDOT, DLAMCH, DLANGE
197: * ..
198: * .. External Subroutines ..
199: EXTERNAL DLABAD, DLALN2, DLASY2, DSCAL, XERBLA
200: * ..
201: * .. Intrinsic Functions ..
202: INTRINSIC ABS, DBLE, MAX, MIN
203: * ..
204: * .. Executable Statements ..
205: *
206: * Decode and Test input parameters
207: *
208: NOTRNA = LSAME( TRANA, 'N' )
209: NOTRNB = LSAME( TRANB, 'N' )
210: *
211: INFO = 0
212: IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT.
213: $ LSAME( TRANA, 'C' ) ) THEN
214: INFO = -1
215: ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT.
216: $ LSAME( TRANB, 'C' ) ) THEN
217: INFO = -2
218: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
219: INFO = -3
220: ELSE IF( M.LT.0 ) THEN
221: INFO = -4
222: ELSE IF( N.LT.0 ) THEN
223: INFO = -5
224: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
225: INFO = -7
226: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
227: INFO = -9
228: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
229: INFO = -11
230: END IF
231: IF( INFO.NE.0 ) THEN
232: CALL XERBLA( 'DTRSYL', -INFO )
233: RETURN
234: END IF
235: *
236: * Quick return if possible
237: *
238: SCALE = ONE
239: IF( M.EQ.0 .OR. N.EQ.0 )
240: $ RETURN
241: *
242: * Set constants to control overflow
243: *
244: EPS = DLAMCH( 'P' )
245: SMLNUM = DLAMCH( 'S' )
246: BIGNUM = ONE / SMLNUM
247: CALL DLABAD( SMLNUM, BIGNUM )
248: SMLNUM = SMLNUM*DBLE( M*N ) / EPS
249: BIGNUM = ONE / SMLNUM
250: *
251: SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ),
252: $ EPS*DLANGE( 'M', N, N, B, LDB, DUM ) )
253: *
254: SGN = ISGN
255: *
256: IF( NOTRNA .AND. NOTRNB ) THEN
257: *
258: * Solve A*X + ISGN*X*B = scale*C.
259: *
260: * The (K,L)th block of X is determined starting from
261: * bottom-left corner column by column by
262: *
263: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
264: *
265: * Where
266: * M L-1
267: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
268: * I=K+1 J=1
269: *
270: * Start column loop (index = L)
271: * L1 (L2) : column index of the first (first) row of X(K,L).
272: *
273: LNEXT = 1
274: DO 60 L = 1, N
275: IF( L.LT.LNEXT )
276: $ GO TO 60
277: IF( L.EQ.N ) THEN
278: L1 = L
279: L2 = L
280: ELSE
281: IF( B( L+1, L ).NE.ZERO ) THEN
282: L1 = L
283: L2 = L + 1
284: LNEXT = L + 2
285: ELSE
286: L1 = L
287: L2 = L
288: LNEXT = L + 1
289: END IF
290: END IF
291: *
292: * Start row loop (index = K)
293: * K1 (K2): row index of the first (last) row of X(K,L).
294: *
295: KNEXT = M
296: DO 50 K = M, 1, -1
297: IF( K.GT.KNEXT )
298: $ GO TO 50
299: IF( K.EQ.1 ) THEN
300: K1 = K
301: K2 = K
302: ELSE
303: IF( A( K, K-1 ).NE.ZERO ) THEN
304: K1 = K - 1
305: K2 = K
306: KNEXT = K - 2
307: ELSE
308: K1 = K
309: K2 = K
310: KNEXT = K - 1
311: END IF
312: END IF
313: *
314: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
315: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
316: $ C( MIN( K1+1, M ), L1 ), 1 )
317: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
318: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
319: SCALOC = ONE
320: *
321: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
322: DA11 = ABS( A11 )
323: IF( DA11.LE.SMIN ) THEN
324: A11 = SMIN
325: DA11 = SMIN
326: INFO = 1
327: END IF
328: DB = ABS( VEC( 1, 1 ) )
329: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
330: IF( DB.GT.BIGNUM*DA11 )
331: $ SCALOC = ONE / DB
332: END IF
333: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
334: *
335: IF( SCALOC.NE.ONE ) THEN
336: DO 10 J = 1, N
337: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
338: 10 CONTINUE
339: SCALE = SCALE*SCALOC
340: END IF
341: C( K1, L1 ) = X( 1, 1 )
342: *
343: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
344: *
345: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
346: $ C( MIN( K2+1, M ), L1 ), 1 )
347: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
348: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
349: *
350: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
351: $ C( MIN( K2+1, M ), L1 ), 1 )
352: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
353: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
354: *
355: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
356: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
357: $ ZERO, X, 2, SCALOC, XNORM, IERR )
358: IF( IERR.NE.0 )
359: $ INFO = 1
360: *
361: IF( SCALOC.NE.ONE ) THEN
362: DO 20 J = 1, N
363: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
364: 20 CONTINUE
365: SCALE = SCALE*SCALOC
366: END IF
367: C( K1, L1 ) = X( 1, 1 )
368: C( K2, L1 ) = X( 2, 1 )
369: *
370: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
371: *
372: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
373: $ C( MIN( K1+1, M ), L1 ), 1 )
374: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
375: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
376: *
377: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
378: $ C( MIN( K1+1, M ), L2 ), 1 )
379: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
380: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
381: *
382: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
383: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
384: $ ZERO, X, 2, SCALOC, XNORM, IERR )
385: IF( IERR.NE.0 )
386: $ INFO = 1
387: *
388: IF( SCALOC.NE.ONE ) THEN
389: DO 30 J = 1, N
390: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
391: 30 CONTINUE
392: SCALE = SCALE*SCALOC
393: END IF
394: C( K1, L1 ) = X( 1, 1 )
395: C( K1, L2 ) = X( 2, 1 )
396: *
397: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
398: *
399: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
400: $ C( MIN( K2+1, M ), L1 ), 1 )
401: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
402: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
403: *
404: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
405: $ C( MIN( K2+1, M ), L2 ), 1 )
406: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
407: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
408: *
409: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
410: $ C( MIN( K2+1, M ), L1 ), 1 )
411: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
412: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
413: *
414: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
415: $ C( MIN( K2+1, M ), L2 ), 1 )
416: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
417: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
418: *
419: CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2,
420: $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC,
421: $ 2, SCALOC, X, 2, XNORM, IERR )
422: IF( IERR.NE.0 )
423: $ INFO = 1
424: *
425: IF( SCALOC.NE.ONE ) THEN
426: DO 40 J = 1, N
427: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
428: 40 CONTINUE
429: SCALE = SCALE*SCALOC
430: END IF
431: C( K1, L1 ) = X( 1, 1 )
432: C( K1, L2 ) = X( 1, 2 )
433: C( K2, L1 ) = X( 2, 1 )
434: C( K2, L2 ) = X( 2, 2 )
435: END IF
436: *
437: 50 CONTINUE
438: *
439: 60 CONTINUE
440: *
441: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
442: *
443: * Solve A**T *X + ISGN*X*B = scale*C.
444: *
445: * The (K,L)th block of X is determined starting from
446: * upper-left corner column by column by
447: *
448: * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
449: *
450: * Where
451: * K-1 T L-1
452: * R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
453: * I=1 J=1
454: *
455: * Start column loop (index = L)
456: * L1 (L2): column index of the first (last) row of X(K,L)
457: *
458: LNEXT = 1
459: DO 120 L = 1, N
460: IF( L.LT.LNEXT )
461: $ GO TO 120
462: IF( L.EQ.N ) THEN
463: L1 = L
464: L2 = L
465: ELSE
466: IF( B( L+1, L ).NE.ZERO ) THEN
467: L1 = L
468: L2 = L + 1
469: LNEXT = L + 2
470: ELSE
471: L1 = L
472: L2 = L
473: LNEXT = L + 1
474: END IF
475: END IF
476: *
477: * Start row loop (index = K)
478: * K1 (K2): row index of the first (last) row of X(K,L)
479: *
480: KNEXT = 1
481: DO 110 K = 1, M
482: IF( K.LT.KNEXT )
483: $ GO TO 110
484: IF( K.EQ.M ) THEN
485: K1 = K
486: K2 = K
487: ELSE
488: IF( A( K+1, K ).NE.ZERO ) THEN
489: K1 = K
490: K2 = K + 1
491: KNEXT = K + 2
492: ELSE
493: K1 = K
494: K2 = K
495: KNEXT = K + 1
496: END IF
497: END IF
498: *
499: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
500: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
501: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
502: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
503: SCALOC = ONE
504: *
505: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
506: DA11 = ABS( A11 )
507: IF( DA11.LE.SMIN ) THEN
508: A11 = SMIN
509: DA11 = SMIN
510: INFO = 1
511: END IF
512: DB = ABS( VEC( 1, 1 ) )
513: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
514: IF( DB.GT.BIGNUM*DA11 )
515: $ SCALOC = ONE / DB
516: END IF
517: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
518: *
519: IF( SCALOC.NE.ONE ) THEN
520: DO 70 J = 1, N
521: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
522: 70 CONTINUE
523: SCALE = SCALE*SCALOC
524: END IF
525: C( K1, L1 ) = X( 1, 1 )
526: *
527: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
528: *
529: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
530: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
531: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
532: *
533: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
534: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
535: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
536: *
537: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
538: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
539: $ ZERO, X, 2, SCALOC, XNORM, IERR )
540: IF( IERR.NE.0 )
541: $ INFO = 1
542: *
543: IF( SCALOC.NE.ONE ) THEN
544: DO 80 J = 1, N
545: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
546: 80 CONTINUE
547: SCALE = SCALE*SCALOC
548: END IF
549: C( K1, L1 ) = X( 1, 1 )
550: C( K2, L1 ) = X( 2, 1 )
551: *
552: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
553: *
554: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
555: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
556: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
557: *
558: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
559: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
560: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
561: *
562: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
563: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
564: $ ZERO, X, 2, SCALOC, XNORM, IERR )
565: IF( IERR.NE.0 )
566: $ INFO = 1
567: *
568: IF( SCALOC.NE.ONE ) THEN
569: DO 90 J = 1, N
570: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
571: 90 CONTINUE
572: SCALE = SCALE*SCALOC
573: END IF
574: C( K1, L1 ) = X( 1, 1 )
575: C( K1, L2 ) = X( 2, 1 )
576: *
577: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
578: *
579: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
580: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
581: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
582: *
583: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
584: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
585: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
586: *
587: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
588: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
589: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
590: *
591: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
592: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
593: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
594: *
595: CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ),
596: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
597: $ 2, XNORM, IERR )
598: IF( IERR.NE.0 )
599: $ INFO = 1
600: *
601: IF( SCALOC.NE.ONE ) THEN
602: DO 100 J = 1, N
603: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
604: 100 CONTINUE
605: SCALE = SCALE*SCALOC
606: END IF
607: C( K1, L1 ) = X( 1, 1 )
608: C( K1, L2 ) = X( 1, 2 )
609: C( K2, L1 ) = X( 2, 1 )
610: C( K2, L2 ) = X( 2, 2 )
611: END IF
612: *
613: 110 CONTINUE
614: 120 CONTINUE
615: *
616: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
617: *
618: * Solve A**T*X + ISGN*X*B**T = scale*C.
619: *
620: * The (K,L)th block of X is determined starting from
621: * top-right corner column by column by
622: *
623: * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
624: *
625: * Where
626: * K-1 N
627: * R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
628: * I=1 J=L+1
629: *
630: * Start column loop (index = L)
631: * L1 (L2): column index of the first (last) row of X(K,L)
632: *
633: LNEXT = N
634: DO 180 L = N, 1, -1
635: IF( L.GT.LNEXT )
636: $ GO TO 180
637: IF( L.EQ.1 ) THEN
638: L1 = L
639: L2 = L
640: ELSE
641: IF( B( L, L-1 ).NE.ZERO ) THEN
642: L1 = L - 1
643: L2 = L
644: LNEXT = L - 2
645: ELSE
646: L1 = L
647: L2 = L
648: LNEXT = L - 1
649: END IF
650: END IF
651: *
652: * Start row loop (index = K)
653: * K1 (K2): row index of the first (last) row of X(K,L)
654: *
655: KNEXT = 1
656: DO 170 K = 1, M
657: IF( K.LT.KNEXT )
658: $ GO TO 170
659: IF( K.EQ.M ) THEN
660: K1 = K
661: K2 = K
662: ELSE
663: IF( A( K+1, K ).NE.ZERO ) THEN
664: K1 = K
665: K2 = K + 1
666: KNEXT = K + 2
667: ELSE
668: K1 = K
669: K2 = K
670: KNEXT = K + 1
671: END IF
672: END IF
673: *
674: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
675: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
676: SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
677: $ B( L1, MIN( L1+1, N ) ), LDB )
678: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
679: SCALOC = ONE
680: *
681: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
682: DA11 = ABS( A11 )
683: IF( DA11.LE.SMIN ) THEN
684: A11 = SMIN
685: DA11 = SMIN
686: INFO = 1
687: END IF
688: DB = ABS( VEC( 1, 1 ) )
689: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
690: IF( DB.GT.BIGNUM*DA11 )
691: $ SCALOC = ONE / DB
692: END IF
693: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
694: *
695: IF( SCALOC.NE.ONE ) THEN
696: DO 130 J = 1, N
697: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
698: 130 CONTINUE
699: SCALE = SCALE*SCALOC
700: END IF
701: C( K1, L1 ) = X( 1, 1 )
702: *
703: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
704: *
705: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
706: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
707: $ B( L1, MIN( L2+1, N ) ), LDB )
708: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
709: *
710: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
711: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
712: $ B( L1, MIN( L2+1, N ) ), LDB )
713: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
714: *
715: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
716: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
717: $ ZERO, X, 2, SCALOC, XNORM, IERR )
718: IF( IERR.NE.0 )
719: $ INFO = 1
720: *
721: IF( SCALOC.NE.ONE ) THEN
722: DO 140 J = 1, N
723: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
724: 140 CONTINUE
725: SCALE = SCALE*SCALOC
726: END IF
727: C( K1, L1 ) = X( 1, 1 )
728: C( K2, L1 ) = X( 2, 1 )
729: *
730: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
731: *
732: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
733: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
734: $ B( L1, MIN( L2+1, N ) ), LDB )
735: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
736: *
737: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
738: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
739: $ B( L2, MIN( L2+1, N ) ), LDB )
740: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
741: *
742: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
743: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
744: $ ZERO, X, 2, SCALOC, XNORM, IERR )
745: IF( IERR.NE.0 )
746: $ INFO = 1
747: *
748: IF( SCALOC.NE.ONE ) THEN
749: DO 150 J = 1, N
750: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
751: 150 CONTINUE
752: SCALE = SCALE*SCALOC
753: END IF
754: C( K1, L1 ) = X( 1, 1 )
755: C( K1, L2 ) = X( 2, 1 )
756: *
757: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
758: *
759: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
760: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
761: $ B( L1, MIN( L2+1, N ) ), LDB )
762: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
763: *
764: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
765: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
766: $ B( L2, MIN( L2+1, N ) ), LDB )
767: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
768: *
769: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
770: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
771: $ B( L1, MIN( L2+1, N ) ), LDB )
772: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
773: *
774: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
775: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
776: $ B( L2, MIN( L2+1, N ) ), LDB )
777: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
778: *
779: CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
780: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
781: $ 2, XNORM, IERR )
782: IF( IERR.NE.0 )
783: $ INFO = 1
784: *
785: IF( SCALOC.NE.ONE ) THEN
786: DO 160 J = 1, N
787: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
788: 160 CONTINUE
789: SCALE = SCALE*SCALOC
790: END IF
791: C( K1, L1 ) = X( 1, 1 )
792: C( K1, L2 ) = X( 1, 2 )
793: C( K2, L1 ) = X( 2, 1 )
794: C( K2, L2 ) = X( 2, 2 )
795: END IF
796: *
797: 170 CONTINUE
798: 180 CONTINUE
799: *
800: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
801: *
802: * Solve A*X + ISGN*X*B**T = scale*C.
803: *
804: * The (K,L)th block of X is determined starting from
805: * bottom-right corner column by column by
806: *
807: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
808: *
809: * Where
810: * M N
811: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
812: * I=K+1 J=L+1
813: *
814: * Start column loop (index = L)
815: * L1 (L2): column index of the first (last) row of X(K,L)
816: *
817: LNEXT = N
818: DO 240 L = N, 1, -1
819: IF( L.GT.LNEXT )
820: $ GO TO 240
821: IF( L.EQ.1 ) THEN
822: L1 = L
823: L2 = L
824: ELSE
825: IF( B( L, L-1 ).NE.ZERO ) THEN
826: L1 = L - 1
827: L2 = L
828: LNEXT = L - 2
829: ELSE
830: L1 = L
831: L2 = L
832: LNEXT = L - 1
833: END IF
834: END IF
835: *
836: * Start row loop (index = K)
837: * K1 (K2): row index of the first (last) row of X(K,L)
838: *
839: KNEXT = M
840: DO 230 K = M, 1, -1
841: IF( K.GT.KNEXT )
842: $ GO TO 230
843: IF( K.EQ.1 ) THEN
844: K1 = K
845: K2 = K
846: ELSE
847: IF( A( K, K-1 ).NE.ZERO ) THEN
848: K1 = K - 1
849: K2 = K
850: KNEXT = K - 2
851: ELSE
852: K1 = K
853: K2 = K
854: KNEXT = K - 1
855: END IF
856: END IF
857: *
858: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
859: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
860: $ C( MIN( K1+1, M ), L1 ), 1 )
861: SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
862: $ B( L1, MIN( L1+1, N ) ), LDB )
863: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
864: SCALOC = ONE
865: *
866: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
867: DA11 = ABS( A11 )
868: IF( DA11.LE.SMIN ) THEN
869: A11 = SMIN
870: DA11 = SMIN
871: INFO = 1
872: END IF
873: DB = ABS( VEC( 1, 1 ) )
874: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
875: IF( DB.GT.BIGNUM*DA11 )
876: $ SCALOC = ONE / DB
877: END IF
878: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
879: *
880: IF( SCALOC.NE.ONE ) THEN
881: DO 190 J = 1, N
882: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
883: 190 CONTINUE
884: SCALE = SCALE*SCALOC
885: END IF
886: C( K1, L1 ) = X( 1, 1 )
887: *
888: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
889: *
890: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
891: $ C( MIN( K2+1, M ), L1 ), 1 )
892: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
893: $ B( L1, MIN( L2+1, N ) ), LDB )
894: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
895: *
896: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
897: $ C( MIN( K2+1, M ), L1 ), 1 )
898: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
899: $ B( L1, MIN( L2+1, N ) ), LDB )
900: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
901: *
902: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
903: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
904: $ ZERO, X, 2, SCALOC, XNORM, IERR )
905: IF( IERR.NE.0 )
906: $ INFO = 1
907: *
908: IF( SCALOC.NE.ONE ) THEN
909: DO 200 J = 1, N
910: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
911: 200 CONTINUE
912: SCALE = SCALE*SCALOC
913: END IF
914: C( K1, L1 ) = X( 1, 1 )
915: C( K2, L1 ) = X( 2, 1 )
916: *
917: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
918: *
919: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
920: $ C( MIN( K1+1, M ), L1 ), 1 )
921: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
922: $ B( L1, MIN( L2+1, N ) ), LDB )
923: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
924: *
925: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
926: $ C( MIN( K1+1, M ), L2 ), 1 )
927: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
928: $ B( L2, MIN( L2+1, N ) ), LDB )
929: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
930: *
931: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
932: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
933: $ ZERO, X, 2, SCALOC, XNORM, IERR )
934: IF( IERR.NE.0 )
935: $ INFO = 1
936: *
937: IF( SCALOC.NE.ONE ) THEN
938: DO 210 J = 1, N
939: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
940: 210 CONTINUE
941: SCALE = SCALE*SCALOC
942: END IF
943: C( K1, L1 ) = X( 1, 1 )
944: C( K1, L2 ) = X( 2, 1 )
945: *
946: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
947: *
948: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
949: $ C( MIN( K2+1, M ), L1 ), 1 )
950: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
951: $ B( L1, MIN( L2+1, N ) ), LDB )
952: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
953: *
954: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
955: $ C( MIN( K2+1, M ), L2 ), 1 )
956: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
957: $ B( L2, MIN( L2+1, N ) ), LDB )
958: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
959: *
960: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
961: $ C( MIN( K2+1, M ), L1 ), 1 )
962: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
963: $ B( L1, MIN( L2+1, N ) ), LDB )
964: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
965: *
966: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
967: $ C( MIN( K2+1, M ), L2 ), 1 )
968: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
969: $ B( L2, MIN( L2+1, N ) ), LDB )
970: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
971: *
972: CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
973: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
974: $ 2, XNORM, IERR )
975: IF( IERR.NE.0 )
976: $ INFO = 1
977: *
978: IF( SCALOC.NE.ONE ) THEN
979: DO 220 J = 1, N
980: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
981: 220 CONTINUE
982: SCALE = SCALE*SCALOC
983: END IF
984: C( K1, L1 ) = X( 1, 1 )
985: C( K1, L2 ) = X( 1, 2 )
986: C( K2, L1 ) = X( 2, 1 )
987: C( K2, L2 ) = X( 2, 2 )
988: END IF
989: *
990: 230 CONTINUE
991: 240 CONTINUE
992: *
993: END IF
994: *
995: RETURN
996: *
997: * End of DTRSYL
998: *
999: 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