dtrsyl3()

subroutine dtrsyl3	(	character	trana,
		character	tranb,
		integer	isgn,
		integer	m,
		integer	n,
		double precision, dimension( lda, * )	a,
		integer	lda,
		double precision, dimension( ldb, * )	b,
		integer	ldb,
		double precision, dimension( ldc, * )	c,
		integer	ldc,
		double precision	scale,
		integer, dimension( * )	iwork,
		integer	liwork,
		double precision, dimension( ldswork, * )	swork,
		integer	ldswork,
		integer	info )

DTRSYL3

Purpose:

!>
!>  DTRSYL3 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 DHSEQR), 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.
!>
!>  This is the block version of the algorithm.
!> 

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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION !> The scale factor, scale, set <= 1 to avoid overflow in X. !>
[out]	IWORK	!> IWORK is INTEGER array, dimension (MAX(1,LIWORK)) !> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK. !>
[in]	LIWORK	!> IWORK is INTEGER !> The dimension of the array IWORK. LIWORK >= ((M + NB - 1) / NB + 1) !> + ((N + NB - 1) / NB + 1), where NB is the optimal block size. !> !> If LIWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal dimension of the IWORK array, !> returns this value as the first entry of the IWORK array, and !> no error message related to LIWORK is issued by XERBLA. !>
[out]	SWORK	!> SWORK is DOUBLE PRECISION array, dimension (MAX(2, ROWS), !> MAX(1,COLS)). !> On exit, if INFO = 0, SWORK(1) returns the optimal value ROWS !> and SWORK(2) returns the optimal COLS. !>
[in]	LDSWORK	!> LDSWORK is INTEGER !> LDSWORK >= MAX(2,ROWS), where ROWS = ((M + NB - 1) / NB + 1) !> and NB is the optimal block size. !> !> If LDSWORK = -1, then a workspace query is assumed; the routine !> only calculates the optimal dimensions of the SWORK matrix, !> returns these values as the first and second entry of the SWORK !> matrix, and no error message related LWORK is issued by XERBLA. !>
[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). !>

Definition at line 197 of file dtrsyl3.f.

      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          TRANA, TRANB
      INTEGER            INFO, ISGN, LDA, LDB, LDC, M, N,
     $                   LIWORK, LDSWORK
      DOUBLE PRECISION   SCALE
*     ..
*     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   A( LDA, * ), B( LDB, * ), C( LDC, * ),
     $                   SWORK( LDSWORK, * )
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d+0, one = 1.0d+0 )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOTRNA, NOTRNB, LQUERY, SKIP
      INTEGER            AWRK, BWRK, I, I1, I2, IINFO, J, J1, J2, JJ,
     $                   K, K1, K2, L, L1, L2, LL, NBA, NB, NBB, PC
      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, CNRM, SCAL, SCALOC,
     $                   SCAMIN, SGN, XNRM, BUF, SMLNUM
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   WNRM( MAX( M, N ) )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      DOUBLE PRECISION   DLANGE, DLAMCH, DLARMM
      EXTERNAL           dlange, dlamch, dlarmm, ilaenv,
     $                   lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           dgemm, dlascl, dscal, dtrsyl,
     $                   xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, exponent, max, min
*     ..
*     .. Executable Statements ..
*
*     Decode and Test input parameters
*
      notrna = lsame( trana, 'N' )
      notrnb = lsame( tranb, 'N' )
*
*     Use the same block size for all matrices.
*
      nb = max(8, ilaenv( 1, 'DTRSYL', '', m, n, -1, -1) )
*
*     Compute number of blocks in A and B
*
      nba = max( 1, (m + nb - 1) / nb )
      nbb = max( 1, (n + nb - 1) / nb )
*
*     Compute workspace
*
      info = 0
      lquery = ( liwork.EQ.-1 .OR. ldswork.EQ.-1 )
      iwork( 1 ) = nba + nbb + 2
      IF( lquery ) THEN
         ldswork = 2
         swork( 1, 1 ) = max( nba, nbb )
         swork( 2, 1 ) = 2 * nbb + nba
      END IF
