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