Annotation of rpl/lapack/lapack/ztrsyl3.f, revision 1.1
1.1 ! bertrand 1: *> \brief \b ZTRSYL3
! 2: *
! 3: * Definition:
! 4: * ===========
! 5: *
! 6: *
! 7: *> \par Purpose
! 8: * =============
! 9: *>
! 10: *> \verbatim
! 11: *>
! 12: *> ZTRSYL3 solves the complex Sylvester matrix equation:
! 13: *>
! 14: *> op(A)*X + X*op(B) = scale*C or
! 15: *> op(A)*X - X*op(B) = scale*C,
! 16: *>
! 17: *> where op(A) = A or A**H, and A and B are both upper triangular. A is
! 18: *> M-by-M and B is N-by-N; the right hand side C and the solution X are
! 19: *> M-by-N; and scale is an output scale factor, set <= 1 to avoid
! 20: *> overflow in X.
! 21: *>
! 22: *> This is the block version of the algorithm.
! 23: *> \endverbatim
! 24: *
! 25: * Arguments
! 26: * =========
! 27: *
! 28: *> \param[in] TRANA
! 29: *> \verbatim
! 30: *> TRANA is CHARACTER*1
! 31: *> Specifies the option op(A):
! 32: *> = 'N': op(A) = A (No transpose)
! 33: *> = 'C': op(A) = A**H (Conjugate transpose)
! 34: *> \endverbatim
! 35: *>
! 36: *> \param[in] TRANB
! 37: *> \verbatim
! 38: *> TRANB is CHARACTER*1
! 39: *> Specifies the option op(B):
! 40: *> = 'N': op(B) = B (No transpose)
! 41: *> = 'C': op(B) = B**H (Conjugate transpose)
! 42: *> \endverbatim
! 43: *>
! 44: *> \param[in] ISGN
! 45: *> \verbatim
! 46: *> ISGN is INTEGER
! 47: *> Specifies the sign in the equation:
! 48: *> = +1: solve op(A)*X + X*op(B) = scale*C
! 49: *> = -1: solve op(A)*X - X*op(B) = scale*C
! 50: *> \endverbatim
! 51: *>
! 52: *> \param[in] M
! 53: *> \verbatim
! 54: *> M is INTEGER
! 55: *> The order of the matrix A, and the number of rows in the
! 56: *> matrices X and C. M >= 0.
! 57: *> \endverbatim
! 58: *>
! 59: *> \param[in] N
! 60: *> \verbatim
! 61: *> N is INTEGER
! 62: *> The order of the matrix B, and the number of columns in the
! 63: *> matrices X and C. N >= 0.
! 64: *> \endverbatim
! 65: *>
! 66: *> \param[in] A
! 67: *> \verbatim
! 68: *> A is COMPLEX*16 array, dimension (LDA,M)
! 69: *> The upper triangular matrix A.
! 70: *> \endverbatim
! 71: *>
! 72: *> \param[in] LDA
! 73: *> \verbatim
! 74: *> LDA is INTEGER
! 75: *> The leading dimension of the array A. LDA >= max(1,M).
! 76: *> \endverbatim
! 77: *>
! 78: *> \param[in] B
! 79: *> \verbatim
! 80: *> B is COMPLEX*16 array, dimension (LDB,N)
! 81: *> The upper triangular matrix B.
! 82: *> \endverbatim
! 83: *>
! 84: *> \param[in] LDB
! 85: *> \verbatim
! 86: *> LDB is INTEGER
! 87: *> The leading dimension of the array B. LDB >= max(1,N).
! 88: *> \endverbatim
! 89: *>
! 90: *> \param[in,out] C
! 91: *> \verbatim
! 92: *> C is COMPLEX*16 array, dimension (LDC,N)
! 93: *> On entry, the M-by-N right hand side matrix C.
! 94: *> On exit, C is overwritten by the solution matrix X.
! 95: *> \endverbatim
! 96: *>
! 97: *> \param[in] LDC
! 98: *> \verbatim
! 99: *> LDC is INTEGER
! 100: *> The leading dimension of the array C. LDC >= max(1,M)
! 101: *> \endverbatim
! 102: *>
! 103: *> \param[out] SCALE
! 104: *> \verbatim
! 105: *> SCALE is DOUBLE PRECISION
! 106: *> The scale factor, scale, set <= 1 to avoid overflow in X.
! 107: *> \endverbatim
! 108: *>
! 109: *> \param[out] SWORK
! 110: *> \verbatim
! 111: *> SWORK is DOUBLE PRECISION array, dimension (MAX(2, ROWS),
! 112: *> MAX(1,COLS)).
! 113: *> On exit, if INFO = 0, SWORK(1) returns the optimal value ROWS
! 114: *> and SWORK(2) returns the optimal COLS.
! 115: *> \endverbatim
! 116: *>
! 117: *> \param[in] LDSWORK
! 118: *> \verbatim
! 119: *> LDSWORK is INTEGER
! 120: *> LDSWORK >= MAX(2,ROWS), where ROWS = ((M + NB - 1) / NB + 1)
! 121: *> and NB is the optimal block size.
! 122: *>
! 123: *> If LDSWORK = -1, then a workspace query is assumed; the routine
! 124: *> only calculates the optimal dimensions of the SWORK matrix,
! 125: *> returns these values as the first and second entry of the SWORK
! 126: *> matrix, and no error message related LWORK is issued by XERBLA.
! 127: *> \endverbatim
! 128: *>
! 129: *> \param[out] INFO
! 130: *> \verbatim
! 131: *> INFO is INTEGER
! 132: *> = 0: successful exit
! 133: *> < 0: if INFO = -i, the i-th argument had an illegal value
! 134: *> = 1: A and B have common or very close eigenvalues; perturbed
! 135: *> values were used to solve the equation (but the matrices
! 136: *> A and B are unchanged).
! 137: *> \endverbatim
! 138: *
! 139: *> \ingroup complex16SYcomputational
! 140: *
! 141: * =====================================================================
! 142: * References:
! 143: * E. S. Quintana-Orti and R. A. Van De Geijn (2003). Formal derivation of
! 144: * algorithms: The triangular Sylvester equation, ACM Transactions
! 145: * on Mathematical Software (TOMS), volume 29, pages 218--243.
! 146: *
! 147: * A. Schwarz and C. C. Kjelgaard Mikkelsen (2020). Robust Task-Parallel
! 148: * Solution of the Triangular Sylvester Equation. Lecture Notes in
! 149: * Computer Science, vol 12043, pages 82--92, Springer.
! 150: *
! 151: * Contributor:
! 152: * Angelika Schwarz, Umea University, Sweden.
! 153: *
! 154: * =====================================================================
! 155: SUBROUTINE ZTRSYL3( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C,
! 156: $ LDC, SCALE, SWORK, LDSWORK, INFO )
! 157: IMPLICIT NONE
! 158: *
! 159: * .. Scalar Arguments ..
! 160: CHARACTER TRANA, TRANB
! 161: INTEGER INFO, ISGN, LDA, LDB, LDC, LDSWORK, M, N
! 162: DOUBLE PRECISION SCALE
! 163: * ..
! 164: * .. Array Arguments ..
! 165: COMPLEX*16 A( LDA, * ), B( LDB, * ), C( LDC, * )
! 166: DOUBLE PRECISION SWORK( LDSWORK, * )
! 167: * ..
! 168: * .. Parameters ..
! 169: DOUBLE PRECISION ZERO, ONE
! 170: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 171: COMPLEX*16 CONE
! 172: PARAMETER ( CONE = ( 1.0D0, 0.0D0 ) )
! 173: * ..
! 174: * .. Local Scalars ..
! 175: LOGICAL NOTRNA, NOTRNB, LQUERY
! 176: INTEGER AWRK, BWRK, I, I1, I2, IINFO, J, J1, J2, JJ,
! 177: $ K, K1, K2, L, L1, L2, LL, NBA, NB, NBB
! 178: DOUBLE PRECISION ANRM, BIGNUM, BNRM, CNRM, SCAL, SCALOC,
! 179: $ SCAMIN, SGN, XNRM, BUF, SMLNUM
! 180: COMPLEX*16 CSGN
! 181: * ..
! 182: * .. Local Arrays ..
! 183: DOUBLE PRECISION WNRM( MAX( M, N ) )
! 184: * ..
! 185: * .. External Functions ..
! 186: LOGICAL LSAME
! 187: INTEGER ILAENV
! 188: DOUBLE PRECISION DLAMCH, DLARMM, ZLANGE
! 189: EXTERNAL DLAMCH, DLARMM, ILAENV, LSAME, ZLANGE
! 190: * ..
! 191: * .. External Subroutines ..
! 192: EXTERNAL XERBLA, ZDSCAL, ZGEMM, ZLASCL, ZTRSYL
! 193: * ..
! 194: * .. Intrinsic Functions ..
! 195: INTRINSIC ABS, DBLE, DIMAG, EXPONENT, MAX, MIN
! 196: * ..
! 197: * .. Executable Statements ..
! 198: *
! 199: * Decode and Test input parameters
! 200: *
! 201: NOTRNA = LSAME( TRANA, 'N' )
! 202: NOTRNB = LSAME( TRANB, 'N' )
! 203: *
! 204: * Use the same block size for all matrices.
! 205: *
! 206: NB = MAX( 8, ILAENV( 1, 'ZTRSYL', '', M, N, -1, -1) )
! 207: *
! 208: * Compute number of blocks in A and B
! 209: *
! 210: NBA = MAX( 1, (M + NB - 1) / NB )
! 211: NBB = MAX( 1, (N + NB - 1) / NB )
! 212: *
! 213: * Compute workspace
! 214: *
! 215: INFO = 0
! 216: LQUERY = ( LDSWORK.EQ.-1 )
! 217: IF( LQUERY ) THEN
! 218: LDSWORK = 2
! 219: SWORK(1,1) = MAX( NBA, NBB )
! 220: SWORK(2,1) = 2 * NBB + NBA
! 221: END IF
! 222: *
! 223: * Test the input arguments
! 224: *
! 225: IF( .NOT.NOTRNA .AND. .NOT. LSAME( TRANA, 'C' ) ) THEN
! 226: INFO = -1
! 227: ELSE IF( .NOT.NOTRNB .AND. .NOT. LSAME( TRANB, 'C' ) ) THEN
! 228: INFO = -2
! 229: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN
! 230: INFO = -3
! 231: ELSE IF( M.LT.0 ) THEN
! 232: INFO = -4
! 233: ELSE IF( N.LT.0 ) THEN
! 234: INFO = -5
! 235: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN
! 236: INFO = -7
! 237: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 238: INFO = -9
! 239: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN
! 240: INFO = -11
! 241: END IF
! 242: IF( INFO.NE.0 ) THEN
! 243: CALL XERBLA( 'ZTRSYL3', -INFO )
! 244: RETURN
! 245: ELSE IF( LQUERY ) THEN
! 246: RETURN
! 247: END IF
! 248: *
! 249: * Quick return if possible
! 250: *
! 251: SCALE = ONE
! 252: IF( M.EQ.0 .OR. N.EQ.0 )
! 253: $ RETURN
! 254: *
! 255: * Use unblocked code for small problems or if insufficient
! 256: * workspace is provided
! 257: *
! 258: IF( MIN( NBA, NBB ).EQ.1 .OR. LDSWORK.LT.MAX( NBA, NBB ) ) THEN
! 259: CALL ZTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB,
! 260: $ C, LDC, SCALE, INFO )
! 261: RETURN
! 262: END IF
! 263: *
! 264: * Set constants to control overflow
! 265: *
! 266: SMLNUM = DLAMCH( 'S' )
! 267: BIGNUM = ONE / SMLNUM
! 268: *
! 269: * Set local scaling factors.
! 270: *
! 271: DO L = 1, NBB
! 272: DO K = 1, NBA
! 273: SWORK( K, L ) = ONE
! 274: END DO
! 275: END DO
! 276: *
! 277: * Fallback scaling factor to prevent flushing of SWORK( K, L ) to zero.
! 278: * This scaling is to ensure compatibility with TRSYL and may get flushed.
! 279: *
! 280: BUF = ONE
! 281: *
! 282: * Compute upper bounds of blocks of A and B
! 283: *
! 284: AWRK = NBB
! 285: DO K = 1, NBA
! 286: K1 = (K - 1) * NB + 1
! 287: K2 = MIN( K * NB, M ) + 1
! 288: DO L = K, NBA
! 289: L1 = (L - 1) * NB + 1
! 290: L2 = MIN( L * NB, M ) + 1
! 291: IF( NOTRNA ) THEN
! 292: SWORK( K, AWRK + L ) = ZLANGE( 'I', K2-K1, L2-L1,
! 293: $ A( K1, L1 ), LDA, WNRM )
! 294: ELSE
! 295: SWORK( L, AWRK + K ) = ZLANGE( '1', K2-K1, L2-L1,
! 296: $ A( K1, L1 ), LDA, WNRM )
! 297: END IF
! 298: END DO
! 299: END DO
! 300: BWRK = NBB + NBA
! 301: DO K = 1, NBB
! 302: K1 = (K - 1) * NB + 1
! 303: K2 = MIN( K * NB, N ) + 1
! 304: DO L = K, NBB
! 305: L1 = (L - 1) * NB + 1
! 306: L2 = MIN( L * NB, N ) + 1
! 307: IF( NOTRNB ) THEN
! 308: SWORK( K, BWRK + L ) = ZLANGE( 'I', K2-K1, L2-L1,
! 309: $ B( K1, L1 ), LDB, WNRM )
! 310: ELSE
! 311: SWORK( L, BWRK + K ) = ZLANGE( '1', K2-K1, L2-L1,
! 312: $ B( K1, L1 ), LDB, WNRM )
! 313: END IF
! 314: END DO
! 315: END DO
! 316: *
! 317: SGN = DBLE( ISGN )
! 318: CSGN = DCMPLX( SGN, ZERO )
! 319: *
! 320: IF( NOTRNA .AND. NOTRNB ) THEN
! 321: *
! 322: * Solve A*X + ISGN*X*B = scale*C.
! 323: *
! 324: * The (K,L)th block of X is determined starting from
! 325: * bottom-left corner column by column by
! 326: *
! 327: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
! 328: *
! 329: * Where
! 330: * M L-1
! 331: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)].
! 332: * I=K+1 J=1
! 333: *
! 334: * Start loop over block rows (index = K) and block columns (index = L)
! 335: *
! 336: DO K = NBA, 1, -1
! 337: *
! 338: * K1: row index of the first row in X( K, L )
! 339: * K2: row index of the first row in X( K+1, L )
! 340: * so the K2 - K1 is the column count of the block X( K, L )
! 341: *
! 342: K1 = (K - 1) * NB + 1
! 343: K2 = MIN( K * NB, M ) + 1
! 344: DO L = 1, NBB
! 345: *
! 346: * L1: column index of the first column in X( K, L )
! 347: * L2: column index of the first column in X( K, L + 1)
! 348: * so that L2 - L1 is the row count of the block X( K, L )
! 349: *
! 350: L1 = (L - 1) * NB + 1
! 351: L2 = MIN( L * NB, N ) + 1
! 352: *
! 353: CALL ZTRSYL( TRANA, TRANB, ISGN, K2-K1, L2-L1,
! 354: $ A( K1, K1 ), LDA,
! 355: $ B( L1, L1 ), LDB,
! 356: $ C( K1, L1 ), LDC, SCALOC, IINFO )
! 357: INFO = MAX( INFO, IINFO )
! 358: *
! 359: IF( SCALOC * SWORK( K, L ) .EQ. ZERO ) THEN
! 360: IF( SCALOC .EQ. ZERO ) THEN
! 361: * The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
! 362: * is larger than the product of BIGNUM**2 and cannot be
! 363: * represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
! 364: * Mark the computation as pointless.
! 365: BUF = ZERO
! 366: ELSE
! 367: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 368: END IF
! 369: DO JJ = 1, NBB
! 370: DO LL = 1, NBA
! 371: * Bound by BIGNUM to not introduce Inf. The value
! 372: * is irrelevant; corresponding entries of the
! 373: * solution will be flushed in consistency scaling.
! 374: SWORK( LL, JJ ) = MIN( BIGNUM,
! 375: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 376: END DO
! 377: END DO
! 378: END IF
! 379: SWORK( K, L ) = SCALOC * SWORK( K, L )
! 380: XNRM = ZLANGE( 'I', K2-K1, L2-L1, C( K1, L1 ), LDC,
! 381: $ WNRM )
! 382: *
! 383: DO I = K - 1, 1, -1
! 384: *
! 385: * C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
! 386: *
! 387: I1 = (I - 1) * NB + 1
! 388: I2 = MIN( I * NB, M ) + 1
! 389: *
! 390: * Compute scaling factor to survive the linear update
! 391: * simulating consistent scaling.
! 392: *
! 393: CNRM = ZLANGE( 'I', I2-I1, L2-L1, C( I1, L1 ),
! 394: $ LDC, WNRM )
! 395: SCAMIN = MIN( SWORK( I, L ), SWORK( K, L ) )
! 396: CNRM = CNRM * ( SCAMIN / SWORK( I, L ) )
! 397: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 398: ANRM = SWORK( I, AWRK + K )
! 399: SCALOC = DLARMM( ANRM, XNRM, CNRM )
! 400: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 401: * Use second scaling factor to prevent flushing to zero.
! 402: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 403: DO JJ = 1, NBB
! 404: DO LL = 1, NBA
! 405: SWORK( LL, JJ ) = MIN( BIGNUM,
! 406: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 407: END DO
! 408: END DO
! 409: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 410: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 411: END IF
! 412: CNRM = CNRM * SCALOC
! 413: XNRM = XNRM * SCALOC
! 414: *
! 415: * Simultaneously apply the robust update factor and the
! 416: * consistency scaling factor to C( I, L ) and C( K, L ).
! 417: *
! 418: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 419: IF( SCAL .NE. ONE ) THEN
! 420: DO JJ = L1, L2-1
! 421: CALL ZDSCAL( K2-K1, SCAL, C( K1, JJ ), 1)
! 422: END DO
! 423: ENDIF
! 424: *
! 425: SCAL = ( SCAMIN / SWORK( I, L ) ) * SCALOC
! 426: IF( SCAL .NE. ONE ) THEN
! 427: DO LL = L1, L2-1
! 428: CALL ZDSCAL( I2-I1, SCAL, C( I1, LL ), 1)
! 429: END DO
! 430: ENDIF
! 431: *
! 432: * Record current scaling factor
! 433: *
! 434: SWORK( K, L ) = SCAMIN * SCALOC
! 435: SWORK( I, L ) = SCAMIN * SCALOC
! 436: *
! 437: CALL ZGEMM( 'N', 'N', I2-I1, L2-L1, K2-K1, -CONE,
! 438: $ A( I1, K1 ), LDA, C( K1, L1 ), LDC,
! 439: $ CONE, C( I1, L1 ), LDC )
! 440: *
! 441: END DO
! 442: *
! 443: DO J = L + 1, NBB
! 444: *
! 445: * C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
! 446: *
! 447: J1 = (J - 1) * NB + 1
! 448: J2 = MIN( J * NB, N ) + 1
! 449: *
! 450: * Compute scaling factor to survive the linear update
! 451: * simulating consistent scaling.
! 452: *
! 453: CNRM = ZLANGE( 'I', K2-K1, J2-J1, C( K1, J1 ),
! 454: $ LDC, WNRM )
! 455: SCAMIN = MIN( SWORK( K, J ), SWORK( K, L ) )
! 456: CNRM = CNRM * ( SCAMIN / SWORK( K, J ) )
! 457: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 458: BNRM = SWORK(L, BWRK + J)
! 459: SCALOC = DLARMM( BNRM, XNRM, CNRM )
! 460: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 461: * Use second scaling factor to prevent flushing to zero.
! 462: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 463: DO JJ = 1, NBB
! 464: DO LL = 1, NBA
! 465: SWORK( LL, JJ ) = MIN( BIGNUM,
! 466: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 467: END DO
! 468: END DO
! 469: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 470: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 471: END IF
! 472: CNRM = CNRM * SCALOC
! 473: XNRM = XNRM * SCALOC
! 474: *
! 475: * Simultaneously apply the robust update factor and the
! 476: * consistency scaling factor to C( K, J ) and C( K, L).
! 477: *
! 478: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 479: IF( SCAL .NE. ONE ) THEN
! 480: DO LL = L1, L2-1
! 481: CALL ZDSCAL( K2-K1, SCAL, C( K1, LL ), 1 )
! 482: END DO
! 483: ENDIF
! 484: *
! 485: SCAL = ( SCAMIN / SWORK( K, J ) ) * SCALOC
! 486: IF( SCAL .NE. ONE ) THEN
! 487: DO JJ = J1, J2-1
! 488: CALL ZDSCAL( K2-K1, SCAL, C( K1, JJ ), 1 )
! 489: END DO
! 490: ENDIF
! 491: *
! 492: * Record current scaling factor
! 493: *
! 494: SWORK( K, L ) = SCAMIN * SCALOC
! 495: SWORK( K, J ) = SCAMIN * SCALOC
! 496: *
! 497: CALL ZGEMM( 'N', 'N', K2-K1, J2-J1, L2-L1, -CSGN,
! 498: $ C( K1, L1 ), LDC, B( L1, J1 ), LDB,
! 499: $ CONE, C( K1, J1 ), LDC )
! 500: END DO
! 501: END DO
! 502: END DO
! 503: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN
! 504: *
! 505: * Solve A**H *X + ISGN*X*B = scale*C.
! 506: *
! 507: * The (K,L)th block of X is determined starting from
! 508: * upper-left corner column by column by
! 509: *
! 510: * A(K,K)**H*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L)
! 511: *
! 512: * Where
! 513: * K-1 L-1
! 514: * R(K,L) = SUM [A(I,K)**H*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)]
! 515: * I=1 J=1
! 516: *
! 517: * Start loop over block rows (index = K) and block columns (index = L)
! 518: *
! 519: DO K = 1, NBA
! 520: *
! 521: * K1: row index of the first row in X( K, L )
! 522: * K2: row index of the first row in X( K+1, L )
! 523: * so the K2 - K1 is the column count of the block X( K, L )
! 524: *
! 525: K1 = (K - 1) * NB + 1
! 526: K2 = MIN( K * NB, M ) + 1
! 527: DO L = 1, NBB
! 528: *
! 529: * L1: column index of the first column in X( K, L )
! 530: * L2: column index of the first column in X( K, L + 1)
! 531: * so that L2 - L1 is the row count of the block X( K, L )
! 532: *
! 533: L1 = (L - 1) * NB + 1
! 534: L2 = MIN( L * NB, N ) + 1
! 535: *
! 536: CALL ZTRSYL( TRANA, TRANB, ISGN, K2-K1, L2-L1,
! 537: $ A( K1, K1 ), LDA,
! 538: $ B( L1, L1 ), LDB,
! 539: $ C( K1, L1 ), LDC, SCALOC, IINFO )
! 540: INFO = MAX( INFO, IINFO )
! 541: *
! 542: IF( SCALOC * SWORK( K, L ) .EQ. ZERO ) THEN
! 543: IF( SCALOC .EQ. ZERO ) THEN
! 544: * The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
! 545: * is larger than the product of BIGNUM**2 and cannot be
! 546: * represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
! 547: * Mark the computation as pointless.
! 548: BUF = ZERO
! 549: ELSE
! 550: * Use second scaling factor to prevent flushing to zero.
! 551: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 552: END IF
! 553: DO JJ = 1, NBB
! 554: DO LL = 1, NBA
! 555: * Bound by BIGNUM to not introduce Inf. The value
! 556: * is irrelevant; corresponding entries of the
! 557: * solution will be flushed in consistency scaling.
! 558: SWORK( LL, JJ ) = MIN( BIGNUM,
! 559: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 560: END DO
! 561: END DO
! 562: END IF
! 563: SWORK( K, L ) = SCALOC * SWORK( K, L )
! 564: XNRM = ZLANGE( 'I', K2-K1, L2-L1, C( K1, L1 ), LDC,
! 565: $ WNRM )
! 566: *
! 567: DO I = K + 1, NBA
! 568: *
! 569: * C( I, L ) := C( I, L ) - A( K, I )**H * C( K, L )
! 570: *
! 571: I1 = (I - 1) * NB + 1
! 572: I2 = MIN( I * NB, M ) + 1
! 573: *
! 574: * Compute scaling factor to survive the linear update
! 575: * simulating consistent scaling.
! 576: *
! 577: CNRM = ZLANGE( 'I', I2-I1, L2-L1, C( I1, L1 ),
! 578: $ LDC, WNRM )
! 579: SCAMIN = MIN( SWORK( I, L ), SWORK( K, L ) )
! 580: CNRM = CNRM * ( SCAMIN / SWORK( I, L ) )
! 581: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 582: ANRM = SWORK( I, AWRK + K )
! 583: SCALOC = DLARMM( ANRM, XNRM, CNRM )
! 584: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 585: * Use second scaling factor to prevent flushing to zero.
! 586: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 587: DO JJ = 1, NBB
! 588: DO LL = 1, NBA
! 589: SWORK( LL, JJ ) = MIN( BIGNUM,
! 590: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 591: END DO
! 592: END DO
! 593: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 594: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 595: END IF
! 596: CNRM = CNRM * SCALOC
! 597: XNRM = XNRM * SCALOC
! 598: *
! 599: * Simultaneously apply the robust update factor and the
! 600: * consistency scaling factor to to C( I, L ) and C( K, L).
! 601: *
! 602: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 603: IF( SCAL .NE. ONE ) THEN
! 604: DO LL = L1, L2-1
! 605: CALL ZDSCAL( K2-K1, SCAL, C( K1, LL ), 1 )
! 606: END DO
! 607: ENDIF
! 608: *
! 609: SCAL = ( SCAMIN / SWORK( I, L ) ) * SCALOC
! 610: IF( SCAL .NE. ONE ) THEN
! 611: DO LL = L1, L2-1
! 612: CALL ZDSCAL( I2-I1, SCAL, C( I1, LL ), 1 )
! 613: END DO
! 614: ENDIF
! 615: *
! 616: * Record current scaling factor
! 617: *
! 618: SWORK( K, L ) = SCAMIN * SCALOC
! 619: SWORK( I, L ) = SCAMIN * SCALOC
! 620: *
! 621: CALL ZGEMM( 'C', 'N', I2-I1, L2-L1, K2-K1, -CONE,
! 622: $ A( K1, I1 ), LDA, C( K1, L1 ), LDC,
! 623: $ CONE, C( I1, L1 ), LDC )
! 624: END DO
! 625: *
! 626: DO J = L + 1, NBB
! 627: *
! 628: * C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( L, J )
! 629: *
! 630: J1 = (J - 1) * NB + 1
! 631: J2 = MIN( J * NB, N ) + 1
! 632: *
! 633: * Compute scaling factor to survive the linear update
! 634: * simulating consistent scaling.
! 635: *
! 636: CNRM = ZLANGE( 'I', K2-K1, J2-J1, C( K1, J1 ),
! 637: $ LDC, WNRM )
! 638: SCAMIN = MIN( SWORK( K, J ), SWORK( K, L ) )
! 639: CNRM = CNRM * ( SCAMIN / SWORK( K, J ) )
! 640: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 641: BNRM = SWORK( L, BWRK + J )
! 642: SCALOC = DLARMM( BNRM, XNRM, CNRM )
! 643: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 644: * Use second scaling factor to prevent flushing to zero.
! 645: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 646: DO JJ = 1, NBB
! 647: DO LL = 1, NBA
! 648: SWORK( LL, JJ ) = MIN( BIGNUM,
! 649: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 650: END DO
! 651: END DO
! 652: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 653: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 654: END IF
! 655: CNRM = CNRM * SCALOC
! 656: XNRM = XNRM * SCALOC
! 657: *
! 658: * Simultaneously apply the robust update factor and the
! 659: * consistency scaling factor to to C( K, J ) and C( K, L).
! 660: *
! 661: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 662: IF( SCAL .NE. ONE ) THEN
! 663: DO LL = L1, L2-1
! 664: CALL ZDSCAL( K2-K1, SCAL, C( K1, LL ), 1 )
! 665: END DO
! 666: ENDIF
! 667: *
! 668: SCAL = ( SCAMIN / SWORK( K, J ) ) * SCALOC
! 669: IF( SCAL .NE. ONE ) THEN
! 670: DO JJ = J1, J2-1
! 671: CALL ZDSCAL( K2-K1, SCAL, C( K1, JJ ), 1 )
! 672: END DO
! 673: ENDIF
! 674: *
! 675: * Record current scaling factor
! 676: *
! 677: SWORK( K, L ) = SCAMIN * SCALOC
! 678: SWORK( K, J ) = SCAMIN * SCALOC
! 679: *
! 680: CALL ZGEMM( 'N', 'N', K2-K1, J2-J1, L2-L1, -CSGN,
! 681: $ C( K1, L1 ), LDC, B( L1, J1 ), LDB,
! 682: $ CONE, C( K1, J1 ), LDC )
! 683: END DO
! 684: END DO
! 685: END DO
! 686: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN
! 687: *
! 688: * Solve A**H *X + ISGN*X*B**H = scale*C.
! 689: *
! 690: * The (K,L)th block of X is determined starting from
! 691: * top-right corner column by column by
! 692: *
! 693: * A(K,K)**H*X(K,L) + ISGN*X(K,L)*B(L,L)**H = C(K,L) - R(K,L)
! 694: *
! 695: * Where
! 696: * K-1 N
! 697: * R(K,L) = SUM [A(I,K)**H*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**H].
! 698: * I=1 J=L+1
! 699: *
! 700: * Start loop over block rows (index = K) and block columns (index = L)
! 701: *
! 702: DO K = 1, NBA
! 703: *
! 704: * K1: row index of the first row in X( K, L )
! 705: * K2: row index of the first row in X( K+1, L )
! 706: * so the K2 - K1 is the column count of the block X( K, L )
! 707: *
! 708: K1 = (K - 1) * NB + 1
! 709: K2 = MIN( K * NB, M ) + 1
! 710: DO L = NBB, 1, -1
! 711: *
! 712: * L1: column index of the first column in X( K, L )
! 713: * L2: column index of the first column in X( K, L + 1)
! 714: * so that L2 - L1 is the row count of the block X( K, L )
! 715: *
! 716: L1 = (L - 1) * NB + 1
! 717: L2 = MIN( L * NB, N ) + 1
! 718: *
! 719: CALL ZTRSYL( TRANA, TRANB, ISGN, K2-K1, L2-L1,
! 720: $ A( K1, K1 ), LDA,
! 721: $ B( L1, L1 ), LDB,
! 722: $ C( K1, L1 ), LDC, SCALOC, IINFO )
! 723: INFO = MAX( INFO, IINFO )
! 724: *
! 725: IF( SCALOC * SWORK( K, L ) .EQ. ZERO ) THEN
! 726: IF( SCALOC .EQ. ZERO ) THEN
! 727: * The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
! 728: * is larger than the product of BIGNUM**2 and cannot be
! 729: * represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
! 730: * Mark the computation as pointless.
! 731: BUF = ZERO
! 732: ELSE
! 733: * Use second scaling factor to prevent flushing to zero.
! 734: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 735: END IF
! 736: DO JJ = 1, NBB
! 737: DO LL = 1, NBA
! 738: * Bound by BIGNUM to not introduce Inf. The value
! 739: * is irrelevant; corresponding entries of the
! 740: * solution will be flushed in consistency scaling.
! 741: SWORK( LL, JJ ) = MIN( BIGNUM,
! 742: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 743: END DO
! 744: END DO
! 745: END IF
! 746: SWORK( K, L ) = SCALOC * SWORK( K, L )
! 747: XNRM = ZLANGE( 'I', K2-K1, L2-L1, C( K1, L1 ), LDC,
! 748: $ WNRM )
! 749: *
! 750: DO I = K + 1, NBA
! 751: *
! 752: * C( I, L ) := C( I, L ) - A( K, I )**H * C( K, L )
! 753: *
! 754: I1 = (I - 1) * NB + 1
! 755: I2 = MIN( I * NB, M ) + 1
! 756: *
! 757: * Compute scaling factor to survive the linear update
! 758: * simulating consistent scaling.
! 759: *
! 760: CNRM = ZLANGE( 'I', I2-I1, L2-L1, C( I1, L1 ),
! 761: $ LDC, WNRM )
! 762: SCAMIN = MIN( SWORK( I, L ), SWORK( K, L ) )
! 763: CNRM = CNRM * ( SCAMIN / SWORK( I, L ) )
! 764: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 765: ANRM = SWORK( I, AWRK + K )
! 766: SCALOC = DLARMM( ANRM, XNRM, CNRM )
! 767: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 768: * Use second scaling factor to prevent flushing to zero.
! 769: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 770: DO JJ = 1, NBB
! 771: DO LL = 1, NBA
! 772: SWORK( LL, JJ ) = MIN( BIGNUM,
! 773: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 774: END DO
! 775: END DO
! 776: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 777: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 778: END IF
! 779: CNRM = CNRM * SCALOC
! 780: XNRM = XNRM * SCALOC
! 781: *
! 782: * Simultaneously apply the robust update factor and the
! 783: * consistency scaling factor to C( I, L ) and C( K, L).
! 784: *
! 785: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 786: IF( SCAL .NE. ONE ) THEN
! 787: DO LL = L1, L2-1
! 788: CALL ZDSCAL( K2-K1, SCAL, C( K1, LL ), 1 )
! 789: END DO
! 790: ENDIF
! 791: *
! 792: SCAL = ( SCAMIN / SWORK( I, L ) ) * SCALOC
! 793: IF( SCAL .NE. ONE ) THEN
! 794: DO LL = L1, L2-1
! 795: CALL ZDSCAL( I2-I1, SCAL, C( I1, LL ), 1 )
! 796: END DO
! 797: ENDIF
! 798: *
! 799: * Record current scaling factor
! 800: *
! 801: SWORK( K, L ) = SCAMIN * SCALOC
! 802: SWORK( I, L ) = SCAMIN * SCALOC
! 803: *
! 804: CALL ZGEMM( 'C', 'N', I2-I1, L2-L1, K2-K1, -CONE,
! 805: $ A( K1, I1 ), LDA, C( K1, L1 ), LDC,
! 806: $ CONE, C( I1, L1 ), LDC )
! 807: END DO
! 808: *
! 809: DO J = 1, L - 1
! 810: *
! 811: * C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**H
! 812: *
! 813: J1 = (J - 1) * NB + 1
! 814: J2 = MIN( J * NB, N ) + 1
! 815: *
! 816: * Compute scaling factor to survive the linear update
! 817: * simulating consistent scaling.
! 818: *
! 819: CNRM = ZLANGE( 'I', K2-K1, J2-J1, C( K1, J1 ),
! 820: $ LDC, WNRM )
! 821: SCAMIN = MIN( SWORK( K, J ), SWORK( K, L ) )
! 822: CNRM = CNRM * ( SCAMIN / SWORK( K, J ) )
! 823: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 824: BNRM = SWORK( L, BWRK + J )
! 825: SCALOC = DLARMM( BNRM, XNRM, CNRM )
! 826: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 827: * Use second scaling factor to prevent flushing to zero.
! 828: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 829: DO JJ = 1, NBB
! 830: DO LL = 1, NBA
! 831: SWORK( LL, JJ ) = MIN( BIGNUM,
! 832: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 833: END DO
! 834: END DO
! 835: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 836: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 837: END IF
! 838: CNRM = CNRM * SCALOC
! 839: XNRM = XNRM * SCALOC
! 840: *
! 841: * Simultaneously apply the robust update factor and the
! 842: * consistency scaling factor to C( K, J ) and C( K, L).
! 843: *
! 844: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 845: IF( SCAL .NE. ONE ) THEN
! 846: DO LL = L1, L2-1
! 847: CALL ZDSCAL( K2-K1, SCAL, C( K1, LL ), 1)
! 848: END DO
! 849: ENDIF
! 850: *
! 851: SCAL = ( SCAMIN / SWORK( K, J ) ) * SCALOC
! 852: IF( SCAL .NE. ONE ) THEN
! 853: DO JJ = J1, J2-1
! 854: CALL ZDSCAL( K2-K1, SCAL, C( K1, JJ ), 1 )
! 855: END DO
! 856: ENDIF
! 857: *
! 858: * Record current scaling factor
! 859: *
! 860: SWORK( K, L ) = SCAMIN * SCALOC
! 861: SWORK( K, J ) = SCAMIN * SCALOC
! 862: *
! 863: CALL ZGEMM( 'N', 'C', K2-K1, J2-J1, L2-L1, -CSGN,
! 864: $ C( K1, L1 ), LDC, B( J1, L1 ), LDB,
! 865: $ CONE, C( K1, J1 ), LDC )
! 866: END DO
! 867: END DO
! 868: END DO
! 869: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN
! 870: *
! 871: * Solve A*X + ISGN*X*B**H = scale*C.
! 872: *
! 873: * The (K,L)th block of X is determined starting from
! 874: * bottom-right corner column by column by
! 875: *
! 876: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)**H = C(K,L) - R(K,L)
! 877: *
! 878: * Where
! 879: * M N
! 880: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)**H].
! 881: * I=K+1 J=L+1
! 882: *
! 883: * Start loop over block rows (index = K) and block columns (index = L)
! 884: *
! 885: DO K = NBA, 1, -1
! 886: *
! 887: * K1: row index of the first row in X( K, L )
! 888: * K2: row index of the first row in X( K+1, L )
! 889: * so the K2 - K1 is the column count of the block X( K, L )
! 890: *
! 891: K1 = (K - 1) * NB + 1
! 892: K2 = MIN( K * NB, M ) + 1
! 893: DO L = NBB, 1, -1
! 894: *
! 895: * L1: column index of the first column in X( K, L )
! 896: * L2: column index of the first column in X( K, L + 1)
! 897: * so that L2 - L1 is the row count of the block X( K, L )
! 898: *
! 899: L1 = (L - 1) * NB + 1
! 900: L2 = MIN( L * NB, N ) + 1
! 901: *
! 902: CALL ZTRSYL( TRANA, TRANB, ISGN, K2-K1, L2-L1,
! 903: $ A( K1, K1 ), LDA,
! 904: $ B( L1, L1 ), LDB,
! 905: $ C( K1, L1 ), LDC, SCALOC, IINFO )
! 906: INFO = MAX( INFO, IINFO )
! 907: *
! 908: IF( SCALOC * SWORK( K, L ) .EQ. ZERO ) THEN
! 909: IF( SCALOC .EQ. ZERO ) THEN
! 910: * The magnitude of the largest entry of X(K1:K2-1, L1:L2-1)
! 911: * is larger than the product of BIGNUM**2 and cannot be
! 912: * represented in the form (1/SCALE)*X(K1:K2-1, L1:L2-1).
! 913: * Mark the computation as pointless.
! 914: BUF = ZERO
! 915: ELSE
! 916: * Use second scaling factor to prevent flushing to zero.
! 917: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 918: END IF
! 919: DO JJ = 1, NBB
! 920: DO LL = 1, NBA
! 921: * Bound by BIGNUM to not introduce Inf. The value
! 922: * is irrelevant; corresponding entries of the
! 923: * solution will be flushed in consistency scaling.
! 924: SWORK( LL, JJ ) = MIN( BIGNUM,
! 925: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 926: END DO
! 927: END DO
! 928: END IF
! 929: SWORK( K, L ) = SCALOC * SWORK( K, L )
! 930: XNRM = ZLANGE( 'I', K2-K1, L2-L1, C( K1, L1 ), LDC,
! 931: $ WNRM )
! 932: *
! 933: DO I = 1, K - 1
! 934: *
! 935: * C( I, L ) := C( I, L ) - A( I, K ) * C( K, L )
! 936: *
! 937: I1 = (I - 1) * NB + 1
! 938: I2 = MIN( I * NB, M ) + 1
! 939: *
! 940: * Compute scaling factor to survive the linear update
! 941: * simulating consistent scaling.
! 942: *
! 943: CNRM = ZLANGE( 'I', I2-I1, L2-L1, C( I1, L1 ),
! 944: $ LDC, WNRM )
! 945: SCAMIN = MIN( SWORK( I, L ), SWORK( K, L ) )
! 946: CNRM = CNRM * ( SCAMIN / SWORK( I, L ) )
! 947: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 948: ANRM = SWORK( I, AWRK + K )
! 949: SCALOC = DLARMM( ANRM, XNRM, CNRM )
! 950: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 951: * Use second scaling factor to prevent flushing to zero.
! 952: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 953: DO JJ = 1, NBB
! 954: DO LL = 1, NBA
! 955: SWORK( LL, JJ ) = MIN( BIGNUM,
! 956: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 957: END DO
! 958: END DO
! 959: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 960: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 961: END IF
! 962: CNRM = CNRM * SCALOC
! 963: XNRM = XNRM * SCALOC
! 964: *
! 965: * Simultaneously apply the robust update factor and the
! 966: * consistency scaling factor to C( I, L ) and C( K, L).
! 967: *
! 968: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 969: IF( SCAL .NE. ONE ) THEN
! 970: DO LL = L1, L2-1
! 971: CALL ZDSCAL( K2-K1, SCAL, C( K1, LL ), 1 )
! 972: END DO
! 973: ENDIF
! 974: *
! 975: SCAL = ( SCAMIN / SWORK( I, L ) ) * SCALOC
! 976: IF( SCAL .NE. ONE ) THEN
! 977: DO LL = L1, L2-1
! 978: CALL ZDSCAL( I2-I1, SCAL, C( I1, LL ), 1 )
! 979: END DO
! 980: ENDIF
! 981: *
! 982: * Record current scaling factor
! 983: *
! 984: SWORK( K, L ) = SCAMIN * SCALOC
! 985: SWORK( I, L ) = SCAMIN * SCALOC
! 986: *
! 987: CALL ZGEMM( 'N', 'N', I2-I1, L2-L1, K2-K1, -CONE,
! 988: $ A( I1, K1 ), LDA, C( K1, L1 ), LDC,
! 989: $ CONE, C( I1, L1 ), LDC )
! 990: *
! 991: END DO
! 992: *
! 993: DO J = 1, L - 1
! 994: *
! 995: * C( K, J ) := C( K, J ) - SGN * C( K, L ) * B( J, L )**H
! 996: *
! 997: J1 = (J - 1) * NB + 1
! 998: J2 = MIN( J * NB, N ) + 1
! 999: *
! 1000: * Compute scaling factor to survive the linear update
! 1001: * simulating consistent scaling.
! 1002: *
! 1003: CNRM = ZLANGE( 'I', K2-K1, J2-J1, C( K1, J1 ),
! 1004: $ LDC, WNRM )
! 1005: SCAMIN = MIN( SWORK( K, J ), SWORK( K, L ) )
! 1006: CNRM = CNRM * ( SCAMIN / SWORK( K, J ) )
! 1007: XNRM = XNRM * ( SCAMIN / SWORK( K, L ) )
! 1008: BNRM = SWORK( L, BWRK + J )
! 1009: SCALOC = DLARMM( BNRM, XNRM, CNRM )
! 1010: IF( SCALOC * SCAMIN .EQ. ZERO ) THEN
! 1011: * Use second scaling factor to prevent flushing to zero.
! 1012: BUF = BUF*2.D0**EXPONENT( SCALOC )
! 1013: DO JJ = 1, NBB
! 1014: DO LL = 1, NBA
! 1015: SWORK( LL, JJ ) = MIN( BIGNUM,
! 1016: $ SWORK( LL, JJ ) / 2.D0**EXPONENT( SCALOC ) )
! 1017: END DO
! 1018: END DO
! 1019: SCAMIN = SCAMIN / 2.D0**EXPONENT( SCALOC )
! 1020: SCALOC = SCALOC / 2.D0**EXPONENT( SCALOC )
! 1021: END IF
! 1022: CNRM = CNRM * SCALOC
! 1023: XNRM = XNRM * SCALOC
! 1024: *
! 1025: * Simultaneously apply the robust update factor and the
! 1026: * consistency scaling factor to C( K, J ) and C( K, L).
! 1027: *
! 1028: SCAL = ( SCAMIN / SWORK( K, L ) ) * SCALOC
! 1029: IF( SCAL .NE. ONE ) THEN
! 1030: DO JJ = L1, L2-1
! 1031: CALL ZDSCAL( K2-K1, SCAL, C( K1, JJ ), 1 )
! 1032: END DO
! 1033: ENDIF
! 1034: *
! 1035: SCAL = ( SCAMIN / SWORK( K, J ) ) * SCALOC
! 1036: IF( SCAL .NE. ONE ) THEN
! 1037: DO JJ = J1, J2-1
! 1038: CALL ZDSCAL( K2-K1, SCAL, C( K1, JJ ), 1 )
! 1039: END DO
! 1040: ENDIF
! 1041: *
! 1042: * Record current scaling factor
! 1043: *
! 1044: SWORK( K, L ) = SCAMIN * SCALOC
! 1045: SWORK( K, J ) = SCAMIN * SCALOC
! 1046: *
! 1047: CALL ZGEMM( 'N', 'C', K2-K1, J2-J1, L2-L1, -CSGN,
! 1048: $ C( K1, L1 ), LDC, B( J1, L1 ), LDB,
! 1049: $ CONE, C( K1, J1 ), LDC )
! 1050: END DO
! 1051: END DO
! 1052: END DO
! 1053: *
! 1054: END IF
! 1055: *
! 1056: * Reduce local scaling factors
! 1057: *
! 1058: SCALE = SWORK( 1, 1 )
! 1059: DO K = 1, NBA
! 1060: DO L = 1, NBB
! 1061: SCALE = MIN( SCALE, SWORK( K, L ) )
! 1062: END DO
! 1063: END DO
! 1064: IF( SCALE .EQ. ZERO ) THEN
! 1065: *
! 1066: * The magnitude of the largest entry of the solution is larger
! 1067: * than the product of BIGNUM**2 and cannot be represented in the
! 1068: * form (1/SCALE)*X if SCALE is DOUBLE PRECISION. Set SCALE to
! 1069: * zero and give up.
! 1070: *
! 1071: SWORK(1,1) = MAX( NBA, NBB )
! 1072: SWORK(2,1) = 2 * NBB + NBA
! 1073: RETURN
! 1074: END IF
! 1075: *
! 1076: * Realize consistent scaling
! 1077: *
! 1078: DO K = 1, NBA
! 1079: K1 = (K - 1) * NB + 1
! 1080: K2 = MIN( K * NB, M ) + 1
! 1081: DO L = 1, NBB
! 1082: L1 = (L - 1) * NB + 1
! 1083: L2 = MIN( L * NB, N ) + 1
! 1084: SCAL = SCALE / SWORK( K, L )
! 1085: IF( SCAL .NE. ONE ) THEN
! 1086: DO LL = L1, L2-1
! 1087: CALL ZDSCAL( K2-K1, SCAL, C( K1, LL ), 1 )
! 1088: END DO
! 1089: ENDIF
! 1090: END DO
! 1091: END DO
! 1092: *
! 1093: IF( BUF .NE. ONE .AND. BUF.GT.ZERO ) THEN
! 1094: *
! 1095: * Decrease SCALE as much as possible.
! 1096: *
! 1097: SCALOC = MIN( SCALE / SMLNUM, ONE / BUF )
! 1098: BUF = BUF * SCALOC
! 1099: SCALE = SCALE / SCALOC
! 1100: END IF
! 1101: *
! 1102: IF( BUF.NE.ONE .AND. BUF.GT.ZERO ) THEN
! 1103: *
! 1104: * In case of overly aggressive scaling during the computation,
! 1105: * flushing of the global scale factor may be prevented by
! 1106: * undoing some of the scaling. This step is to ensure that
! 1107: * this routine flushes only scale factors that TRSYL also
! 1108: * flushes and be usable as a drop-in replacement.
! 1109: *
! 1110: * How much can the normwise largest entry be upscaled?
! 1111: *
! 1112: SCAL = MAX( ABS( DBLE( C( 1, 1 ) ) ),
! 1113: $ ABS( DIMAG( C ( 1, 1 ) ) ) )
! 1114: DO K = 1, M
! 1115: DO L = 1, N
! 1116: SCAL = MAX( SCAL, ABS( DBLE ( C( K, L ) ) ),
! 1117: $ ABS( DIMAG ( C( K, L ) ) ) )
! 1118: END DO
! 1119: END DO
! 1120: *
! 1121: * Increase BUF as close to 1 as possible and apply scaling.
! 1122: *
! 1123: SCALOC = MIN( BIGNUM / SCAL, ONE / BUF )
! 1124: BUF = BUF * SCALOC
! 1125: CALL ZLASCL( 'G', -1, -1, ONE, SCALOC, M, N, C, LDC, IINFO )
! 1126: END IF
! 1127: *
! 1128: * Combine with buffer scaling factor. SCALE will be flushed if
! 1129: * BUF is less than one here.
! 1130: *
! 1131: SCALE = SCALE * BUF
! 1132: *
! 1133: * Restore workspace dimensions
! 1134: *
! 1135: SWORK(1,1) = MAX( NBA, NBB )
! 1136: SWORK(2,1) = 2 * NBB + NBA
! 1137: *
! 1138: RETURN
! 1139: *
! 1140: * End of ZTRSYL3
! 1141: *
! 1142: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to ztrsyl3.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack