LAPACK  3.10.0
LAPACK: Linear Algebra PACKage

◆ strsyl()

subroutine strsyl ( character  TRANA,
character  TRANB,
integer  ISGN,
integer  M,
integer  N,
real, dimension( lda, * )  A,
integer  LDA,
real, dimension( ldb, * )  B,
integer  LDB,
real, dimension( ldc, * )  C,
integer  LDC,
real  SCALE,
integer  INFO 
)

STRSYL

Download STRSYL + dependencies [TGZ] [ZIP] [TXT]

Purpose:
 STRSYL solves the real Sylvester matrix equation:

    op(A)*X + X*op(B) = scale*C or
    op(A)*X - X*op(B) = scale*C,

 where op(A) = A or A**T, and  A and B are both upper quasi-
 triangular. A is M-by-M and B is N-by-N; the right hand side C and
 the solution X are M-by-N; and scale is an output scale factor, set
 <= 1 to avoid overflow in X.

 A and B must be in Schur canonical form (as returned by SHSEQR), that
 is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks;
 each 2-by-2 diagonal block has its diagonal elements equal and its
 off-diagonal elements of opposite sign.
Parameters
[in]TRANA
          TRANA is CHARACTER*1
          Specifies the option op(A):
          = 'N': op(A) = A    (No transpose)
          = 'T': op(A) = A**T (Transpose)
          = 'C': op(A) = A**H (Conjugate transpose = Transpose)
[in]TRANB
          TRANB is CHARACTER*1
          Specifies the option op(B):
          = 'N': op(B) = B    (No transpose)
          = 'T': op(B) = B**T (Transpose)
          = 'C': op(B) = B**H (Conjugate transpose = Transpose)
[in]ISGN
          ISGN is INTEGER
          Specifies the sign in the equation:
          = +1: solve op(A)*X + X*op(B) = scale*C
          = -1: solve op(A)*X - X*op(B) = scale*C
[in]M
          M is INTEGER
          The order of the matrix A, and the number of rows in the
          matrices X and C. M >= 0.
[in]N
          N is INTEGER
          The order of the matrix B, and the number of columns in the
          matrices X and C. N >= 0.
[in]A
          A is REAL array, dimension (LDA,M)
          The upper quasi-triangular matrix A, in Schur canonical form.
[in]LDA
          LDA is INTEGER
          The leading dimension of the array A. LDA >= max(1,M).
[in]B
          B is REAL array, dimension (LDB,N)
          The upper quasi-triangular matrix B, in Schur canonical form.
[in]LDB
          LDB is INTEGER
          The leading dimension of the array B. LDB >= max(1,N).
[in,out]C
          C is REAL array, dimension (LDC,N)
          On entry, the M-by-N right hand side matrix C.
          On exit, C is overwritten by the solution matrix X.
[in]LDC
          LDC is INTEGER
          The leading dimension of the array C. LDC >= max(1,M)
[out]SCALE
          SCALE is REAL
          The scale factor, scale, set <= 1 to avoid overflow in X.
[out]INFO
          INFO is INTEGER
          = 0: successful exit
          < 0: if INFO = -i, the i-th argument had an illegal value
          = 1: A and B have common or very close eigenvalues; perturbed
               values were used to solve the equation (but the matrices
               A and B are unchanged).
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.

