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 178 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: