ztrsyl3()

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

ZTRSYL3

Purpose:

!>
!>  ZTRSYL3 solves the complex 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**H, and  A and B are both upper 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.
!>
!>  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) !> = 'C': op(A) = A**H (Conjugate transpose) !>
[in]	TRANB	!> TRANB is CHARACTER*1 !> Specifies the option op(B): !> = 'N': op(B) = B (No transpose) !> = 'C': op(B) = B**H (Conjugate 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 COMPLEX*16 array, dimension (LDA,M) !> The upper triangular matrix A. !>
[in]	LDA	!> LDA is INTEGER !> The leading dimension of the array A. LDA >= max(1,M). !>
[in]	B	!> B is COMPLEX*16 array, dimension (LDB,N) !> The upper triangular matrix B. !>
[in]	LDB	!> LDB is INTEGER !> The leading dimension of the array B. LDB >= max(1,N). !>
[in,out]	C	!> C is COMPLEX*16 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]	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 168 of file ztrsyl3.f.

      IMPLICIT NONE
*
*     .. Scalar Arguments ..
      CHARACTER          TRANA, TRANB
      INTEGER            INFO, ISGN, LDA, LDB, LDC, LDSWORK, M, N
      DOUBLE PRECISION   SCALE
*     ..
*     .. Array Arguments ..
      COMPLEX*16         A( LDA, * ), B( LDB, * ), C( LDC, * )
      DOUBLE PRECISION   SWORK( LDSWORK, * )
*     ..
*     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE
      parameter( zero = 0.0d0, one = 1.0d0 )
      COMPLEX*16         CONE
      parameter( cone = ( 1.0d0, 0.0d0 ) )
*     ..
*     .. Local Scalars ..
      LOGICAL            NOTRNA, NOTRNB, LQUERY
      INTEGER            AWRK, BWRK, I, I1, I2, IINFO, J, J1, J2, JJ,
     $                   K, K1, K2, L, L1, L2, LL, NBA, NB, NBB
      DOUBLE PRECISION   ANRM, BIGNUM, BNRM, CNRM, SCAL, SCALOC,
     $                   SCAMIN, SGN, XNRM, BUF, SMLNUM
      COMPLEX*16         CSGN
*     ..
*     .. Local Arrays ..
      DOUBLE PRECISION   WNRM( MAX( M, N ) )
*     ..
*     .. External Functions ..
      LOGICAL            LSAME
      INTEGER            ILAENV
      DOUBLE PRECISION   DLAMCH, DLARMM, ZLANGE
      EXTERNAL           dlamch, dlarmm, ilaenv, lsame,
     $                   zlange
