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: *> \date December 2016
160: *
161: *> \ingroup doubleSYcomputational
162: *
163: * =====================================================================
164: SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
165: $ LDC, SCALE, INFO )
166: *
167: * -- LAPACK computational routine (version 3.7.0) --
168: * -- LAPACK is a software package provided by Univ. of Tennessee, --
169: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
170: * December 2016
171: *
172: * .. Scalar Arguments ..
173: CHARACTER TRANA, TRANB
174: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
175: DOUBLE PRECISION SCALE
176: * ..
177: * .. Array Arguments ..
178: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
179: * ..
180: *
181: * =====================================================================
182: *
183: * .. Parameters ..
184: DOUBLE PRECISION ZERO, ONE
185: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
186: * ..
187: * .. Local Scalars ..
188: LOGICAL NOTRNA, NOTRNB
189: INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
190: DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
191: $ SMLNUM, SUML, SUMR, XNORM
192: * ..
193: * .. Local Arrays ..
194: DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
195: * ..
196: * .. External Functions ..
197: LOGICAL LSAME
198: DOUBLE PRECISION DDOT, DLAMCH, DLANGE
199: EXTERNAL LSAME, DDOT, DLAMCH, DLANGE
200: * ..
201: * .. External Subroutines ..
202: EXTERNAL DLABAD, DLALN2, DLASY2, DSCAL, XERBLA
203: * ..
204: * .. Intrinsic Functions ..
205: INTRINSIC ABS, DBLE, MAX, MIN
206: * ..
207: * .. Executable Statements ..
208: *
209: * Decode and Test input parameters
210: *
211: NOTRNA = LSAME( TRANA, 'N' )
212: NOTRNB = LSAME( TRANB, 'N' )
213: *
214: INFO = 0
215: IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT.
216: $ LSAME( TRANA, 'C' ) ) THEN
217: INFO = -1
218: ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT.
219: $ LSAME( TRANB, 'C' ) ) THEN
220: INFO = -2
221: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
222: INFO = -3
223: ELSE IF( M.LT.0 ) THEN
224: INFO = -4
225: ELSE IF( N.LT.0 ) THEN
226: INFO = -5
227: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
228: INFO = -7
229: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
230: INFO = -9
231: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
232: INFO = -11
233: END IF
234: IF( INFO.NE.0 ) THEN
235: CALL XERBLA( 'DTRSYL', -INFO )
236: RETURN
237: END IF
238: *
239: * Quick return if possible
240: *
241: SCALE = ONE
242: IF( M.EQ.0 .OR. N.EQ.0 )
243: $ RETURN
244: *
245: * Set constants to control overflow
246: *
247: EPS = DLAMCH( 'P' )
248: SMLNUM = DLAMCH( 'S' )
249: BIGNUM = ONE / SMLNUM
250: CALL DLABAD( SMLNUM, BIGNUM )
251: SMLNUM = SMLNUM*DBLE( M*N ) / EPS
252: BIGNUM = ONE / SMLNUM
253: *
254: SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ),
255: $ EPS*DLANGE( 'M', N, N, B, LDB, DUM ) )
256: *
257: SGN = ISGN
258: *
259: IF( NOTRNA .AND. NOTRNB ) THEN
260: *
261: * Solve A*X + ISGN*X*B = scale*C.
262: *
263: * The (K,L)th block of X is determined starting from
264: * bottom-left corner column by column by
265: *
266: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
267: *
268: * Where
269: * M L-1
270: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
271: * I=K+1 J=1
272: *
273: * Start column loop (index = L)
274: * L1 (L2) : column index of the first (first) row of X(K,L).
275: *
276: LNEXT = 1
277: DO 60 L = 1, N
278: IF( L.LT.LNEXT )
279: $ GO TO 60
280: IF( L.EQ.N ) THEN
281: L1 = L
282: L2 = L
283: ELSE
284: IF( B( L+1, L ).NE.ZERO ) THEN
285: L1 = L
286: L2 = L + 1
287: LNEXT = L + 2
288: ELSE
289: L1 = L
290: L2 = L
291: LNEXT = L + 1
292: END IF
293: END IF
294: *
295: * Start row loop (index = K)
296: * K1 (K2): row index of the first (last) row of X(K,L).
297: *
298: KNEXT = M
299: DO 50 K = M, 1, -1
300: IF( K.GT.KNEXT )
301: $ GO TO 50
302: IF( K.EQ.1 ) THEN
303: K1 = K
304: K2 = K
305: ELSE
306: IF( A( K, K-1 ).NE.ZERO ) THEN
307: K1 = K - 1
308: K2 = K
309: KNEXT = K - 2
310: ELSE
311: K1 = K
312: K2 = K
313: KNEXT = K - 1
314: END IF
315: END IF
316: *
317: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
318: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
319: $ C( MIN( K1+1, M ), L1 ), 1 )
320: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
321: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
322: SCALOC = ONE
323: *
324: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
325: DA11 = ABS( A11 )
326: IF( DA11.LE.SMIN ) THEN
327: A11 = SMIN
328: DA11 = SMIN
329: INFO = 1
330: END IF
331: DB = ABS( VEC( 1, 1 ) )
332: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
333: IF( DB.GT.BIGNUM*DA11 )
334: $ SCALOC = ONE / DB
335: END IF
336: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
337: *
338: IF( SCALOC.NE.ONE ) THEN
339: DO 10 J = 1, N
340: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
341: 10 CONTINUE
342: SCALE = SCALE*SCALOC
343: END IF
344: C( K1, L1 ) = X( 1, 1 )
345: *
346: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
347: *
348: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
349: $ C( MIN( K2+1, M ), L1 ), 1 )
350: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
351: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
352: *
353: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
354: $ C( MIN( K2+1, M ), L1 ), 1 )
355: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
356: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
357: *
358: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
359: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
360: $ ZERO, X, 2, SCALOC, XNORM, IERR )
361: IF( IERR.NE.0 )
362: $ INFO = 1
363: *
364: IF( SCALOC.NE.ONE ) THEN
365: DO 20 J = 1, N
366: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
367: 20 CONTINUE
368: SCALE = SCALE*SCALOC
369: END IF
370: C( K1, L1 ) = X( 1, 1 )
371: C( K2, L1 ) = X( 2, 1 )
372: *
373: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
374: *
375: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
376: $ C( MIN( K1+1, M ), L1 ), 1 )
377: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
378: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
379: *
380: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
381: $ C( MIN( K1+1, M ), L2 ), 1 )
382: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
383: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
384: *
385: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
386: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
387: $ ZERO, X, 2, SCALOC, XNORM, IERR )
388: IF( IERR.NE.0 )
389: $ INFO = 1
390: *
391: IF( SCALOC.NE.ONE ) THEN
392: DO 30 J = 1, N
393: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
394: 30 CONTINUE
395: SCALE = SCALE*SCALOC
396: END IF
397: C( K1, L1 ) = X( 1, 1 )
398: C( K1, L2 ) = X( 2, 1 )
399: *
400: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
401: *
402: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
403: $ C( MIN( K2+1, M ), L1 ), 1 )
404: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
405: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
406: *
407: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
408: $ C( MIN( K2+1, M ), L2 ), 1 )
409: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
410: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
411: *
412: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
413: $ C( MIN( K2+1, M ), L1 ), 1 )
414: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
415: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
416: *
417: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
418: $ C( MIN( K2+1, M ), L2 ), 1 )
419: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
420: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
421: *
422: CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2,
423: $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC,
424: $ 2, SCALOC, X, 2, XNORM, IERR )
425: IF( IERR.NE.0 )
426: $ INFO = 1
427: *
428: IF( SCALOC.NE.ONE ) THEN
429: DO 40 J = 1, N
430: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
431: 40 CONTINUE
432: SCALE = SCALE*SCALOC
433: END IF
434: C( K1, L1 ) = X( 1, 1 )
435: C( K1, L2 ) = X( 1, 2 )
436: C( K2, L1 ) = X( 2, 1 )
437: C( K2, L2 ) = X( 2, 2 )
438: END IF
439: *
440: 50 CONTINUE
441: *
442: 60 CONTINUE
443: *
444: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
445: *
446: * Solve A**T *X + ISGN*X*B = scale*C.
447: *
448: * The (K,L)th block of X is determined starting from
449: * upper-left corner column by column by
450: *
451: * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
452: *
453: * Where
454: * K-1 T L-1
455: * R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
456: * I=1 J=1
457: *
458: * Start column loop (index = L)
459: * L1 (L2): column index of the first (last) row of X(K,L)
460: *
461: LNEXT = 1
462: DO 120 L = 1, N
463: IF( L.LT.LNEXT )
464: $ GO TO 120
465: IF( L.EQ.N ) THEN
466: L1 = L
467: L2 = L
468: ELSE
469: IF( B( L+1, L ).NE.ZERO ) THEN
470: L1 = L
471: L2 = L + 1
472: LNEXT = L + 2
473: ELSE
474: L1 = L
475: L2 = L
476: LNEXT = L + 1
477: END IF
478: END IF
479: *
480: * Start row loop (index = K)
481: * K1 (K2): row index of the first (last) row of X(K,L)
482: *
483: KNEXT = 1
484: DO 110 K = 1, M
485: IF( K.LT.KNEXT )
486: $ GO TO 110
487: IF( K.EQ.M ) THEN
488: K1 = K
489: K2 = K
490: ELSE
491: IF( A( K+1, K ).NE.ZERO ) THEN
492: K1 = K
493: K2 = K + 1
494: KNEXT = K + 2
495: ELSE
496: K1 = K
497: K2 = K
498: KNEXT = K + 1
499: END IF
500: END IF
501: *
502: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
503: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
504: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
505: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
506: SCALOC = ONE
507: *
508: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
509: DA11 = ABS( A11 )
510: IF( DA11.LE.SMIN ) THEN
511: A11 = SMIN
512: DA11 = SMIN
513: INFO = 1
514: END IF
515: DB = ABS( VEC( 1, 1 ) )
516: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
517: IF( DB.GT.BIGNUM*DA11 )
518: $ SCALOC = ONE / DB
519: END IF
520: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
521: *
522: IF( SCALOC.NE.ONE ) THEN
523: DO 70 J = 1, N
524: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
525: 70 CONTINUE
526: SCALE = SCALE*SCALOC
527: END IF
528: C( K1, L1 ) = X( 1, 1 )
529: *
530: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
531: *
532: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
533: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
534: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
535: *
536: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
537: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
538: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
539: *
540: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
541: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
542: $ ZERO, X, 2, SCALOC, XNORM, IERR )
543: IF( IERR.NE.0 )
544: $ INFO = 1
545: *
546: IF( SCALOC.NE.ONE ) THEN
547: DO 80 J = 1, N
548: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
549: 80 CONTINUE
550: SCALE = SCALE*SCALOC
551: END IF
552: C( K1, L1 ) = X( 1, 1 )
553: C( K2, L1 ) = X( 2, 1 )
554: *
555: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
556: *
557: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
558: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
559: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
560: *
561: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
562: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
563: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
564: *
565: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ),
566: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
567: $ ZERO, X, 2, SCALOC, XNORM, IERR )
568: IF( IERR.NE.0 )
569: $ INFO = 1
570: *
571: IF( SCALOC.NE.ONE ) THEN
572: DO 90 J = 1, N
573: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
574: 90 CONTINUE
575: SCALE = SCALE*SCALOC
576: END IF
577: C( K1, L1 ) = X( 1, 1 )
578: C( K1, L2 ) = X( 2, 1 )
579: *
580: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
581: *
582: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
583: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 )
584: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
585: *
586: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
587: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 )
588: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
589: *
590: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
591: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 )
592: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
593: *
594: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
595: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 )
596: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
597: *
598: CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ),
599: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
600: $ 2, XNORM, IERR )
601: IF( IERR.NE.0 )
602: $ INFO = 1
603: *
604: IF( SCALOC.NE.ONE ) THEN
605: DO 100 J = 1, N
606: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
607: 100 CONTINUE
608: SCALE = SCALE*SCALOC
609: END IF
610: C( K1, L1 ) = X( 1, 1 )
611: C( K1, L2 ) = X( 1, 2 )
612: C( K2, L1 ) = X( 2, 1 )
613: C( K2, L2 ) = X( 2, 2 )
614: END IF
615: *
616: 110 CONTINUE
617: 120 CONTINUE
618: *
619: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
620: *
621: * Solve A**T*X + ISGN*X*B**T = scale*C.
622: *
623: * The (K,L)th block of X is determined starting from
624: * top-right corner column by column by
625: *
626: * A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
627: *
628: * Where
629: * K-1 N
630: * R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
631: * I=1 J=L+1
632: *
633: * Start column loop (index = L)
634: * L1 (L2): column index of the first (last) row of X(K,L)
635: *
636: LNEXT = N
637: DO 180 L = N, 1, -1
638: IF( L.GT.LNEXT )
639: $ GO TO 180
640: IF( L.EQ.1 ) THEN
641: L1 = L
642: L2 = L
643: ELSE
644: IF( B( L, L-1 ).NE.ZERO ) THEN
645: L1 = L - 1
646: L2 = L
647: LNEXT = L - 2
648: ELSE
649: L1 = L
650: L2 = L
651: LNEXT = L - 1
652: END IF
653: END IF
654: *
655: * Start row loop (index = K)
656: * K1 (K2): row index of the first (last) row of X(K,L)
657: *
658: KNEXT = 1
659: DO 170 K = 1, M
660: IF( K.LT.KNEXT )
661: $ GO TO 170
662: IF( K.EQ.M ) THEN
663: K1 = K
664: K2 = K
665: ELSE
666: IF( A( K+1, K ).NE.ZERO ) THEN
667: K1 = K
668: K2 = K + 1
669: KNEXT = K + 2
670: ELSE
671: K1 = K
672: K2 = K
673: KNEXT = K + 1
674: END IF
675: END IF
676: *
677: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
678: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
679: SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
680: $ B( L1, MIN( L1+1, N ) ), LDB )
681: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
682: SCALOC = ONE
683: *
684: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
685: DA11 = ABS( A11 )
686: IF( DA11.LE.SMIN ) THEN
687: A11 = SMIN
688: DA11 = SMIN
689: INFO = 1
690: END IF
691: DB = ABS( VEC( 1, 1 ) )
692: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
693: IF( DB.GT.BIGNUM*DA11 )
694: $ SCALOC = ONE / DB
695: END IF
696: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
697: *
698: IF( SCALOC.NE.ONE ) THEN
699: DO 130 J = 1, N
700: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
701: 130 CONTINUE
702: SCALE = SCALE*SCALOC
703: END IF
704: C( K1, L1 ) = X( 1, 1 )
705: *
706: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
707: *
708: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
709: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
710: $ B( L1, MIN( L2+1, N ) ), LDB )
711: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
712: *
713: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
714: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
715: $ B( L1, MIN( L2+1, N ) ), LDB )
716: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
717: *
718: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ),
719: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
720: $ ZERO, X, 2, SCALOC, XNORM, IERR )
721: IF( IERR.NE.0 )
722: $ INFO = 1
723: *
724: IF( SCALOC.NE.ONE ) THEN
725: DO 140 J = 1, N
726: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
727: 140 CONTINUE
728: SCALE = SCALE*SCALOC
729: END IF
730: C( K1, L1 ) = X( 1, 1 )
731: C( K2, L1 ) = X( 2, 1 )
732: *
733: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
734: *
735: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
736: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
737: $ B( L1, MIN( L2+1, N ) ), LDB )
738: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
739: *
740: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
741: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
742: $ B( L2, MIN( L2+1, N ) ), LDB )
743: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
744: *
745: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
746: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
747: $ ZERO, X, 2, SCALOC, XNORM, IERR )
748: IF( IERR.NE.0 )
749: $ INFO = 1
750: *
751: IF( SCALOC.NE.ONE ) THEN
752: DO 150 J = 1, N
753: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
754: 150 CONTINUE
755: SCALE = SCALE*SCALOC
756: END IF
757: C( K1, L1 ) = X( 1, 1 )
758: C( K1, L2 ) = X( 2, 1 )
759: *
760: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
761: *
762: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 )
763: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
764: $ B( L1, MIN( L2+1, N ) ), LDB )
765: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
766: *
767: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 )
768: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
769: $ B( L2, MIN( L2+1, N ) ), LDB )
770: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
771: *
772: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 )
773: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
774: $ B( L1, MIN( L2+1, N ) ), LDB )
775: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
776: *
777: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 )
778: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
779: $ B( L2, MIN( L2+1, N ) ), LDB )
780: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
781: *
782: CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
783: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
784: $ 2, XNORM, IERR )
785: IF( IERR.NE.0 )
786: $ INFO = 1
787: *
788: IF( SCALOC.NE.ONE ) THEN
789: DO 160 J = 1, N
790: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
791: 160 CONTINUE
792: SCALE = SCALE*SCALOC
793: END IF
794: C( K1, L1 ) = X( 1, 1 )
795: C( K1, L2 ) = X( 1, 2 )
796: C( K2, L1 ) = X( 2, 1 )
797: C( K2, L2 ) = X( 2, 2 )
798: END IF
799: *
800: 170 CONTINUE
801: 180 CONTINUE
802: *
803: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
804: *
805: * Solve A*X + ISGN*X*B**T = scale*C.
806: *
807: * The (K,L)th block of X is determined starting from
808: * bottom-right corner column by column by
809: *
810: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
811: *
812: * Where
813: * M N
814: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
815: * I=K+1 J=L+1
816: *
817: * Start column loop (index = L)
818: * L1 (L2): column index of the first (last) row of X(K,L)
819: *
820: LNEXT = N
821: DO 240 L = N, 1, -1
822: IF( L.GT.LNEXT )
823: $ GO TO 240
824: IF( L.EQ.1 ) THEN
825: L1 = L
826: L2 = L
827: ELSE
828: IF( B( L, L-1 ).NE.ZERO ) THEN
829: L1 = L - 1
830: L2 = L
831: LNEXT = L - 2
832: ELSE
833: L1 = L
834: L2 = L
835: LNEXT = L - 1
836: END IF
837: END IF
838: *
839: * Start row loop (index = K)
840: * K1 (K2): row index of the first (last) row of X(K,L)
841: *
842: KNEXT = M
843: DO 230 K = M, 1, -1
844: IF( K.GT.KNEXT )
845: $ GO TO 230
846: IF( K.EQ.1 ) THEN
847: K1 = K
848: K2 = K
849: ELSE
850: IF( A( K, K-1 ).NE.ZERO ) THEN
851: K1 = K - 1
852: K2 = K
853: KNEXT = K - 2
854: ELSE
855: K1 = K
856: K2 = K
857: KNEXT = K - 1
858: END IF
859: END IF
860: *
861: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN
862: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
863: $ C( MIN( K1+1, M ), L1 ), 1 )
864: SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC,
865: $ B( L1, MIN( L1+1, N ) ), LDB )
866: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
867: SCALOC = ONE
868: *
869: A11 = A( K1, K1 ) + SGN*B( L1, L1 )
870: DA11 = ABS( A11 )
871: IF( DA11.LE.SMIN ) THEN
872: A11 = SMIN
873: DA11 = SMIN
874: INFO = 1
875: END IF
876: DB = ABS( VEC( 1, 1 ) )
877: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN
878: IF( DB.GT.BIGNUM*DA11 )
879: $ SCALOC = ONE / DB
880: END IF
881: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11
882: *
883: IF( SCALOC.NE.ONE ) THEN
884: DO 190 J = 1, N
885: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
886: 190 CONTINUE
887: SCALE = SCALE*SCALOC
888: END IF
889: C( K1, L1 ) = X( 1, 1 )
890: *
891: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN
892: *
893: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
894: $ C( MIN( K2+1, M ), L1 ), 1 )
895: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
896: $ B( L1, MIN( L2+1, N ) ), LDB )
897: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
898: *
899: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
900: $ C( MIN( K2+1, M ), L1 ), 1 )
901: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
902: $ B( L1, MIN( L2+1, N ) ), LDB )
903: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
904: *
905: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ),
906: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ),
907: $ ZERO, X, 2, SCALOC, XNORM, IERR )
908: IF( IERR.NE.0 )
909: $ INFO = 1
910: *
911: IF( SCALOC.NE.ONE ) THEN
912: DO 200 J = 1, N
913: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
914: 200 CONTINUE
915: SCALE = SCALE*SCALOC
916: END IF
917: C( K1, L1 ) = X( 1, 1 )
918: C( K2, L1 ) = X( 2, 1 )
919: *
920: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN
921: *
922: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
923: $ C( MIN( K1+1, M ), L1 ), 1 )
924: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
925: $ B( L1, MIN( L2+1, N ) ), LDB )
926: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) )
927: *
928: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA,
929: $ C( MIN( K1+1, M ), L2 ), 1 )
930: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
931: $ B( L2, MIN( L2+1, N ) ), LDB )
932: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) )
933: *
934: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ),
935: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ),
936: $ ZERO, X, 2, SCALOC, XNORM, IERR )
937: IF( IERR.NE.0 )
938: $ INFO = 1
939: *
940: IF( SCALOC.NE.ONE ) THEN
941: DO 210 J = 1, N
942: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
943: 210 CONTINUE
944: SCALE = SCALE*SCALOC
945: END IF
946: C( K1, L1 ) = X( 1, 1 )
947: C( K1, L2 ) = X( 2, 1 )
948: *
949: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN
950: *
951: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
952: $ C( MIN( K2+1, M ), L1 ), 1 )
953: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
954: $ B( L1, MIN( L2+1, N ) ), LDB )
955: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR )
956: *
957: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA,
958: $ C( MIN( K2+1, M ), L2 ), 1 )
959: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC,
960: $ B( L2, MIN( L2+1, N ) ), LDB )
961: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR )
962: *
963: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
964: $ C( MIN( K2+1, M ), L1 ), 1 )
965: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
966: $ B( L1, MIN( L2+1, N ) ), LDB )
967: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR )
968: *
969: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA,
970: $ C( MIN( K2+1, M ), L2 ), 1 )
971: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC,
972: $ B( L2, MIN( L2+1, N ) ), LDB )
973: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR )
974: *
975: CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ),
976: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X,
977: $ 2, XNORM, IERR )
978: IF( IERR.NE.0 )
979: $ INFO = 1
980: *
981: IF( SCALOC.NE.ONE ) THEN
982: DO 220 J = 1, N
983: CALL DSCAL( M, SCALOC, C( 1, J ), 1 )
984: 220 CONTINUE
985: SCALE = SCALE*SCALOC
986: END IF
987: C( K1, L1 ) = X( 1, 1 )
988: C( K1, L2 ) = X( 1, 2 )
989: C( K2, L1 ) = X( 2, 1 )
990: C( K2, L2 ) = X( 2, 2 )
991: END IF
992: *
993: 230 CONTINUE
994: 240 CONTINUE
995: *
996: END IF
997: *
998: RETURN
999: *
1000: * End of DTRSYL
1001: *
1002: 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