LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
strsyl.f
Go to the documentation of this file.
1 *> \brief \b STRSYL
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download STRSYL + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/strsyl.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/strsyl.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/strsyl.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE STRSYL( 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 * REAL SCALE
28 * ..
29 * .. Array Arguments ..
30 * REAL A( LDA, * ), B( LDB, * ), C( LDC, * )
31 * ..
32 *
33 *
34 *> \par Purpose:
35 * =============
36 *>
37 *> \verbatim
38 *>
39 *> STRSYL 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 SHSEQR), 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 REAL 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 REAL 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 REAL 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 REAL
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 realSYcomputational
162 *
163 * =====================================================================
164  SUBROUTINE strsyl( 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  REAL scale
176 * ..
177 * .. Array Arguments ..
178  REAL a( lda, * ), b( ldb, * ), c( ldc, * )
179 * ..
180 *
181 * =====================================================================
182 *
183 * .. Parameters ..
184  REAL zero, one
185  parameter( zero = 0.0e+0, one = 1.0e+0 )
186 * ..
187 * .. Local Scalars ..
188  LOGICAL notrna, notrnb
189  INTEGER ierr, j, k, k1, k2, knext, l, l1, l2, lnext
190  REAL a11, bignum, da11, db, eps, scaloc, sgn, smin,
191  $ smlnum, suml, sumr, xnorm
192 * ..
193 * .. Local Arrays ..
194  REAL dum( 1 ), vec( 2, 2 ), x( 2, 2 )
195 * ..
196 * .. External Functions ..
197  LOGICAL lsame
198  REAL sdot, slamch, slange
199  EXTERNAL lsame, sdot, slamch, slange
200 * ..
201 * .. External Subroutines ..
202  EXTERNAL slabad, slaln2, slasy2, sscal, xerbla
203 * ..
204 * .. Intrinsic Functions ..
205  INTRINSIC abs, max, min, real
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( 'STRSYL', -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 = slamch( 'P' )
248  smlnum = slamch( 'S' )
249  bignum = one / smlnum
250  CALL slabad( smlnum, bignum )
251  smlnum = smlnum*REAL( M*N ) / eps
252  bignum = one / smlnum
253 *
254  smin = max( smlnum, eps*slange( 'M', m, m, a, lda, dum ),
255  $ eps*slange( '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 70 l = 1, n
278  IF( l.LT.lnext )
279  $ go to 70
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 60 k = m, 1, -1
300  IF( k.GT.knext )
301  $ go to 60
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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
319  $ c( min( k1+1, m ), l1 ), 1 )
320  sumr = sdot( 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 sscal( 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
349  $ c( min( k2+1, m ), l1 ), 1 )
350  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
351  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
352 *
353  suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
354  $ c( min( k2+1, m ), l1 ), 1 )
355  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
356  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
357 *
358  CALL slaln2( .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 sscal( 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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
376  $ c( min( k1+1, m ), l1 ), 1 )
377  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
378  vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
379 *
380  suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
381  $ c( min( k1+1, m ), l2 ), 1 )
382  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
383  vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
384 *
385  CALL slaln2( .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 40 j = 1, n
393  CALL sscal( m, scaloc, c( 1, j ), 1 )
394  40 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
403  $ c( min( k2+1, m ), l1 ), 1 )
404  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
405  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
406 *
407  suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
408  $ c( min( k2+1, m ), l2 ), 1 )
409  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
410  vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
411 *
412  suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
413  $ c( min( k2+1, m ), l1 ), 1 )
414  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
415  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
416 *
417  suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
418  $ c( min( k2+1, m ), l2 ), 1 )
419  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
420  vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
421 *
422  CALL slasy2( .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 50 j = 1, n
430  CALL sscal( m, scaloc, c( 1, j ), 1 )
431  50 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  60 continue
441 *
442  70 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 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 130 l = 1, n
463  IF( l.LT.lnext )
464  $ go to 130
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 120 k = 1, m
485  IF( k.LT.knext )
486  $ go to 120
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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
504  sumr = sdot( 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 80 j = 1, n
524  CALL sscal( m, scaloc, c( 1, j ), 1 )
525  80 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
533  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
534  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
535 *
536  suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
537  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
538  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
539 *
540  CALL slaln2( .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 90 j = 1, n
548  CALL sscal( m, scaloc, c( 1, j ), 1 )
549  90 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
558  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
559  vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
560 *
561  suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
562  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
563  vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
564 *
565  CALL slaln2( .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 100 j = 1, n
573  CALL sscal( m, scaloc, c( 1, j ), 1 )
574  100 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
583  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
584  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
585 *
586  suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
587  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
588  vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
589 *
590  suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
591  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
592  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
593 *
594  suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
595  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
596  vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
597 *
598  CALL slasy2( .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 110 j = 1, n
606  CALL sscal( m, scaloc, c( 1, j ), 1 )
607  110 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  120 continue
617  130 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 190 l = n, 1, -1
638  IF( l.GT.lnext )
639  $ go to 190
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 180 k = 1, m
660  IF( k.LT.knext )
661  $ go to 180
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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
679  sumr = sdot( 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 140 j = 1, n
700  CALL sscal( m, scaloc, c( 1, j ), 1 )
701  140 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
709  sumr = sdot( 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 = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
714  sumr = sdot( 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 slaln2( .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 150 j = 1, n
726  CALL sscal( m, scaloc, c( 1, j ), 1 )
727  150 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
736  sumr = sdot( 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
741  sumr = sdot( 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 slaln2( .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 160 j = 1, n
753  CALL sscal( m, scaloc, c( 1, j ), 1 )
754  160 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
763  sumr = sdot( 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
768  sumr = sdot( 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 = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
773  sumr = sdot( 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 = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
778  sumr = sdot( 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 slasy2( .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 170 j = 1, n
790  CALL sscal( m, scaloc, c( 1, j ), 1 )
791  170 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  180 continue
801  190 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 250 l = n, 1, -1
822  IF( l.GT.lnext )
823  $ go to 250
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 240 k = m, 1, -1
844  IF( k.GT.knext )
845  $ go to 240
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 = sdot( m-k1, a( k1, min(k1+1, m ) ), lda,
863  $ c( min( k1+1, m ), l1 ), 1 )
864  sumr = sdot( 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 200 j = 1, n
885  CALL sscal( m, scaloc, c( 1, j ), 1 )
886  200 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
894  $ c( min( k2+1, m ), l1 ), 1 )
895  sumr = sdot( 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 = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
900  $ c( min( k2+1, m ), l1 ), 1 )
901  sumr = sdot( 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 slaln2( .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 210 j = 1, n
913  CALL sscal( m, scaloc, c( 1, j ), 1 )
914  210 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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
923  $ c( min( k1+1, m ), l1 ), 1 )
924  sumr = sdot( 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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
929  $ c( min( k1+1, m ), l2 ), 1 )
930  sumr = sdot( 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 slaln2( .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 220 j = 1, n
942  CALL sscal( m, scaloc, c( 1, j ), 1 )
943  220 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
952  $ c( min( k2+1, m ), l1 ), 1 )
953  sumr = sdot( 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
958  $ c( min( k2+1, m ), l2 ), 1 )
959  sumr = sdot( 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 = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
964  $ c( min( k2+1, m ), l1 ), 1 )
965  sumr = sdot( 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 = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
970  $ c( min( k2+1, m ), l2 ), 1 )
971  sumr = sdot( 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 slasy2( .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 230 j = 1, n
983  CALL sscal( m, scaloc, c( 1, j ), 1 )
984  230 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  240 continue
994  250 continue
995 *
996  END IF
997 *
998  return
999 *
1000 * End of STRSYL
1001 *
1002  END