*
*     Test the input arguments
*
      IF( .NOT.notrna .AND. .NOT.lsame( trana, 'T' ) .AND. .NOT.
     $    lsame( trana, 'C' ) ) THEN
         info = -1
      ELSE IF( .NOT.notrnb .AND. .NOT.lsame( tranb, 'T' ) .AND. .NOT.
     $         lsame( tranb, 'C' ) ) THEN
         info = -2
      ELSE IF( isgn.NE.1 .AND. isgn.NE.-1 ) THEN
         info = -3
      ELSE IF( m.LT.0 ) THEN
         info = -4
      ELSE IF( n.LT.0 ) THEN
         info = -5
      ELSE IF( lda.LT.max( 1, m ) ) THEN
         info = -7
      ELSE IF( ldb.LT.max( 1, n ) ) THEN
         info = -9
      ELSE IF( ldc.LT.max( 1, m ) ) THEN
         info = -11
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'DTRSYL3', -info )
         RETURN
      ELSE IF( lquery ) THEN
         RETURN
      END IF
*
*     Quick return if possible
*
      scale = one
      IF( m.EQ.0 .OR. n.EQ.0 )
     $   RETURN
*
*     Use unblocked code for small problems or if insufficient
*     workspaces are provided
*
      IF( min( nba, nbb ).EQ.1 .OR. ldswork.LT.max( nba, nbb ) .OR.
     $    liwork.LT.iwork(1) ) THEN
        CALL dtrsyl( trana, tranb, isgn, m, n, a, lda, b, ldb,
     $               c, ldc, scale, info )
        RETURN
      END IF
*
*     Set constants to control overflow
*
      smlnum = dlamch( 'S' )
      bignum = one / smlnum
*
*      Partition A such that 2-by-2 blocks on the diagonal are not split
*
       skip = .false.
       DO i = 1, nba
          iwork( i ) = ( i - 1 ) * nb + 1
       END DO
       iwork( nba + 1 ) = m + 1
       DO k = 1, nba
          l1 = iwork( k )
          l2 = iwork( k + 1 ) - 1
          DO l = l1, l2
             IF( skip ) THEN
                skip = .false.
                cycle
             END IF
             IF( l.GE.m ) THEN
*               A( M, M ) is a 1-by-1 block
                cycle
             END IF
             IF( a( l, l+1 ).NE.zero .AND. a( l+1, l ).NE.zero ) THEN
*               Check if 2-by-2 block is split
                IF( l + 1 .EQ. iwork( k + 1 ) ) THEN
                   iwork( k + 1 ) = iwork( k + 1 ) + 1
                   cycle
                END IF
                skip = .true.
             END IF
          END DO
       END DO
       iwork( nba + 1 ) = m + 1
       IF( iwork( nba ).GE.iwork( nba + 1 ) ) THEN
          iwork( nba ) = iwork( nba + 1 )
          nba = nba - 1
       END IF
*
*      Partition B such that 2-by-2 blocks on the diagonal are not split
*
       pc = nba + 1
       skip = .false.
       DO i = 1, nbb
          iwork( pc + i ) = ( i - 1 ) * nb + 1
       END DO
       iwork( pc + nbb + 1 ) = n + 1
       DO k = 1, nbb
          l1 = iwork( pc + k )
          l2 = iwork( pc + k + 1 ) - 1
          DO l = l1, l2
             IF( skip ) THEN
                skip = .false.
                cycle
             END IF
             IF( l.GE.n ) THEN
*               B( N, N ) is a 1-by-1 block
                cycle
             END IF
             IF( b( l, l+1 ).NE.zero .AND. b( l+1, l ).NE.zero ) THEN
*               Check if 2-by-2 block is split
                IF( l + 1 .EQ. iwork( pc + k + 1 ) ) THEN
                   iwork( pc + k + 1 ) = iwork( pc + k + 1 ) + 1
                   cycle
                END IF
                skip = .true.
             END IF
          END DO
       END DO
       iwork( pc + nbb + 1 ) = n + 1
       IF( iwork( pc + nbb ).GE.iwork( pc + nbb + 1 ) ) THEN
          iwork( pc + nbb ) = iwork( pc + nbb + 1 )
          nbb = nbb - 1
       END IF
*
*     Set local scaling factors - must never attain zero.
*
      DO l = 1, nbb
         DO k = 1, nba
            swork( k, l ) = one
         END DO
      END DO
*
*     Fallback scaling factor to prevent flushing of SWORK( K, L ) to zero.
*     This scaling is to ensure compatibility with TRSYL and may get flushed.
*
      buf = one
*
*     Compute upper bounds of blocks of A and B
*
      awrk = nbb
      DO k = 1, nba
         k1 = iwork( k )
         k2 = iwork( k + 1 )
         DO l = k, nba
            l1 = iwork( l )
            l2 = iwork( l + 1 )
            IF( notrna ) THEN
               swork( k, awrk + l ) = dlange( 'I', k2-k1, l2-l1,
     $                                        a( k1, l1 ), lda, wnrm )
            ELSE
               swork( l, awrk + k ) = dlange( '1', k2-k1, l2-l1,
     $                                        a( k1, l1 ), lda, wnrm )
            END IF
         END DO
      END DO
      bwrk = nbb + nba
      DO k = 1, nbb
         k1 = iwork( pc + k )
         k2 = iwork( pc + k + 1 )
         DO l = k, nbb
            l1 = iwork( pc + l )
            l2 = iwork( pc + l + 1 )
            IF( notrnb ) THEN
               swork( k, bwrk + l ) = dlange( 'I', k2-k1, l2-l1,
     $                                        b( k1, l1 ), ldb, wnrm )
            ELSE
               swork( l, bwrk + k ) = dlange( '1', k2-k1, l2-l1,
     $                                        b( k1, l1 ), ldb, wnrm )
            END IF
         END DO
      END DO
*
      sgn = dble( isgn )
*
      IF( notrna .AND. notrnb ) THEN
*
*        Solve    A*X + ISGN*X*B = scale*C.
*
*        The (K,L)th block of X is determined starting from
*        bottom-left corner column by column by
*
*         A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
*
*        Where
*                  M                         L-1
*        R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
*                I=K+1                       J=1
*
*        Start loop over block rows (index = K) and block columns (index = L)
*
         DO k = nba, 1, -1
*
*           K1: row index of the first row in X( K, L )
*           K2: row index of the first row in X( K+1, L )
*           so the K2 - K1 is the column count of the block X( K, L )
*
            k1 = iwork( k )
            k2 = iwork( k + 1 )
            DO l = 1, nbb
*
*              L1: column index of the first column in X( K, L )
*              L2: column index of the first column in X( K, L + 1)
*              so that L2 - L1 is the row count of the block X( K, L )
*
               l1 = iwork( pc + l )
               l2 = iwork( pc + l + 1 )
*
               CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
     $                      a( k1, k1 ), lda,
     $                      b( l1, l1 ), ldb,
     $                      c( k1, l1 ), ldc, scaloc, iinfo )
               info = max( info, iinfo )
*
               IF ( scaloc * swork( k, l ) .EQ. zero ) THEN
                  IF( scaloc .EQ. zero ) THEN
*                    The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
*                    is larger than the product of BIGNUM**2 and cannot be
*                    represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
*                    Mark the computation as pointless.
                     buf = zero
                  ELSE
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                  END IF
                  DO jj = 1, nbb
                     DO ll = 1, nba
*                       Bound by BIGNUM to not introduce Inf. The value
*                       is irrelevant; corresponding entries of the
*                       solution will be flushed in consistency scaling.
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               swork( k, l ) = scaloc * swork( k, l )
               xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
     $                        wnrm )
*
               DO i = k - 1, 1, -1
