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