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