*     ..
*     .. External Subroutines ..
      EXTERNAL           xerbla, zdscal, zgemm, zlascl,
     $                   ztrsyl
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          abs, dble, dimag, 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, 'ZTRSYL', '', 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 = ( ldswork.EQ.-1 )
      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, 'C' ) ) THEN
         info = -1
      ELSE IF( .NOT.notrnb .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( 'ZTRSYL3', -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
*     workspace is provided
*
      IF( min( nba, nbb ).EQ.1 .OR. ldswork.LT.max( nba, nbb ) ) THEN
        CALL ztrsyl( 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
*
*     Set local scaling factors.
*
      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 = (k - 1) * nb + 1
         k2 = min( k * nb, m ) + 1
         DO l = k, nba
            l1 = (l - 1) * nb + 1
            l2 = min( l * nb, m ) + 1
            IF( notrna ) THEN
               swork( k, awrk + l ) = zlange( 'I', k2-k1, l2-l1,
     $                                        a( k1, l1 ), lda, wnrm )
            ELSE
               swork( l, awrk + k ) = zlange( '1', k2-k1, l2-l1,
     $                                        a( k1, l1 ), lda, wnrm )
            END IF
         END DO
      END DO
      bwrk = nbb + nba
      DO k = 1, nbb
         k1 = (k - 1) * nb + 1
         k2 = min( k * nb, n ) + 1
         DO l = k, nbb
            l1 = (l - 1) * nb + 1
            l2 = min( l * nb, n ) + 1
            IF( notrnb ) THEN
               swork( k, bwrk + l ) = zlange( 'I', k2-k1, l2-l1,
     $                                        b( k1, l1 ), ldb, wnrm )
            ELSE
               swork( l, bwrk + k ) = zlange( '1', k2-k1, l2-l1,
     $                                        b( k1, l1 ), ldb, wnrm )
            END IF
         END DO
      END DO
*
      sgn = dble( isgn )
      csgn = dcmplx( sgn, zero )
*
      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 = (k - 1) * nb + 1
            k2 = min( k * nb, m ) + 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 = (l - 1) * nb + 1
               l2 = min( l * nb, n ) + 1
*
               CALL ztrsyl( 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
                     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 = zlange( '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 = (i - 1) * nb + 1
                  i2 = min( i * nb, m ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -cone,
     $                        a( i1, k1 ), lda, c( k1, l1 ), ldc,
     $                        cone, 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 = (j - 1) * nb + 1
                  j2 = min( j * nb, n ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -csgn,
     $                        c( k1, l1 ), ldc, b( l1, j1 ), ldb,
     $                        cone, c( k1, j1 ), ldc )
               END DO
            END DO
         END DO
      ELSE IF( .NOT.notrna .AND. notrnb ) THEN
*
*        Solve    A**H *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)**H*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)**H*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 = (k - 1) * nb + 1
            k2 = min( k * nb, m ) + 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 = (l - 1) * nb + 1
               l2 = min( l * nb, n ) + 1
*
               CALL ztrsyl( 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 = zlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
     $                        wnrm )
*
               DO i = k + 1, nba
*
*                 C( I, L ) := C( I, L ) - A( K, I )**H * C( K, L )
*
                  i1 = (i - 1) * nb + 1
                  i2 = min( i * nb, m ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'C', 'N', i2-i1, l2-l1, k2-k1, -cone,
     $                        a( k1, i1 ), lda, c( k1, l1 ), ldc,
     $                        cone, 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 = (j - 1) * nb + 1
                  j2 = min( j * nb, n ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'N', 'N', k2-k1, j2-j1, l2-l1, -csgn,
     $                        c( k1, l1 ), ldc, b( l1, j1 ), ldb,
     $                        cone, c( k1, j1 ), ldc )
               END DO
            END DO
         END DO
      ELSE IF( .NOT.notrna .AND. .NOT.notrnb ) THEN
*
*        Solve    A**H *X + ISGN*X*B**H = scale*C.
*
*        The (K,L)th block of X is determined starting from
*        top-right corner column by column by
*
*           A(K,K)**H*X(K,L) + ISGN*X(K,L)*B(L,L)**H = C(K,L) - R(K,L)
*
*        Where
*                     K-1                          N
*            R(K,L) = SUM [A(I,K)**H*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**H].
*                     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 = (k - 1) * nb + 1
            k2 = min( k * nb, m ) + 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 = (l - 1) * nb + 1
               l2 = min( l * nb, n ) + 1
*
               CALL ztrsyl( 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 = zlange( 'I', k2-k1, l2-l1, c( k1, l1 ), ldc,
     $                        wnrm )
*
               DO i = k + 1, nba
*
*                 C( I, L ) := C( I, L ) - A( K, I )**H * C( K, L )
*
                  i1 = (i - 1) * nb + 1
                  i2 = min( i * nb, m ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'C', 'N', i2-i1, l2-l1, k2-k1, -cone,
     $                        a( k1, i1 ), lda, c( k1, l1 ), ldc,
     $                        cone, c( i1, l1 ), ldc )
               END DO
*
               DO j = 1, l - 1
*
*                 C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**H
*
                  j1 = (j - 1) * nb + 1
                  j2 = min( j * nb, n ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'N', 'C', k2-k1, j2-j1, l2-l1, -csgn,
     $                        c( k1, l1 ), ldc, b( j1, l1 ), ldb,
     $                        cone, c( k1, j1 ), ldc )
               END DO
            END DO
         END DO
      ELSE IF( notrna .AND. .NOT.notrnb ) THEN
*
*        Solve    A*X + ISGN*X*B**H = 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)**H = 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)**H].
*                    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 = (k - 1) * nb + 1
            k2 = min( k * nb, m ) + 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 = (l - 1) * nb + 1
               l2 = min( l * nb, n ) + 1
*
               CALL ztrsyl( 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 = zlange( '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 = (i - 1) * nb + 1
                  i2 = min( i * nb, m ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'N', 'N', i2-i1, l2-l1, k2-k1, -cone,
     $                        a( i1, k1 ), lda, c( k1, l1 ), ldc,
     $                        cone, c( i1, l1 ), ldc )
*
               END DO
*
               DO j = 1, l - 1
*
*                 C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**H
*
                  j1 = (j - 1) * nb + 1
                  j2 = min( j * nb, n ) + 1
*
*                 Compute scaling factor to survive the linear update
*                 simulating consistent scaling.
*
                  cnrm = zlange( '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 zdscal( 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 zdscal( 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 zgemm( 'N', 'C', k2-k1, j2-j1, l2-l1, -csgn,
     $                        c( k1, l1 ), ldc, b( j1, l1 ), ldb,
     $                        cone, 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.
*
         swork(1,1) = max( nba, nbb )
         swork(2,1) = 2 * nbb + nba
         RETURN
      END IF
*
*     Realize consistent scaling
*
      DO k = 1, nba
         k1 = (k - 1) * nb + 1
         k2 = min( k * nb, m ) + 1
         DO l = 1, nbb
            l1 = (l - 1) * nb + 1
            l2 = min( l * nb, n ) + 1
            scal = scale / swork( k, l )
            IF( scal .NE. one ) THEN
               DO ll = l1, l2-1
                  CALL zdscal( 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 = max( abs( dble( c( 1, 1 ) ) ),
     $               abs( dimag( c( 1, 1 ) ) ) )
         DO k = 1, m
            DO l = 1, n
               scal = max( scal, abs( dble( c( k, l ) ) ),
     $                     abs( dimag( 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 zlascl( 'G', -1, -1, one, scaloc, m, n, c, ldc, iinfo )
      END IF
*
*     Combine with buffer scaling factor. SCALE will be flushed if
*     BUF is less than one here.
*
      scale = scale * buf
*
*     Restore workspace dimensions
*
      swork(1,1) = max( nba, nbb )
      swork(2,1) = 2 * nbb + nba
*
      RETURN
*
*     End of ZTRSYL3
*

Here is the call graph for this function:

Here is the caller graph for this function: