1 SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
2 $ LDC, SCALE, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**T *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)**T*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)**T*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**T*X + ISGN*X*B**T = 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)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
539 *
540 * Where
541 * K-1 N
542 * R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
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**T = 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)**T = 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)**T].
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
2 $ LDC, SCALE, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
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**T *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)**T*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)**T*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**T*X + ISGN*X*B**T = 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)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
539 *
540 * Where
541 * K-1 N
542 * R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
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**T = 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)**T = 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)**T].
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