LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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*> \htmlonly
9*> Download DTRSYL + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrsyl.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrsyl.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrsyl.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
22* LDC, SCALE, INFO )
23*
24* .. Scalar Arguments ..
25* CHARACTER TRANA, TRANB
26* INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
27* DOUBLE PRECISION SCALE
28* ..
29* .. Array Arguments ..
30* DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
31* ..
32*
33*
34*> \par Purpose:
35* =============
36*>
37*> \verbatim
38*>
39*> DTRSYL solves the real Sylvester matrix equation:
40*>
41*> op(A)*X + X*op(B) = scale*C or
42*> op(A)*X - X*op(B) = scale*C,
43*>
44*> where op(A) = A or A**T, and A and B are both upper quasi-
45*> triangular. A is M-by-M and B is N-by-N; the right hand side C and
46*> the solution X are M-by-N; and scale is an output scale factor, set
47*> <= 1 to avoid overflow in X.
48*>
49*> A and B must be in Schur canonical form (as returned by DHSEQR), that
50*> is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
51*> each 2-by-2 diagonal block has its diagonal elements equal and its
52*> off-diagonal elements of opposite sign.
53*> \endverbatim
54*
55* Arguments:
56* ==========
57*
58*> \param[in] TRANA
59*> \verbatim
60*> TRANA is CHARACTER*1
61*> Specifies the option op(A):
62*> = 'N': op(A) = A (No transpose)
63*> = 'T': op(A) = A**T (Transpose)
64*> = 'C': op(A) = A**H (Conjugate transpose = Transpose)
65*> \endverbatim
66*>
67*> \param[in] TRANB
68*> \verbatim
69*> TRANB is CHARACTER*1
70*> Specifies the option op(B):
71*> = 'N': op(B) = B (No transpose)
72*> = 'T': op(B) = B**T (Transpose)
73*> = 'C': op(B) = B**H (Conjugate transpose = Transpose)
74*> \endverbatim
75*>
76*> \param[in] ISGN
77*> \verbatim
78*> ISGN is INTEGER
79*> Specifies the sign in the equation:
80*> = +1: solve op(A)*X + X*op(B) = scale*C
81*> = -1: solve op(A)*X - X*op(B) = scale*C
82*> \endverbatim
83*>
84*> \param[in] M
85*> \verbatim
86*> M is INTEGER
87*> The order of the matrix A, and the number of rows in the
88*> matrices X and C. M >= 0.
89*> \endverbatim
90*>
91*> \param[in] N
92*> \verbatim
93*> N is INTEGER
94*> The order of the matrix B, and the number of columns in the
95*> matrices X and C. N >= 0.
96*> \endverbatim
97*>
98*> \param[in] A
99*> \verbatim
100*> A is DOUBLE PRECISION array, dimension (LDA,M)
101*> The upper quasi-triangular matrix A, in Schur canonical form.
102*> \endverbatim
103*>
104*> \param[in] LDA
105*> \verbatim
106*> LDA is INTEGER
107*> The leading dimension of the array A. LDA >= max(1,M).
108*> \endverbatim
109*>
110*> \param[in] B
111*> \verbatim
112*> B is DOUBLE PRECISION array, dimension (LDB,N)
113*> The upper quasi-triangular matrix B, in Schur canonical form.
114*> \endverbatim
115*>
116*> \param[in] LDB
117*> \verbatim
118*> LDB is INTEGER
119*> The leading dimension of the array B. LDB >= max(1,N).
120*> \endverbatim
121*>
122*> \param[in,out] C
123*> \verbatim
124*> C is DOUBLE PRECISION array, dimension (LDC,N)
125*> On entry, the M-by-N right hand side matrix C.
126*> On exit, C is overwritten by the solution matrix X.
127*> \endverbatim
128*>
129*> \param[in] LDC
130*> \verbatim
131*> LDC is INTEGER
132*> The leading dimension of the array C. LDC >= max(1,M)
133*> \endverbatim
134*>
135*> \param[out] SCALE
136*> \verbatim
137*> SCALE is DOUBLE PRECISION
138*> The scale factor, scale, set <= 1 to avoid overflow in X.
139*> \endverbatim
140*>
141*> \param[out] INFO
142*> \verbatim
143*> INFO is INTEGER
144*> = 0: successful exit
145*> < 0: if INFO = -i, the i-th argument had an illegal value
146*> = 1: A and B have common or very close eigenvalues; perturbed
147*> values were used to solve the equation (but the matrices
148*> A and B are unchanged).
149*> \endverbatim
150*
151* Authors:
152* ========
153*
154*> \author Univ. of Tennessee
155*> \author Univ. of California Berkeley
156*> \author Univ. of Colorado Denver
157*> \author NAG Ltd.
158*
159*> \ingroup trsyl
160*
161* =====================================================================
162 SUBROUTINE dtrsyl( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
163 $ LDC, SCALE, INFO )
164*
165* -- LAPACK computational routine --
166* -- LAPACK is a software package provided by Univ. of Tennessee, --
167* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
168*
169* .. Scalar Arguments ..
170 CHARACTER TRANA, TRANB
171 INTEGER INFO, ISGN, LDA, LDB, LDC, M, N
172 DOUBLE PRECISION SCALE
173* ..
174* .. Array Arguments ..
175 DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * )
176* ..
177*
178* =====================================================================
179*
180* .. Parameters ..
181 DOUBLE PRECISION ZERO, ONE
182 parameter( zero = 0.0d+0, one = 1.0d+0 )
183* ..
184* .. Local Scalars ..
185 LOGICAL NOTRNA, NOTRNB
186 INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
187 DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
188 $ smlnum, suml, sumr, xnorm
189* ..
190* .. Local Arrays ..
191 DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
192* ..
193* .. External Functions ..
194 LOGICAL LSAME
195 DOUBLE PRECISION DDOT, DLAMCH, DLANGE
196 EXTERNAL lsame, ddot, dlamch, dlange
197* ..
198* .. External Subroutines ..
199 EXTERNAL dlaln2, dlasy2, dscal, xerbla
200* ..
201* .. Intrinsic Functions ..
202 INTRINSIC abs, dble, max, min
203* ..
204* .. Executable Statements ..
205*
206* Decode and Test input parameters
207*
208 notrna = lsame( trana, 'N' )
209 notrnb = lsame( tranb, 'N' )
210*
211 info = 0
212 IF( .NOT.notrna .AND. .NOT.lsame( trana, 'T' ) .AND. .NOT.
213 $ lsame( trana, 'C' ) ) THEN
214 info = -1
215 ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'T' ) .AND. .NOT.
216 $ lsame( tranb, 'C' ) ) THEN
217 info = -2
218 ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
219 info = -3
220 ELSE IF( m.LT.0 ) THEN
221 info = -4
222 ELSE IF( n.LT.0 ) THEN
223 info = -5
224 ELSE IF( lda.LT.max( 1, m ) ) THEN
225 info = -7
226 ELSE IF( ldb.LT.max( 1, n ) ) THEN
227 info = -9
228 ELSE IF( ldc.LT.max( 1, m ) ) THEN
229 info = -11
230 END IF
231 IF( info.NE.0 ) THEN
232 CALL xerbla( 'DTRSYL', -info )
233 RETURN
234 END IF
235*
236* Quick return if possible
237*
238 scale = one
239 IF( m.EQ.0 .OR. n.EQ.0 )
240 $ RETURN
241*
242* Set constants to control overflow
243*
244 eps = dlamch( 'P' )
245 smlnum = dlamch( 'S' )
246 bignum = one / smlnum
247 smlnum = smlnum*dble( m*n ) / eps
248 bignum = one / smlnum
249*
250 smin = max( smlnum, eps*dlange( 'M', m, m, a, lda, dum ),
251 $ eps*dlange( 'M', n, n, b, ldb, dum ) )
252*
253 sgn = isgn
254*
255 IF( notrna .AND. notrnb ) THEN
256*
257* Solve A*X + ISGN*X*B = scale*C.
258*
259* The (K,L)th block of X is determined starting from
260* bottom-left corner column by column by
261*
262* A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
263*
264* Where
265* M L-1
266* R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
267* I=K+1 J=1
268*
269* Start column loop (index = L)
270* L1 (L2) : column index of the first (first) row of X(K,L).
271*
272 lnext = 1
273 DO 60 l = 1, n
274 IF( l.LT.lnext )
275 $ GO TO 60
276 IF( l.EQ.n ) THEN
277 l1 = l
278 l2 = l
279 ELSE
280 IF( b( l+1, l ).NE.zero ) THEN
281 l1 = l
282 l2 = l + 1
283 lnext = l + 2
284 ELSE
285 l1 = l
286 l2 = l
287 lnext = l + 1
288 END IF
289 END IF
290*
291* Start row loop (index = K)
292* K1 (K2): row index of the first (last) row of X(K,L).
293*
294 knext = m
295 DO 50 k = m, 1, -1
296 IF( k.GT.knext )
297 $ GO TO 50
298 IF( k.EQ.1 ) THEN
299 k1 = k
300 k2 = k
301 ELSE
302 IF( a( k, k-1 ).NE.zero ) THEN
303 k1 = k - 1
304 k2 = k
305 knext = k - 2
306 ELSE
307 k1 = k
308 k2 = k
309 knext = k - 1
310 END IF
311 END IF
312*
313 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
314 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
315 $ c( min( k1+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 scaloc = one
319*
320 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
321 da11 = abs( a11 )
322 IF( da11.LE.smin ) THEN
323 a11 = smin
324 da11 = smin
325 info = 1
326 END IF
327 db = abs( vec( 1, 1 ) )
328 IF( da11.LT.one .AND. db.GT.one ) THEN
329 IF( db.GT.bignum*da11 )
330 $ scaloc = one / db
331 END IF
332 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
333*
334 IF( scaloc.NE.one ) THEN
335 DO 10 j = 1, n
336 CALL dscal( m, scaloc, c( 1, j ), 1 )
337 10 CONTINUE
338 scale = scale*scaloc
339 END IF
340 c( k1, l1 ) = x( 1, 1 )
341*
342 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
343*
344 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
345 $ c( min( k2+1, m ), l1 ), 1 )
346 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
347 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
348*
349 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
350 $ c( min( k2+1, m ), l1 ), 1 )
351 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
352 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
353*
354 CALL dlaln2( .false., 2, 1, smin, one, a( k1, k1 ),
355 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
356 $ zero, x, 2, scaloc, xnorm, ierr )
357 IF( ierr.NE.0 )
358 $ info = 1
359*
360 IF( scaloc.NE.one ) THEN
361 DO 20 j = 1, n
362 CALL dscal( m, scaloc, c( 1, j ), 1 )
363 20 CONTINUE
364 scale = scale*scaloc
365 END IF
366 c( k1, l1 ) = x( 1, 1 )
367 c( k2, l1 ) = x( 2, 1 )
368*
369 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
370*
371 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
372 $ c( min( k1+1, m ), l1 ), 1 )
373 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
374 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
375*
376 suml = ddot( m-k1, a( k1, min( k1+1, m ) ), lda,
377 $ c( min( k1+1, m ), l2 ), 1 )
378 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
379 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
380*
381 CALL dlaln2( .true., 2, 1, smin, one, b( l1, l1 ),
382 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
383 $ zero, x, 2, scaloc, xnorm, ierr )
384 IF( ierr.NE.0 )
385 $ info = 1
386*
387 IF( scaloc.NE.one ) THEN
388 DO 30 j = 1, n
389 CALL dscal( m, scaloc, c( 1, j ), 1 )
390 30 CONTINUE
391 scale = scale*scaloc
392 END IF
393 c( k1, l1 ) = x( 1, 1 )
394 c( k1, l2 ) = x( 2, 1 )
395*
396 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
397*
398 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
399 $ c( min( k2+1, m ), l1 ), 1 )
400 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
401 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
402*
403 suml = ddot( m-k2, a( k1, min( k2+1, m ) ), lda,
404 $ c( min( k2+1, m ), l2 ), 1 )
405 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
406 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
407*
408 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
409 $ c( min( k2+1, m ), l1 ), 1 )
410 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
411 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
412*
413 suml = ddot( m-k2, a( k2, min( k2+1, m ) ), lda,
414 $ c( min( k2+1, m ), l2 ), 1 )
415 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
416 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
417*
418 CALL dlasy2( .false., .false., isgn, 2, 2,
419 $ a( k1, k1 ), lda, b( l1, l1 ), ldb, vec,
420 $ 2, scaloc, x, 2, xnorm, ierr )
421 IF( ierr.NE.0 )
422 $ info = 1
423*
424 IF( scaloc.NE.one ) THEN
425 DO 40 j = 1, n
426 CALL dscal( m, scaloc, c( 1, j ), 1 )
427 40 CONTINUE
428 scale = scale*scaloc
429 END IF
430 c( k1, l1 ) = x( 1, 1 )
431 c( k1, l2 ) = x( 1, 2 )
432 c( k2, l1 ) = x( 2, 1 )
433 c( k2, l2 ) = x( 2, 2 )
434 END IF
435*
436 50 CONTINUE
437*
438 60 CONTINUE
439*
440 ELSE IF( .NOT.notrna .AND. notrnb ) THEN
441*
442* Solve A**T *X + ISGN*X*B = scale*C.
443*
444* The (K,L)th block of X is determined starting from
445* upper-left corner column by column by
446*
447* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
448*
449* Where
450* K-1 T L-1
451* R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
452* I=1 J=1
453*
454* Start column loop (index = L)
455* L1 (L2): column index of the first (last) row of X(K,L)
456*
457 lnext = 1
458 DO 120 l = 1, n
459 IF( l.LT.lnext )
460 $ GO TO 120
461 IF( l.EQ.n ) THEN
462 l1 = l
463 l2 = l
464 ELSE
465 IF( b( l+1, l ).NE.zero ) THEN
466 l1 = l
467 l2 = l + 1
468 lnext = l + 2
469 ELSE
470 l1 = l
471 l2 = l
472 lnext = l + 1
473 END IF
474 END IF
475*
476* Start row loop (index = K)
477* K1 (K2): row index of the first (last) row of X(K,L)
478*
479 knext = 1
480 DO 110 k = 1, m
481 IF( k.LT.knext )
482 $ GO TO 110
483 IF( k.EQ.m ) THEN
484 k1 = k
485 k2 = k
486 ELSE
487 IF( a( k+1, k ).NE.zero ) THEN
488 k1 = k
489 k2 = k + 1
490 knext = k + 2
491 ELSE
492 k1 = k
493 k2 = k
494 knext = k + 1
495 END IF
496 END IF
497*
498 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
499 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
500 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
501 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
502 scaloc = one
503*
504 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
505 da11 = abs( a11 )
506 IF( da11.LE.smin ) THEN
507 a11 = smin
508 da11 = smin
509 info = 1
510 END IF
511 db = abs( vec( 1, 1 ) )
512 IF( da11.LT.one .AND. db.GT.one ) THEN
513 IF( db.GT.bignum*da11 )
514 $ scaloc = one / db
515 END IF
516 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
517*
518 IF( scaloc.NE.one ) THEN
519 DO 70 j = 1, n
520 CALL dscal( m, scaloc, c( 1, j ), 1 )
521 70 CONTINUE
522 scale = scale*scaloc
523 END IF
524 c( k1, l1 ) = x( 1, 1 )
525*
526 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
527*
528 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
529 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
530 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
531*
532 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
533 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
534 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
535*
536 CALL dlaln2( .true., 2, 1, smin, one, a( k1, k1 ),
537 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
538 $ zero, x, 2, scaloc, xnorm, ierr )
539 IF( ierr.NE.0 )
540 $ info = 1
541*
542 IF( scaloc.NE.one ) THEN
543 DO 80 j = 1, n
544 CALL dscal( m, scaloc, c( 1, j ), 1 )
545 80 CONTINUE
546 scale = scale*scaloc
547 END IF
548 c( k1, l1 ) = x( 1, 1 )
549 c( k2, l1 ) = x( 2, 1 )
550*
551 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
552*
553 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
554 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
555 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
556*
557 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
558 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
559 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
560*
561 CALL dlaln2( .true., 2, 1, smin, one, b( l1, l1 ),
562 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
563 $ zero, x, 2, scaloc, xnorm, ierr )
564 IF( ierr.NE.0 )
565 $ info = 1
566*
567 IF( scaloc.NE.one ) THEN
568 DO 90 j = 1, n
569 CALL dscal( m, scaloc, c( 1, j ), 1 )
570 90 CONTINUE
571 scale = scale*scaloc
572 END IF
573 c( k1, l1 ) = x( 1, 1 )
574 c( k1, l2 ) = x( 2, 1 )
575*
576 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
577*
578 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
579 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
580 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
581*
582 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
583 sumr = ddot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
584 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
585*
586 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
587 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
588 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
589*
590 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
591 sumr = ddot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
592 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
593*
594 CALL dlasy2( .true., .false., isgn, 2, 2, a( k1, k1 ),
595 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
596 $ 2, xnorm, ierr )
597 IF( ierr.NE.0 )
598 $ info = 1
599*
600 IF( scaloc.NE.one ) THEN
601 DO 100 j = 1, n
602 CALL dscal( m, scaloc, c( 1, j ), 1 )
603 100 CONTINUE
604 scale = scale*scaloc
605 END IF
606 c( k1, l1 ) = x( 1, 1 )
607 c( k1, l2 ) = x( 1, 2 )
608 c( k2, l1 ) = x( 2, 1 )
609 c( k2, l2 ) = x( 2, 2 )
610 END IF
611*
612 110 CONTINUE
613 120 CONTINUE
614*
615 ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
616*
617* Solve A**T*X + ISGN*X*B**T = scale*C.
618*
619* The (K,L)th block of X is determined starting from
620* top-right corner column by column by
621*
622* A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
623*
624* Where
625* K-1 N
626* R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
627* I=1 J=L+1
628*
629* Start column loop (index = L)
630* L1 (L2): column index of the first (last) row of X(K,L)
631*
632 lnext = n
633 DO 180 l = n, 1, -1
634 IF( l.GT.lnext )
635 $ GO TO 180
636 IF( l.EQ.1 ) THEN
637 l1 = l
638 l2 = l
639 ELSE
640 IF( b( l, l-1 ).NE.zero ) THEN
641 l1 = l - 1
642 l2 = l
643 lnext = l - 2
644 ELSE
645 l1 = l
646 l2 = l
647 lnext = l - 1
648 END IF
649 END IF
650*
651* Start row loop (index = K)
652* K1 (K2): row index of the first (last) row of X(K,L)
653*
654 knext = 1
655 DO 170 k = 1, m
656 IF( k.LT.knext )
657 $ GO TO 170
658 IF( k.EQ.m ) THEN
659 k1 = k
660 k2 = k
661 ELSE
662 IF( a( k+1, k ).NE.zero ) THEN
663 k1 = k
664 k2 = k + 1
665 knext = k + 2
666 ELSE
667 k1 = k
668 k2 = k
669 knext = k + 1
670 END IF
671 END IF
672*
673 IF( l1.EQ.l2 .AND. k1.EQ.k2 ) THEN
674 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
675 sumr = ddot( n-l1, c( k1, min( l1+1, n ) ), ldc,
676 $ b( l1, min( l1+1, n ) ), ldb )
677 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
678 scaloc = one
679*
680 a11 = a( k1, k1 ) + sgn*b( l1, l1 )
681 da11 = abs( a11 )
682 IF( da11.LE.smin ) THEN
683 a11 = smin
684 da11 = smin
685 info = 1
686 END IF
687 db = abs( vec( 1, 1 ) )
688 IF( da11.LT.one .AND. db.GT.one ) THEN
689 IF( db.GT.bignum*da11 )
690 $ scaloc = one / db
691 END IF
692 x( 1, 1 ) = ( vec( 1, 1 )*scaloc ) / a11
693*
694 IF( scaloc.NE.one ) THEN
695 DO 130 j = 1, n
696 CALL dscal( m, scaloc, c( 1, j ), 1 )
697 130 CONTINUE
698 scale = scale*scaloc
699 END IF
700 c( k1, l1 ) = x( 1, 1 )
701*
702 ELSE IF( l1.EQ.l2 .AND. k1.NE.k2 ) THEN
703*
704 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
705 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
706 $ b( l1, min( l2+1, n ) ), ldb )
707 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
708*
709 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
710 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
711 $ b( l1, min( l2+1, n ) ), ldb )
712 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
713*
714 CALL dlaln2( .true., 2, 1, smin, one, a( k1, k1 ),
715 $ lda, one, one, vec, 2, -sgn*b( l1, l1 ),
716 $ zero, x, 2, scaloc, xnorm, ierr )
717 IF( ierr.NE.0 )
718 $ info = 1
719*
720 IF( scaloc.NE.one ) THEN
721 DO 140 j = 1, n
722 CALL dscal( m, scaloc, c( 1, j ), 1 )
723 140 CONTINUE
724 scale = scale*scaloc
725 END IF
726 c( k1, l1 ) = x( 1, 1 )
727 c( k2, l1 ) = x( 2, 1 )
728*
729 ELSE IF( l1.NE.l2 .AND. k1.EQ.k2 ) THEN
730*
731 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
732 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
733 $ b( l1, min( l2+1, n ) ), ldb )
734 vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
735*
736 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
737 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
738 $ b( l2, min( l2+1, n ) ), ldb )
739 vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
740*
741 CALL dlaln2( .false., 2, 1, smin, one, b( l1, l1 ),
742 $ ldb, one, one, vec, 2, -sgn*a( k1, k1 ),
743 $ zero, x, 2, scaloc, xnorm, ierr )
744 IF( ierr.NE.0 )
745 $ info = 1
746*
747 IF( scaloc.NE.one ) THEN
748 DO 150 j = 1, n
749 CALL dscal( m, scaloc, c( 1, j ), 1 )
750 150 CONTINUE
751 scale = scale*scaloc
752 END IF
753 c( k1, l1 ) = x( 1, 1 )
754 c( k1, l2 ) = x( 2, 1 )
755*
756 ELSE IF( l1.NE.l2 .AND. k1.NE.k2 ) THEN
757*
758 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
759 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
760 $ b( l1, min( l2+1, n ) ), ldb )
761 vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
762*
763 suml = ddot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
764 sumr = ddot( n-l2, c( k1, min( l2+1, n ) ), ldc,
765 $ b( l2, min( l2+1, n ) ), ldb )
766 vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
767*
768 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
769 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
770 $ b( l1, min( l2+1, n ) ), ldb )
771 vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
772*
773 suml = ddot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
774 sumr = ddot( n-l2, c( k2, min( l2+1, n ) ), ldc,
775 $ b( l2, min( l2+1, n ) ), ldb )
776 vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
777*
778 CALL dlasy2( .true., .true., isgn, 2, 2, a( k1, 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, k1 ),
972 $ lda, b( l1, l1 ), ldb, vec, 2, scaloc, x,
973 $ 2, xnorm, ierr )
974 IF( ierr.NE.0 )
975 $ info = 1
976*
977 IF( scaloc.NE.one ) THEN
978 DO 220 j = 1, n
979 CALL dscal( m, scaloc, c( 1, j ), 1 )
980 220 CONTINUE
981 scale = scale*scaloc
982 END IF
983 c( k1, l1 ) = x( 1, 1 )
984 c( k1, l2 ) = x( 1, 2 )
985 c( k2, l1 ) = x( 2, 1 )
986 c( k2, l2 ) = x( 2, 2 )
987 END IF
988*
989 230 CONTINUE
990 240 CONTINUE
991*
992 END IF
993*
994 RETURN
995*
996* End of DTRSYL
997*
998 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:218
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:174
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:164