strsyl3.f

*> \brief \b STRSYL3
*
* Definition:
* ===========
*
*
*>  \par Purpose
*  =============
*>
*> \verbatim
*>
*>  STRSYL3 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.
*>
*>  This is the block version of the algorithm.
*> \endverbatim
*
*  Arguments
*  =========
*
*> \param[in] TRANA
*> \verbatim
*>          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)
*> \endverbatim
*>
*> \param[in] TRANB
*> \verbatim
*>          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)
*> \endverbatim
*>
*> \param[in] ISGN
*> \verbatim
*>          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
*> \endverbatim
*>
*> \param[in] M
*> \verbatim
*>          M is INTEGER
*>          The order of the matrix A, and the number of rows in the
*>          matrices X and C. M >= 0.
*> \endverbatim
*>
*> \param[in] N
*> \verbatim
*>          N is INTEGER
*>          The order of the matrix B, and the number of columns in the
*>          matrices X and C. N >= 0.
*> \endverbatim
*>
*> \param[in] A
*> \verbatim
*>          A is REAL array, dimension (LDA,M)
*>          The upper quasi-triangular matrix A, in Schur canonical form.
*> \endverbatim
*>
*> \param[in] LDA
*> \verbatim
*>          LDA is INTEGER
*>          The leading dimension of the array A. LDA >= max(1,M).
*> \endverbatim
*>
*> \param[in] B
*> \verbatim
*>          B is REAL array, dimension (LDB,N)
*>          The upper quasi-triangular matrix B, in Schur canonical form.
*> \endverbatim
*>
*> \param[in] LDB
*> \verbatim
*>          LDB is INTEGER
*>          The leading dimension of the array B. LDB >= max(1,N).
*> \endverbatim
*>
*> \param[in,out] C
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[in] LDC
*> \verbatim
*>          LDC is INTEGER
*>          The leading dimension of the array C. LDC >= max(1,M)
*> \endverbatim
*>
*> \param[out] SCALE
*> \verbatim
*>          SCALE is REAL
*>          The scale factor, scale, set <= 1 to avoid overflow in X.
*> \endverbatim
*>
*> \param[out] IWORK
*> \verbatim
*>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
*>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
*> \endverbatim
*>
*> \param[in] LIWORK
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] SWORK
*> \verbatim
*>          SWORK is REAL 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.
*> \endverbatim
*>
*> \param[in] LDSWORK
*> \verbatim
*>          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.
*> \endverbatim
*>
*> \param[out] INFO
*> \verbatim
*>          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).
*> \endverbatim
*
*  =====================================================================
*  References:
*   E. S. Quintana-Orti and R. A. Van De Geijn (2003). Formal derivation of
*   algorithms: The triangular Sylvester equation, ACM Transactions
*   on Mathematical Software (TOMS), volume 29, pages 218--243.
*
*   A. Schwarz and C. C. Kjelgaard Mikkelsen (2020). Robust Task-Parallel
*   Solution of the Triangular Sylvester Equation. Lecture Notes in
*   Computer Science, vol 12043, pages 82--92, Springer.
*
*  Contributor:
*   Angelika Schwarz, Umea University, Sweden.
*
*  =====================================================================
      SUBROUTINE strsyl3( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
     $                    LDC, SCALE, IWORK, LIWORK, SWORK, LDSWORK,
     $                    INFO )
      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          trana, tranb
      INTEGER            info, isgn, lda, ldb, ldc, m, n,
     $                   liwork, ldswork
      REAL               scale
*     ..
*     .. Array Arguments ..
      INTEGER            iwork( * )
      REAL               a( lda, * ), b( ldb, * ), c( ldc, * ),
     $                   swork( ldswork, * )
*     ..
*     .. Parameters ..
      REAL               zero, one
      parameter( zero = 0.0e+0, one = 1.0e+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
      REAL               anrm, bignum, bnrm, cnrm, scal, scaloc,
     $                   scamin, sgn, xnrm, buf, smlnum
*     ..
*     .. Local Arrays ..
      REAL               wnrm( max( m, n ) )
*     ..
*     .. External Functions ..
      LOGICAL            lsame
      INTEGER            ilaenv
      REAL               slange, slamch, slarmm
      EXTERNAL           slange, slamch, slarmm, ilaenv, lsame
*     ..
*     .. External Subroutines ..
      EXTERNAL           sgemm, slascl, sscal, strsyl, xerbla
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, exponent, max, min, real
*     ..
*     .. 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, 'STRSYL', '', 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
      ELSE IF( .NOT.lquery .AND. liwork.LT.iwork(1) ) THEN
         info = -14
      ELSE IF( .NOT.lquery .AND. ldswork.LT.max( nba, nbb ) ) THEN
         info = -16
      END IF
      IF( info.NE.0 ) THEN
         CALL xerbla( 'STRSYL3', -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 strsyl( trana, tranb, isgn, m, n, a, lda, b, ldb,
     $               c, ldc, scale, info )
        RETURN
      END IF
*
*     Set constants to control overflow
*
      smlnum = slamch( '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 ) = slange( 'I', k2-k1, l2-l1,
     $                                        a( k1, l1 ), lda, wnrm )
            ELSE
               swork( l, awrk + k ) = slange( '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 ) = slange( 'I', k2-k1, l2-l1,
     $                                        b( k1, l1 ), ldb, wnrm )
            ELSE
               swork( l, bwrk + k ) = slange( '1', k2-k1, l2-l1,
     $                                        b( k1, l1 ), ldb, wnrm )
            END IF
         END DO
      END DO
*
      sgn = REAL( 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 strsyl( 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.e0**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.e0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               swork( k, l ) = scaloc * swork( k, l )
               xnrm = slange( '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 = slange( '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 = slarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 = slange( '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 = slarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 strsyl( 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.e0**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.e0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               swork( k, l ) = scaloc * swork( k, l )
               xnrm = slange( '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 = slange( '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 = slarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 = slange( '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 = slarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 strsyl( 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.e0**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.e0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               swork( k, l ) = scaloc * swork( k, l )
               xnrm = slange( '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 = slange( '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 = slarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 = slange( '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 = slarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 strsyl( 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.e0**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.e0**exponent( scaloc ) )
                     END DO
                  END DO
               END IF
               swork( k, l ) = scaloc * swork( k, l )
               xnrm = slange( '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 = slange( '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 = slarmm( anrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 = slange( '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 = slarmm( bnrm, xnrm, cnrm )
                  IF( scaloc * scamin .EQ. zero ) THEN
*                    Use second scaling factor to prevent flushing to zero.
                     buf = buf*2.e0**exponent( scaloc )
                     DO jj = 1, nbb
                        DO ll = 1, nba
                        swork( ll, jj ) = min( bignum,
     $                     swork( ll, jj ) / 2.e0**exponent( scaloc ) )
                        END DO
                     END DO
                     scamin = scamin / 2.e0**exponent( scaloc )
                     scaloc = scaloc / 2.e0**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 sscal( 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 sscal( 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 sgemm( '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 REAL. 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 sscal( 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 slascl( 'G', -1, -1, one, scaloc, m, n, c, ldc, iwork )
      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 STRSYL3
*
      END