*
*                 C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
*
                  i1 = iwork( i )
                  i2 = iwork( i + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( i, l ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( i, l ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  anrm = swork( i, awrk + k )
                  scaloc = dlarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to C( I, L ) and C( K, L ).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                      DO jj = l1, l2-1
                         CALL dscal( k2-k1, scal, c( k1, jj ), 1)
                      END DO
                  ENDIF
*
                  scal = ( scamin / swork( i, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                      DO ll = l1, l2-1
                         CALL dscal( i2-i1, scal, c( i1, ll ), 1)
                      END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( i, l ) = scamin * scaloc
*
                  CALL dgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -one,
     $                        a( i1, k1 ), lda, c( k1, l1 ), ldc,
     $                        one, c( i1, l1 ), ldc )
*
               END DO
*
               DO j = l + 1, nbb
*
*                 C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
*
                  j1 = iwork( pc + j )
                  j2 = iwork( pc + j + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( k, j ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( k, j ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  bnrm = swork(l, bwrk + j)
                  scaloc = dlarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to C( K, J ) and C( K, L).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF( scal .NE. one ) THEN
                     DO ll = l1, l2-1
                        CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
                     END DO
                  ENDIF
*
                  scal = ( scamin / swork( k, j ) ) * scaloc
                  IF( scal .NE. one ) THEN
                      DO jj = j1, j2-1
                         CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
                      END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( k, j ) = scamin * scaloc
*
                  CALL dgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -sgn,
     $                        c( k1, l1 ), ldc, b( l1, j1 ), ldb,
     $                        one, c( k1, j1 ), ldc )
               END DO
            END DO
         END DO
      ELSE IF( .NOT.notrna .AND. notrnb ) THEN
*
*        Solve    A**T*X + ISGN*X*B = scale*C.
*
*        The (K,L)th block of X is determined starting from
*        upper-left corner column by column by
*
*          A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
*
*        Where
*                   K-1                        L-1
*          R(K,L) = SUM [A(I,K)**T*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
*                   I=1                        J=1
*
*        Start loop over block rows (index = K) and block columns (index = L)
*
         DO k = 1, nba
*
*           K1: row index of the first row in X( K, L )
*           K2: row index of the first row in X( K+1, L )
*           so the K2 - K1 is the column count of the block X( K, L )
*
            k1 = iwork( k )
            k2 = iwork( k + 1 )
            DO l = 1, nbb
*
*              L1: column index of the first column in X( K, L )
*              L2: column index of the first column in X( K, L + 1)
*              so that L2 - L1 is the row count of the block X( K, L )
*
               l1 = iwork( pc + l )
               l2 = iwork( pc + l + 1 )
*
               CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
     $                      a( k1, k1 ), lda,
     $                      b( l1, l1 ), ldb,
     $                      c( k1, l1 ), ldc, scaloc, iinfo )
               info = max( info, iinfo )
*
               IF( scaloc * swork( k, l ) .EQ. zero ) THEN
                  IF( scaloc .EQ. zero ) THEN
*                    The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
*                    is larger than the product of BIGNUM**2 and cannot be
*                    represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
*                    Mark the computation as pointless.
                     buf = zero
                  ELSE
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                  END IF
                  DO jj = 1, nbb
                     DO ll = 1, nba
*                       Bound by BIGNUM to not introduce Inf. The value
*                       is irrelevant; corresponding entries of the
*                       solution will be flushed in consistency scaling.
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               swork( k, l ) = scaloc * swork( k, l )
               xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
     $                        wnrm )
*
               DO i = k + 1, nba
*
*                 C( I, L ) := C( I, L ) - A( K, I )**T * C( K, L )
*
                  i1 = iwork( i )
                  i2 = iwork( i + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( i, l ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( i, l ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  anrm = swork( i, awrk + k )
                  scaloc = dlarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to to C( I, L ) and C( K, L ).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                     DO ll = l1, l2-1
                        CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
                     END DO
                  ENDIF
*
                  scal = ( scamin / swork( i, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                     DO ll = l1, l2-1
                        CALL dscal( i2-i1, scal, c( i1, ll ), 1 )
                     END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( i, l ) = scamin * scaloc
*
                  CALL dgemm( 'T', 'N', i2-i1, l2-l1, k2-k1, -one,
     $                        a( k1, i1 ), lda, c( k1, l1 ), ldc,
     $                        one, c( i1, l1 ), ldc )
               END DO
*
               DO j = l + 1, nbb
*
*                 C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
*
                  j1 = iwork( pc + j )
                  j2 = iwork( pc + j + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( k, j ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( k, j ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  bnrm = swork( l, bwrk + j )
                  scaloc = dlarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to to C( K, J ) and C( K, L ).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF( scal .NE. one ) THEN
                      DO ll = l1, l2-1
                         CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
                      END DO
                  ENDIF
*
                  scal = ( scamin / swork( k, j ) ) * scaloc
                  IF( scal .NE. one ) THEN
                     DO jj = j1, j2-1
                        CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
                     END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( k, j ) = scamin * scaloc
*
                  CALL dgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -sgn,
     $                        c( k1, l1 ), ldc, b( l1, j1 ), ldb,
     $                        one, c( k1, j1 ), ldc )
               END DO
            END DO
         END DO
      ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
*
*        Solve    A**T*X + ISGN*X*B**T = scale*C.
*
*        The (K,L)th block of X is determined starting from
*        top-right corner column by column by
*
*           A(K,K)**T*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
*
*        Where
*                     K-1                          N
*            R(K,L) = SUM [A(I,K)**T*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
*                     I=1                        J=L+1
*
*        Start loop over block rows (index = K) and block columns (index = L)
*
         DO k = 1, nba
*
*           K1: row index of the first row in X( K, L )
*           K2: row index of the first row in X( K+1, L )
*           so the K2 - K1 is the column count of the block X( K, L )
*
            k1 = iwork( k )
            k2 = iwork( k + 1 )
            DO l = nbb, 1, -1
*
*              L1: column index of the first column in X( K, L )
*              L2: column index of the first column in X( K, L + 1)
*              so that L2 - L1 is the row count of the block X( K, L )
*
               l1 = iwork( pc + l )
               l2 = iwork( pc + l + 1 )
*
               CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
     $                      a( k1, k1 ), lda,
     $                      b( l1, l1 ), ldb,
     $                      c( k1, l1 ), ldc, scaloc, iinfo )
               info = max( info, iinfo )
*
               swork( k, l ) = scaloc * swork( k, l )
               IF( scaloc * swork( k, l ) .EQ. zero ) THEN
                  IF( scaloc .EQ. zero ) THEN
*                    The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
*                    is larger than the product of BIGNUM**2 and cannot be
*                    represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
*                    Mark the computation as pointless.
                     buf = zero
                  ELSE
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                  END IF
                  DO jj = 1, nbb
                     DO ll = 1, nba
*                       Bound by BIGNUM to not introduce Inf. The value
*                       is irrelevant; corresponding entries of the
*                       solution will be flushed in consistency scaling.
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
     $                        wnrm )
*
               DO i = k + 1, nba
*
*                 C( I, L ) := C( I, L ) - A( K, I )**T * C( K, L )
*
                  i1 = iwork( i )
                  i2 = iwork( i + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( i, l ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( i, l ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  anrm = swork( i, awrk + k )
                  scaloc = dlarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to C( I, L ) and C( K, L ).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                     DO ll = l1, l2-1
                        CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
                     END DO
                  ENDIF
*
                  scal = ( scamin / swork( i, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                     DO ll = l1, l2-1
                        CALL dscal( i2-i1, scal, c( i1, ll ), 1 )
                     END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( i, l ) = scamin * scaloc
*
                  CALL dgemm( 'T', 'N', i2-i1, l2-l1, k2-k1, -one,
     $                        a( k1, i1 ), lda, c( k1, l1 ), ldc,
     $                        one, c( i1, l1 ), ldc )
               END DO
*
               DO j = 1, l - 1
*
*                 C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**T
*
                  j1 = iwork( pc + j )
                  j2 = iwork( pc + j + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( k, j ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( k, j ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  bnrm = swork( l, bwrk + j )
                  scaloc = dlarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to C( K, J ) and C( K, L ).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF( scal .NE. one ) THEN
                     DO ll = l1, l2-1
                        CALL dscal( k2-k1, scal, c( k1, ll ), 1)
                     END DO
                  ENDIF
*
                  scal = ( scamin / swork( k, j ) ) * scaloc
                  IF( scal .NE. one ) THEN
                     DO jj = j1, j2-1
                        CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
                     END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( k, j ) = scamin * scaloc
*
                  CALL dgemm( 'N', 'T', k2-k1, j2-j1, l2-l1, -sgn,
     $                        c( k1, l1 ), ldc, b( j1, l1 ), ldb,
     $                        one, c( k1, j1 ), ldc )
               END DO
            END DO
         END DO
      ELSE IF( notrna .AND. .NOT.notrnb ) THEN
*
*        Solve    A*X + ISGN*X*B**T = scale*C.
*
*        The (K,L)th block of X is determined starting from
*        bottom-right corner column by column by
*
*            A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**T = C(K,L) - R(K,L)
*
*        Where
*                      M                          N
*            R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**T].
*                    I=K+1                      J=L+1
*
*        Start loop over block rows (index = K) and block columns (index = L)
*
         DO k = nba, 1, -1
*
*           K1: row index of the first row in X( K, L )
*           K2: row index of the first row in X( K+1, L )
*           so the K2 - K1 is the column count of the block X( K, L )
*
            k1 = iwork( k )
            k2 = iwork( k + 1 )
            DO l = nbb, 1, -1
*
*              L1: column index of the first column in X( K, L )
*              L2: column index of the first column in X( K, L + 1)
*              so that L2 - L1 is the row count of the block X( K, L )
*
               l1 = iwork( pc + l )
               l2 = iwork( pc + l + 1 )
*
               CALL dtrsyl( trana, tranb, isgn, k2-k1, l2-l1,
     $                      a( k1, k1 ), lda,
     $                      b( l1, l1 ), ldb,
     $                      c( k1, l1 ), ldc, scaloc, iinfo )
               info = max( info, iinfo )
*
               IF( scaloc * swork( k, l ) .EQ. zero ) THEN
                  IF( scaloc .EQ. zero ) THEN
*                    The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
*                    is larger than the product of BIGNUM**2 and cannot be
*                    represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
*                    Mark the computation as pointless.
                     buf = zero
                  ELSE
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                  END IF
                  DO jj = 1, nbb
                     DO ll = 1, nba
*                       Bound by BIGNUM to not introduce Inf. The value
*                       is irrelevant; corresponding entries of the
*                       solution will be flushed in consistency scaling.
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               swork( k, l ) = scaloc * swork( k, l )
               xnrm = dlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
     $                        wnrm )
*
               DO i = 1, k - 1
*
*                 C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
*
                  i1 = iwork( i )
                  i2 = iwork( i + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', i2-i1, l2-l1, c( i1, l1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( i, l ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( i, l ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  anrm = swork( i, awrk + k )
                  scaloc = dlarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to C( I, L ) and C( K, L ).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                     DO ll = l1, l2-1
                        CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
                     END DO
                  ENDIF
*
                  scal = ( scamin / swork( i, l ) ) * scaloc
                  IF (scal .NE. one) THEN
                     DO ll = l1, l2-1
                        CALL dscal( i2-i1, scal, c( i1, ll ), 1 )
                     END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( i, l ) = scamin * scaloc
*
                  CALL dgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -one,
     $                        a( i1, k1 ), lda, c( k1, l1 ), ldc,
     $                        one, c( i1, l1 ), ldc )
*
               END DO
*
               DO j = 1, l - 1
*
*                 C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**T
*
                  j1 = iwork( pc + j )
                  j2 = iwork( pc + j + 1 )
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = dlange( 'I', k2-k1, j2-j1, c( k1, j1 ),
     $                           ldc, wnrm )
                  scamin = min( swork( k, j ), swork( k, l ) )
                  cnrm = cnrm * ( scamin / swork( k, j ) )
                  xnrm = xnrm * ( scamin / swork( k, l ) )
                  bnrm = swork( l, bwrk + j )
                  scaloc = dlarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.d0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.d0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.d0**exponent( scaloc )
                     scaloc = scaloc / 2.d0**exponent( scaloc )
                  END IF
                  cnrm = cnrm * scaloc
                  xnrm = xnrm * scaloc
*
*                 Simultaneously apply the robust update factor and the
*                 consistency scaling factor to C( K, J ) and C( K, L ).
*
                  scal = ( scamin / swork( k, l ) ) * scaloc
                  IF( scal .NE. one ) THEN
                     DO jj = l1, l2-1
                        CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
                     END DO
                  ENDIF
*
                  scal = ( scamin / swork( k, j ) ) * scaloc
                  IF( scal .NE. one ) THEN
                     DO jj = j1, j2-1
                        CALL dscal( k2-k1, scal, c( k1, jj ), 1 )
                     END DO
                  ENDIF
*
*                 Record current scaling factor
*
                  swork( k, l ) = scamin * scaloc
                  swork( k, j ) = scamin * scaloc
*
                  CALL dgemm( 'N', 'T', k2-k1, j2-j1, l2-l1, -sgn,
     $                        c( k1, l1 ), ldc, b( j1, l1 ), ldb,
     $                        one, c( k1, j1 ), ldc )
               END DO
            END DO
         END DO
*
      END IF
*
*     Reduce local scaling factors
*
      scale = swork( 1, 1 )
      DO k = 1, nba
         DO l = 1, nbb
            scale = min( scale, swork( k, l ) )
         END DO
      END DO
*
      IF( scale .EQ. zero ) THEN
*
*        The magnitude of the largest entry of the solution is larger
*        than the product of BIGNUM**2 and cannot be represented in the
*        form (1/SCALE)*X if SCALE is DOUBLE PRECISION. Set SCALE to
*        zero and give up.
*
         iwork(1) = nba + nbb + 2
         swork(1,1) = max( nba, nbb )
         swork(2,1) = 2 * nbb + nba
         RETURN
      END IF
*
*     Realize consistent scaling
*
      DO k = 1, nba
         k1 = iwork( k )
         k2 = iwork( k + 1 )
         DO l = 1, nbb
            l1 = iwork( pc + l )
            l2 = iwork( pc + l + 1 )
            scal = scale / swork( k, l )
            IF( scal .NE. one ) THEN
               DO ll = l1, l2-1
                  CALL dscal( k2-k1, scal, c( k1, ll ), 1 )
               END DO
            ENDIF
         END DO
      END DO
*
      IF( buf .NE. one .AND. buf.GT.zero ) THEN
*
*        Decrease SCALE as much as possible.
*
         scaloc = min( scale / smlnum, one / buf )
         buf = buf * scaloc
         scale = scale / scaloc
      END IF
 
      IF( buf.NE.one .AND. buf.GT.zero ) THEN
*
*        In case of overly aggressive scaling during the computation,
*        flushing of the global scale factor may be prevented by
*        undoing some of the scaling. This step is to ensure that
*        this routine flushes only scale factors that TRSYL also
*        flushes and be usable as a drop-in replacement.
*
*        How much can the normwise largest entry be upscaled?
*
         scal = c( 1, 1 )
         DO k = 1, m
            DO l = 1, n
               scal = max( scal, abs( c( k, l ) ) )
            END DO
         END DO
*
*        Increase BUF as close to 1 as possible and apply scaling.
*
         scaloc = min( bignum / scal, one / buf )
         buf = buf * scaloc
         CALL dlascl( 'G', -1, -1, one, scaloc, m, n, c, ldc,
     $               iwork(1) )
      END IF
*
*     Combine with buffer scaling factor. SCALE will be flushed if
*     BUF is less than one here.
*
      scale = scale * buf
*
*     Restore workspace dimensions
*
      iwork(1) = nba + nbb + 2
      swork(1,1) = max( nba, nbb )
      swork(2,1) = 2 * nbb + nba
*
      RETURN
*
*     End of DTRSYL3
*

Here is the call graph for this function:

Here is the caller graph for this function: