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