Definition at line 162 of file strsyl.f.

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  REAL SCALE
173 * ..
174 * .. Array Arguments ..
175  REAL A( LDA, * ), B( LDB, * ), C( LDC, * )
176 * ..
177 *
178 * =====================================================================
179 *
180 * .. Parameters ..
181  REAL ZERO, ONE
182  parameter( zero = 0.0e+0, one = 1.0e+0 )
183 * ..
184 * .. Local Scalars ..
185  LOGICAL NOTRNA, NOTRNB
186  INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT
187  REAL A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN,
188  $ SMLNUM, SUML, SUMR, XNORM
189 * ..
190 * .. Local Arrays ..
191  REAL DUM( 1 ), VEC( 2, 2 ), X( 2, 2 )
192 * ..
193 * .. External Functions ..
194  LOGICAL LSAME
195  REAL SDOT, SLAMCH, SLANGE
196  EXTERNAL lsame, sdot, slamch, slange
197 * ..
198 * .. External Subroutines ..
199  EXTERNAL slabad, slaln2, slasy2, sscal, xerbla
200 * ..
201 * .. Intrinsic Functions ..
202  INTRINSIC abs, max, min, real
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( 'STRSYL', -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 = slamch( 'P' )
245  smlnum = slamch( 'S' )
246  bignum = one / smlnum
247  CALL slabad( smlnum, bignum )
248  smlnum = smlnum*real( m*n ) / eps
249  bignum = one / smlnum
250 *
251  smin = max( smlnum, eps*slange( 'M', m, m, a, lda, dum ),
252  $ eps*slange( '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 70 l = 1, n
275  IF( l.LT.lnext )
276  $ GO TO 70
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 60 k = m, 1, -1
297  IF( k.GT.knext )
298  $ GO TO 60
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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
316  $ c( min( k1+1, m ), l1 ), 1 )
317  sumr = sdot( 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 sscal( 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
346  $ c( min( k2+1, m ), l1 ), 1 )
347  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
348  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
349 *
350  suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
351  $ c( min( k2+1, m ), l1 ), 1 )
352  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
353  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
354 *
355  CALL slaln2( .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 sscal( 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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
373  $ c( min( k1+1, m ), l1 ), 1 )
374  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
375  vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
376 *
377  suml = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
378  $ c( min( k1+1, m ), l2 ), 1 )
379  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
380  vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
381 *
382  CALL slaln2( .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 40 j = 1, n
390  CALL sscal( m, scaloc, c( 1, j ), 1 )
391  40 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
400  $ c( min( k2+1, m ), l1 ), 1 )
401  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
402  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
403 *
404  suml = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
405  $ c( min( k2+1, m ), l2 ), 1 )
406  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
407  vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
408 *
409  suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
410  $ c( min( k2+1, m ), l1 ), 1 )
411  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
412  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
413 *
414  suml = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
415  $ c( min( k2+1, m ), l2 ), 1 )
416  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
417  vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
418 *
419  CALL slasy2( .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 50 j = 1, n
427  CALL sscal( m, scaloc, c( 1, j ), 1 )
428  50 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  60 CONTINUE
438 *
439  70 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 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 130 l = 1, n
460  IF( l.LT.lnext )
461  $ GO TO 130
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 120 k = 1, m
482  IF( k.LT.knext )
483  $ GO TO 120
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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
501  sumr = sdot( 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 80 j = 1, n
521  CALL sscal( m, scaloc, c( 1, j ), 1 )
522  80 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
530  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
531  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
532 *
533  suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
534  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
535  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
536 *
537  CALL slaln2( .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 90 j = 1, n
545  CALL sscal( m, scaloc, c( 1, j ), 1 )
546  90 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
555  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
556  vec( 1, 1 ) = sgn*( c( k1, l1 )-( suml+sgn*sumr ) )
557 *
558  suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
559  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
560  vec( 2, 1 ) = sgn*( c( k1, l2 )-( suml+sgn*sumr ) )
561 *
562  CALL slaln2( .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 100 j = 1, n
570  CALL sscal( m, scaloc, c( 1, j ), 1 )
571  100 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
580  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l1 ), 1 )
581  vec( 1, 1 ) = c( k1, l1 ) - ( suml+sgn*sumr )
582 *
583  suml = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
584  sumr = sdot( l1-1, c( k1, 1 ), ldc, b( 1, l2 ), 1 )
585  vec( 1, 2 ) = c( k1, l2 ) - ( suml+sgn*sumr )
586 *
587  suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
588  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l1 ), 1 )
589  vec( 2, 1 ) = c( k2, l1 ) - ( suml+sgn*sumr )
590 *
591  suml = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
592  sumr = sdot( l1-1, c( k2, 1 ), ldc, b( 1, l2 ), 1 )
593  vec( 2, 2 ) = c( k2, l2 ) - ( suml+sgn*sumr )
594 *
595  CALL slasy2( .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 110 j = 1, n
603  CALL sscal( m, scaloc, c( 1, j ), 1 )
604  110 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  120 CONTINUE
614  130 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 190 l = n, 1, -1
635  IF( l.GT.lnext )
636  $ GO TO 190
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 180 k = 1, m
657  IF( k.LT.knext )
658  $ GO TO 180
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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
676  sumr = sdot( 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 140 j = 1, n
697  CALL sscal( m, scaloc, c( 1, j ), 1 )
698  140 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
706  sumr = sdot( 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 = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
711  sumr = sdot( 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 slaln2( .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 150 j = 1, n
723  CALL sscal( m, scaloc, c( 1, j ), 1 )
724  150 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
733  sumr = sdot( 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
738  sumr = sdot( 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 slaln2( .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 160 j = 1, n
750  CALL sscal( m, scaloc, c( 1, j ), 1 )
751  160 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l1 ), 1 )
760  sumr = sdot( 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 = sdot( k1-1, a( 1, k1 ), 1, c( 1, l2 ), 1 )
765  sumr = sdot( 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 = sdot( k1-1, a( 1, k2 ), 1, c( 1, l1 ), 1 )
770  sumr = sdot( 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 = sdot( k1-1, a( 1, k2 ), 1, c( 1, l2 ), 1 )
775  sumr = sdot( 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 slasy2( .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 170 j = 1, n
787  CALL sscal( m, scaloc, c( 1, j ), 1 )
788  170 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  180 CONTINUE
798  190 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 250 l = n, 1, -1
819  IF( l.GT.lnext )
820  $ GO TO 250
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 240 k = m, 1, -1
841  IF( k.GT.knext )
842  $ GO TO 240
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 = sdot( m-k1, a( k1, min(k1+1, m ) ), lda,
860  $ c( min( k1+1, m ), l1 ), 1 )
861  sumr = sdot( 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 200 j = 1, n
882  CALL sscal( m, scaloc, c( 1, j ), 1 )
883  200 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
891  $ c( min( k2+1, m ), l1 ), 1 )
892  sumr = sdot( 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 = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
897  $ c( min( k2+1, m ), l1 ), 1 )
898  sumr = sdot( 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 slaln2( .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 210 j = 1, n
910  CALL sscal( m, scaloc, c( 1, j ), 1 )
911  210 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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
920  $ c( min( k1+1, m ), l1 ), 1 )
921  sumr = sdot( 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 = sdot( m-k1, a( k1, min( k1+1, m ) ), lda,
926  $ c( min( k1+1, m ), l2 ), 1 )
927  sumr = sdot( 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 slaln2( .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 220 j = 1, n
939  CALL sscal( m, scaloc, c( 1, j ), 1 )
940  220 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
949  $ c( min( k2+1, m ), l1 ), 1 )
950  sumr = sdot( 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 = sdot( m-k2, a( k1, min( k2+1, m ) ), lda,
955  $ c( min( k2+1, m ), l2 ), 1 )
956  sumr = sdot( 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 = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
961  $ c( min( k2+1, m ), l1 ), 1 )
962  sumr = sdot( 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 = sdot( m-k2, a( k2, min( k2+1, m ) ), lda,
967  $ c( min( k2+1, m ), l2 ), 1 )
968  sumr = sdot( 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 slasy2( .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 230 j = 1, n
980  CALL sscal( m, scaloc, c( 1, j ), 1 )
981  230 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  240 CONTINUE
991  250 CONTINUE
992 *
993  END IF
994 *
995  RETURN
996 *
997 * End of STRSYL
998 *
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:74
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: slange.f:114
subroutine slaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
SLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition: slaln2.f:218
subroutine slasy2(LTRANL, LTRANR, ISGN, N1, N2, TL, LDTL, TR, LDTR, B, LDB, SCALE, X, LDX, XNORM, INFO)
SLASY2 solves the Sylvester matrix equation where the matrices are of order 1 or 2.
Definition: slasy2.f:174
subroutine sscal(N, SA, SX, INCX)
SSCAL
Definition: sscal.f:79
real function sdot(N, SX, INCX, SY, INCY)
SDOT
Definition: sdot.f:82
real function slamch(CMACH)
SLAMCH
Definition: slamch.f:68
Here is the call graph for this function:
Here is the caller graph for this function: