![]() ![]() | ![]() |
Mise à jour de lapack vers la version 3.3.0.
1: SUBROUTINE DTRSYL( TRANA, TRANB, ISGN, M, N, A, LDA, B, LDB, C, 2: $ LDC, SCALE, INFO ) 3: * 4: * -- LAPACK routine (version 3.2) -- 5: * -- LAPACK is a software package provided by Univ. of Tennessee, -- 6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 7: * November 2006 8: * 9: * .. Scalar Arguments .. 10: CHARACTER TRANA, TRANB 11: INTEGER INFO, ISGN, LDA, LDB, LDC, M, N 12: DOUBLE PRECISION SCALE 13: * .. 14: * .. Array Arguments .. 15: DOUBLE PRECISION A( LDA, * ), B( LDB, * ), C( LDC, * ) 16: * .. 17: * 18: * Purpose 19: * ======= 20: * 21: * DTRSYL solves the real Sylvester matrix equation: 22: * 23: * op(A)*X + X*op(B) = scale*C or 24: * op(A)*X - X*op(B) = scale*C, 25: * 26: * where op(A) = A or A**T, and A and B are both upper quasi- 27: * triangular. A is M-by-M and B is N-by-N; the right hand side C and 28: * the solution X are M-by-N; and scale is an output scale factor, set 29: * <= 1 to avoid overflow in X. 30: * 31: * A and B must be in Schur canonical form (as returned by DHSEQR), that 32: * is, block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; 33: * each 2-by-2 diagonal block has its diagonal elements equal and its 34: * off-diagonal elements of opposite sign. 35: * 36: * Arguments 37: * ========= 38: * 39: * TRANA (input) CHARACTER*1 40: * Specifies the option op(A): 41: * = 'N': op(A) = A (No transpose) 42: * = 'T': op(A) = A**T (Transpose) 43: * = 'C': op(A) = A**H (Conjugate transpose = Transpose) 44: * 45: * TRANB (input) CHARACTER*1 46: * Specifies the option op(B): 47: * = 'N': op(B) = B (No transpose) 48: * = 'T': op(B) = B**T (Transpose) 49: * = 'C': op(B) = B**H (Conjugate transpose = Transpose) 50: * 51: * ISGN (input) INTEGER 52: * Specifies the sign in the equation: 53: * = +1: solve op(A)*X + X*op(B) = scale*C 54: * = -1: solve op(A)*X - X*op(B) = scale*C 55: * 56: * M (input) INTEGER 57: * The order of the matrix A, and the number of rows in the 58: * matrices X and C. M >= 0. 59: * 60: * N (input) INTEGER 61: * The order of the matrix B, and the number of columns in the 62: * matrices X and C. N >= 0. 63: * 64: * A (input) DOUBLE PRECISION array, dimension (LDA,M) 65: * The upper quasi-triangular matrix A, in Schur canonical form. 66: * 67: * LDA (input) INTEGER 68: * The leading dimension of the array A. LDA >= max(1,M). 69: * 70: * B (input) DOUBLE PRECISION array, dimension (LDB,N) 71: * The upper quasi-triangular matrix B, in Schur canonical form. 72: * 73: * LDB (input) INTEGER 74: * The leading dimension of the array B. LDB >= max(1,N). 75: * 76: * C (input/output) DOUBLE PRECISION array, dimension (LDC,N) 77: * On entry, the M-by-N right hand side matrix C. 78: * On exit, C is overwritten by the solution matrix X. 79: * 80: * LDC (input) INTEGER 81: * The leading dimension of the array C. LDC >= max(1,M) 82: * 83: * SCALE (output) DOUBLE PRECISION 84: * The scale factor, scale, set <= 1 to avoid overflow in X. 85: * 86: * INFO (output) INTEGER 87: * = 0: successful exit 88: * < 0: if INFO = -i, the i-th argument had an illegal value 89: * = 1: A and B have common or very close eigenvalues; perturbed 90: * values were used to solve the equation (but the matrices 91: * A and B are unchanged). 92: * 93: * ===================================================================== 94: * 95: * .. Parameters .. 96: DOUBLE PRECISION ZERO, ONE 97: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 98: * .. 99: * .. Local Scalars .. 100: LOGICAL NOTRNA, NOTRNB 101: INTEGER IERR, J, K, K1, K2, KNEXT, L, L1, L2, LNEXT 102: DOUBLE PRECISION A11, BIGNUM, DA11, DB, EPS, SCALOC, SGN, SMIN, 103: $ SMLNUM, SUML, SUMR, XNORM 104: * .. 105: * .. Local Arrays .. 106: DOUBLE PRECISION DUM( 1 ), VEC( 2, 2 ), X( 2, 2 ) 107: * .. 108: * .. External Functions .. 109: LOGICAL LSAME 110: DOUBLE PRECISION DDOT, DLAMCH, DLANGE 111: EXTERNAL LSAME, DDOT, DLAMCH, DLANGE 112: * .. 113: * .. External Subroutines .. 114: EXTERNAL DLABAD, DLALN2, DLASY2, DSCAL, XERBLA 115: * .. 116: * .. Intrinsic Functions .. 117: INTRINSIC ABS, DBLE, MAX, MIN 118: * .. 119: * .. Executable Statements .. 120: * 121: * Decode and Test input parameters 122: * 123: NOTRNA = LSAME( TRANA, 'N' ) 124: NOTRNB = LSAME( TRANB, 'N' ) 125: * 126: INFO = 0 127: IF( .NOT.NOTRNA .AND. .NOT.LSAME( TRANA, 'T' ) .AND. .NOT. 128: $ LSAME( TRANA, 'C' ) ) THEN 129: INFO = -1 130: ELSE IF( .NOT.NOTRNB .AND. .NOT.LSAME( TRANB, 'T' ) .AND. .NOT. 131: $ LSAME( TRANB, 'C' ) ) THEN 132: INFO = -2 133: ELSE IF( ISGN.NE.1 .AND. ISGN.NE.-1 ) THEN 134: INFO = -3 135: ELSE IF( M.LT.0 ) THEN 136: INFO = -4 137: ELSE IF( N.LT.0 ) THEN 138: INFO = -5 139: ELSE IF( LDA.LT.MAX( 1, M ) ) THEN 140: INFO = -7 141: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 142: INFO = -9 143: ELSE IF( LDC.LT.MAX( 1, M ) ) THEN 144: INFO = -11 145: END IF 146: IF( INFO.NE.0 ) THEN 147: CALL XERBLA( 'DTRSYL', -INFO ) 148: RETURN 149: END IF 150: * 151: * Quick return if possible 152: * 153: SCALE = ONE 154: IF( M.EQ.0 .OR. N.EQ.0 ) 155: $ RETURN 156: * 157: * Set constants to control overflow 158: * 159: EPS = DLAMCH( 'P' ) 160: SMLNUM = DLAMCH( 'S' ) 161: BIGNUM = ONE / SMLNUM 162: CALL DLABAD( SMLNUM, BIGNUM ) 163: SMLNUM = SMLNUM*DBLE( M*N ) / EPS 164: BIGNUM = ONE / SMLNUM 165: * 166: SMIN = MAX( SMLNUM, EPS*DLANGE( 'M', M, M, A, LDA, DUM ), 167: $ EPS*DLANGE( 'M', N, N, B, LDB, DUM ) ) 168: * 169: SGN = ISGN 170: * 171: IF( NOTRNA .AND. NOTRNB ) THEN 172: * 173: * Solve A*X + ISGN*X*B = scale*C. 174: * 175: * The (K,L)th block of X is determined starting from 176: * bottom-left corner column by column by 177: * 178: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 179: * 180: * Where 181: * M L-1 182: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(J,L)]. 183: * I=K+1 J=1 184: * 185: * Start column loop (index = L) 186: * L1 (L2) : column index of the first (first) row of X(K,L). 187: * 188: LNEXT = 1 189: DO 60 L = 1, N 190: IF( L.LT.LNEXT ) 191: $ GO TO 60 192: IF( L.EQ.N ) THEN 193: L1 = L 194: L2 = L 195: ELSE 196: IF( B( L+1, L ).NE.ZERO ) THEN 197: L1 = L 198: L2 = L + 1 199: LNEXT = L + 2 200: ELSE 201: L1 = L 202: L2 = L 203: LNEXT = L + 1 204: END IF 205: END IF 206: * 207: * Start row loop (index = K) 208: * K1 (K2): row index of the first (last) row of X(K,L). 209: * 210: KNEXT = M 211: DO 50 K = M, 1, -1 212: IF( K.GT.KNEXT ) 213: $ GO TO 50 214: IF( K.EQ.1 ) THEN 215: K1 = K 216: K2 = K 217: ELSE 218: IF( A( K, K-1 ).NE.ZERO ) THEN 219: K1 = K - 1 220: K2 = K 221: KNEXT = K - 2 222: ELSE 223: K1 = K 224: K2 = K 225: KNEXT = K - 1 226: END IF 227: END IF 228: * 229: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 230: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 231: $ C( MIN( K1+1, M ), L1 ), 1 ) 232: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 233: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 234: SCALOC = ONE 235: * 236: A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 237: DA11 = ABS( A11 ) 238: IF( DA11.LE.SMIN ) THEN 239: A11 = SMIN 240: DA11 = SMIN 241: INFO = 1 242: END IF 243: DB = ABS( VEC( 1, 1 ) ) 244: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 245: IF( DB.GT.BIGNUM*DA11 ) 246: $ SCALOC = ONE / DB 247: END IF 248: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 249: * 250: IF( SCALOC.NE.ONE ) THEN 251: DO 10 J = 1, N 252: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 253: 10 CONTINUE 254: SCALE = SCALE*SCALOC 255: END IF 256: C( K1, L1 ) = X( 1, 1 ) 257: * 258: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 259: * 260: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 261: $ C( MIN( K2+1, M ), L1 ), 1 ) 262: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 263: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 264: * 265: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 266: $ C( MIN( K2+1, M ), L1 ), 1 ) 267: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 268: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 269: * 270: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 271: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 272: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 273: IF( IERR.NE.0 ) 274: $ INFO = 1 275: * 276: IF( SCALOC.NE.ONE ) THEN 277: DO 20 J = 1, N 278: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 279: 20 CONTINUE 280: SCALE = SCALE*SCALOC 281: END IF 282: C( K1, L1 ) = X( 1, 1 ) 283: C( K2, L1 ) = X( 2, 1 ) 284: * 285: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 286: * 287: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 288: $ C( MIN( K1+1, M ), L1 ), 1 ) 289: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 290: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 291: * 292: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 293: $ C( MIN( K1+1, M ), L2 ), 1 ) 294: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 295: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 296: * 297: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 298: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 299: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 300: IF( IERR.NE.0 ) 301: $ INFO = 1 302: * 303: IF( SCALOC.NE.ONE ) THEN 304: DO 30 J = 1, N 305: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 306: 30 CONTINUE 307: SCALE = SCALE*SCALOC 308: END IF 309: C( K1, L1 ) = X( 1, 1 ) 310: C( K1, L2 ) = X( 2, 1 ) 311: * 312: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 313: * 314: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 315: $ C( MIN( K2+1, M ), L1 ), 1 ) 316: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 317: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 318: * 319: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 320: $ C( MIN( K2+1, M ), L2 ), 1 ) 321: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 322: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 323: * 324: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 325: $ C( MIN( K2+1, M ), L1 ), 1 ) 326: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 327: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 328: * 329: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 330: $ C( MIN( K2+1, M ), L2 ), 1 ) 331: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 332: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 333: * 334: CALL DLASY2( .FALSE., .FALSE., ISGN, 2, 2, 335: $ A( K1, K1 ), LDA, B( L1, L1 ), LDB, VEC, 336: $ 2, SCALOC, X, 2, XNORM, IERR ) 337: IF( IERR.NE.0 ) 338: $ INFO = 1 339: * 340: IF( SCALOC.NE.ONE ) THEN 341: DO 40 J = 1, N 342: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 343: 40 CONTINUE 344: SCALE = SCALE*SCALOC 345: END IF 346: C( K1, L1 ) = X( 1, 1 ) 347: C( K1, L2 ) = X( 1, 2 ) 348: C( K2, L1 ) = X( 2, 1 ) 349: C( K2, L2 ) = X( 2, 2 ) 350: END IF 351: * 352: 50 CONTINUE 353: * 354: 60 CONTINUE 355: * 356: ELSE IF( .NOT.NOTRNA .AND. NOTRNB ) THEN 357: * 358: * Solve A' *X + ISGN*X*B = scale*C. 359: * 360: * The (K,L)th block of X is determined starting from 361: * upper-left corner column by column by 362: * 363: * A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L) = C(K,L) - R(K,L) 364: * 365: * Where 366: * K-1 L-1 367: * R(K,L) = SUM [A(I,K)'*X(I,L)] +ISGN*SUM [X(K,J)*B(J,L)] 368: * I=1 J=1 369: * 370: * Start column loop (index = L) 371: * L1 (L2): column index of the first (last) row of X(K,L) 372: * 373: LNEXT = 1 374: DO 120 L = 1, N 375: IF( L.LT.LNEXT ) 376: $ GO TO 120 377: IF( L.EQ.N ) THEN 378: L1 = L 379: L2 = L 380: ELSE 381: IF( B( L+1, L ).NE.ZERO ) THEN 382: L1 = L 383: L2 = L + 1 384: LNEXT = L + 2 385: ELSE 386: L1 = L 387: L2 = L 388: LNEXT = L + 1 389: END IF 390: END IF 391: * 392: * Start row loop (index = K) 393: * K1 (K2): row index of the first (last) row of X(K,L) 394: * 395: KNEXT = 1 396: DO 110 K = 1, M 397: IF( K.LT.KNEXT ) 398: $ GO TO 110 399: IF( K.EQ.M ) THEN 400: K1 = K 401: K2 = K 402: ELSE 403: IF( A( K+1, K ).NE.ZERO ) THEN 404: K1 = K 405: K2 = K + 1 406: KNEXT = K + 2 407: ELSE 408: K1 = K 409: K2 = K 410: KNEXT = K + 1 411: END IF 412: END IF 413: * 414: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 415: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 416: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 417: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 418: SCALOC = ONE 419: * 420: A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 421: DA11 = ABS( A11 ) 422: IF( DA11.LE.SMIN ) THEN 423: A11 = SMIN 424: DA11 = SMIN 425: INFO = 1 426: END IF 427: DB = ABS( VEC( 1, 1 ) ) 428: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 429: IF( DB.GT.BIGNUM*DA11 ) 430: $ SCALOC = ONE / DB 431: END IF 432: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 433: * 434: IF( SCALOC.NE.ONE ) THEN 435: DO 70 J = 1, N 436: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 437: 70 CONTINUE 438: SCALE = SCALE*SCALOC 439: END IF 440: C( K1, L1 ) = X( 1, 1 ) 441: * 442: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 443: * 444: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 445: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 446: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 447: * 448: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 449: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 450: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 451: * 452: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 453: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 454: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 455: IF( IERR.NE.0 ) 456: $ INFO = 1 457: * 458: IF( SCALOC.NE.ONE ) THEN 459: DO 80 J = 1, N 460: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 461: 80 CONTINUE 462: SCALE = SCALE*SCALOC 463: END IF 464: C( K1, L1 ) = X( 1, 1 ) 465: C( K2, L1 ) = X( 2, 1 ) 466: * 467: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 468: * 469: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 470: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 471: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 472: * 473: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 474: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 475: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 476: * 477: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, B( L1, L1 ), 478: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 479: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 480: IF( IERR.NE.0 ) 481: $ INFO = 1 482: * 483: IF( SCALOC.NE.ONE ) THEN 484: DO 90 J = 1, N 485: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 486: 90 CONTINUE 487: SCALE = SCALE*SCALOC 488: END IF 489: C( K1, L1 ) = X( 1, 1 ) 490: C( K1, L2 ) = X( 2, 1 ) 491: * 492: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 493: * 494: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 495: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L1 ), 1 ) 496: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 497: * 498: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 499: SUMR = DDOT( L1-1, C( K1, 1 ), LDC, B( 1, L2 ), 1 ) 500: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 501: * 502: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 503: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L1 ), 1 ) 504: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 505: * 506: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 507: SUMR = DDOT( L1-1, C( K2, 1 ), LDC, B( 1, L2 ), 1 ) 508: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 509: * 510: CALL DLASY2( .TRUE., .FALSE., ISGN, 2, 2, A( K1, K1 ), 511: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 512: $ 2, XNORM, IERR ) 513: IF( IERR.NE.0 ) 514: $ INFO = 1 515: * 516: IF( SCALOC.NE.ONE ) THEN 517: DO 100 J = 1, N 518: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 519: 100 CONTINUE 520: SCALE = SCALE*SCALOC 521: END IF 522: C( K1, L1 ) = X( 1, 1 ) 523: C( K1, L2 ) = X( 1, 2 ) 524: C( K2, L1 ) = X( 2, 1 ) 525: C( K2, L2 ) = X( 2, 2 ) 526: END IF 527: * 528: 110 CONTINUE 529: 120 CONTINUE 530: * 531: ELSE IF( .NOT.NOTRNA .AND. .NOT.NOTRNB ) THEN 532: * 533: * Solve A'*X + ISGN*X*B' = scale*C. 534: * 535: * The (K,L)th block of X is determined starting from 536: * top-right corner column by column by 537: * 538: * A(K,K)'*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L) 539: * 540: * Where 541: * K-1 N 542: * R(K,L) = SUM [A(I,K)'*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)']. 543: * I=1 J=L+1 544: * 545: * Start column loop (index = L) 546: * L1 (L2): column index of the first (last) row of X(K,L) 547: * 548: LNEXT = N 549: DO 180 L = N, 1, -1 550: IF( L.GT.LNEXT ) 551: $ GO TO 180 552: IF( L.EQ.1 ) THEN 553: L1 = L 554: L2 = L 555: ELSE 556: IF( B( L, L-1 ).NE.ZERO ) THEN 557: L1 = L - 1 558: L2 = L 559: LNEXT = L - 2 560: ELSE 561: L1 = L 562: L2 = L 563: LNEXT = L - 1 564: END IF 565: END IF 566: * 567: * Start row loop (index = K) 568: * K1 (K2): row index of the first (last) row of X(K,L) 569: * 570: KNEXT = 1 571: DO 170 K = 1, M 572: IF( K.LT.KNEXT ) 573: $ GO TO 170 574: IF( K.EQ.M ) THEN 575: K1 = K 576: K2 = K 577: ELSE 578: IF( A( K+1, K ).NE.ZERO ) THEN 579: K1 = K 580: K2 = K + 1 581: KNEXT = K + 2 582: ELSE 583: K1 = K 584: K2 = K 585: KNEXT = K + 1 586: END IF 587: END IF 588: * 589: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 590: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 591: SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 592: $ B( L1, MIN( L1+1, N ) ), LDB ) 593: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 594: SCALOC = ONE 595: * 596: A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 597: DA11 = ABS( A11 ) 598: IF( DA11.LE.SMIN ) THEN 599: A11 = SMIN 600: DA11 = SMIN 601: INFO = 1 602: END IF 603: DB = ABS( VEC( 1, 1 ) ) 604: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 605: IF( DB.GT.BIGNUM*DA11 ) 606: $ SCALOC = ONE / DB 607: END IF 608: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 609: * 610: IF( SCALOC.NE.ONE ) THEN 611: DO 130 J = 1, N 612: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 613: 130 CONTINUE 614: SCALE = SCALE*SCALOC 615: END IF 616: C( K1, L1 ) = X( 1, 1 ) 617: * 618: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 619: * 620: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 621: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 622: $ B( L1, MIN( L2+1, N ) ), LDB ) 623: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 624: * 625: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 626: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 627: $ B( L1, MIN( L2+1, N ) ), LDB ) 628: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 629: * 630: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, A( K1, K1 ), 631: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 632: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 633: IF( IERR.NE.0 ) 634: $ INFO = 1 635: * 636: IF( SCALOC.NE.ONE ) THEN 637: DO 140 J = 1, N 638: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 639: 140 CONTINUE 640: SCALE = SCALE*SCALOC 641: END IF 642: C( K1, L1 ) = X( 1, 1 ) 643: C( K2, L1 ) = X( 2, 1 ) 644: * 645: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 646: * 647: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 648: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 649: $ B( L1, MIN( L2+1, N ) ), LDB ) 650: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 651: * 652: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 653: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 654: $ B( L2, MIN( L2+1, N ) ), LDB ) 655: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 656: * 657: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 658: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 659: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 660: IF( IERR.NE.0 ) 661: $ INFO = 1 662: * 663: IF( SCALOC.NE.ONE ) THEN 664: DO 150 J = 1, N 665: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 666: 150 CONTINUE 667: SCALE = SCALE*SCALOC 668: END IF 669: C( K1, L1 ) = X( 1, 1 ) 670: C( K1, L2 ) = X( 2, 1 ) 671: * 672: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 673: * 674: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L1 ), 1 ) 675: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 676: $ B( L1, MIN( L2+1, N ) ), LDB ) 677: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 678: * 679: SUML = DDOT( K1-1, A( 1, K1 ), 1, C( 1, L2 ), 1 ) 680: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 681: $ B( L2, MIN( L2+1, N ) ), LDB ) 682: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 683: * 684: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L1 ), 1 ) 685: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 686: $ B( L1, MIN( L2+1, N ) ), LDB ) 687: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 688: * 689: SUML = DDOT( K1-1, A( 1, K2 ), 1, C( 1, L2 ), 1 ) 690: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 691: $ B( L2, MIN( L2+1, N ) ), LDB ) 692: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 693: * 694: CALL DLASY2( .TRUE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 695: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 696: $ 2, XNORM, IERR ) 697: IF( IERR.NE.0 ) 698: $ INFO = 1 699: * 700: IF( SCALOC.NE.ONE ) THEN 701: DO 160 J = 1, N 702: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 703: 160 CONTINUE 704: SCALE = SCALE*SCALOC 705: END IF 706: C( K1, L1 ) = X( 1, 1 ) 707: C( K1, L2 ) = X( 1, 2 ) 708: C( K2, L1 ) = X( 2, 1 ) 709: C( K2, L2 ) = X( 2, 2 ) 710: END IF 711: * 712: 170 CONTINUE 713: 180 CONTINUE 714: * 715: ELSE IF( NOTRNA .AND. .NOT.NOTRNB ) THEN 716: * 717: * Solve A*X + ISGN*X*B' = scale*C. 718: * 719: * The (K,L)th block of X is determined starting from 720: * bottom-right corner column by column by 721: * 722: * A(K,K)*X(K,L) + ISGN*X(K,L)*B(L,L)' = C(K,L) - R(K,L) 723: * 724: * Where 725: * M N 726: * R(K,L) = SUM [A(K,I)*X(I,L)] + ISGN*SUM [X(K,J)*B(L,J)']. 727: * I=K+1 J=L+1 728: * 729: * Start column loop (index = L) 730: * L1 (L2): column index of the first (last) row of X(K,L) 731: * 732: LNEXT = N 733: DO 240 L = N, 1, -1 734: IF( L.GT.LNEXT ) 735: $ GO TO 240 736: IF( L.EQ.1 ) THEN 737: L1 = L 738: L2 = L 739: ELSE 740: IF( B( L, L-1 ).NE.ZERO ) THEN 741: L1 = L - 1 742: L2 = L 743: LNEXT = L - 2 744: ELSE 745: L1 = L 746: L2 = L 747: LNEXT = L - 1 748: END IF 749: END IF 750: * 751: * Start row loop (index = K) 752: * K1 (K2): row index of the first (last) row of X(K,L) 753: * 754: KNEXT = M 755: DO 230 K = M, 1, -1 756: IF( K.GT.KNEXT ) 757: $ GO TO 230 758: IF( K.EQ.1 ) THEN 759: K1 = K 760: K2 = K 761: ELSE 762: IF( A( K, K-1 ).NE.ZERO ) THEN 763: K1 = K - 1 764: K2 = K 765: KNEXT = K - 2 766: ELSE 767: K1 = K 768: K2 = K 769: KNEXT = K - 1 770: END IF 771: END IF 772: * 773: IF( L1.EQ.L2 .AND. K1.EQ.K2 ) THEN 774: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 775: $ C( MIN( K1+1, M ), L1 ), 1 ) 776: SUMR = DDOT( N-L1, C( K1, MIN( L1+1, N ) ), LDC, 777: $ B( L1, MIN( L1+1, N ) ), LDB ) 778: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 779: SCALOC = ONE 780: * 781: A11 = A( K1, K1 ) + SGN*B( L1, L1 ) 782: DA11 = ABS( A11 ) 783: IF( DA11.LE.SMIN ) THEN 784: A11 = SMIN 785: DA11 = SMIN 786: INFO = 1 787: END IF 788: DB = ABS( VEC( 1, 1 ) ) 789: IF( DA11.LT.ONE .AND. DB.GT.ONE ) THEN 790: IF( DB.GT.BIGNUM*DA11 ) 791: $ SCALOC = ONE / DB 792: END IF 793: X( 1, 1 ) = ( VEC( 1, 1 )*SCALOC ) / A11 794: * 795: IF( SCALOC.NE.ONE ) THEN 796: DO 190 J = 1, N 797: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 798: 190 CONTINUE 799: SCALE = SCALE*SCALOC 800: END IF 801: C( K1, L1 ) = X( 1, 1 ) 802: * 803: ELSE IF( L1.EQ.L2 .AND. K1.NE.K2 ) THEN 804: * 805: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 806: $ C( MIN( K2+1, M ), L1 ), 1 ) 807: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 808: $ B( L1, MIN( L2+1, N ) ), LDB ) 809: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 810: * 811: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 812: $ C( MIN( K2+1, M ), L1 ), 1 ) 813: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 814: $ B( L1, MIN( L2+1, N ) ), LDB ) 815: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 816: * 817: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, A( K1, K1 ), 818: $ LDA, ONE, ONE, VEC, 2, -SGN*B( L1, L1 ), 819: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 820: IF( IERR.NE.0 ) 821: $ INFO = 1 822: * 823: IF( SCALOC.NE.ONE ) THEN 824: DO 200 J = 1, N 825: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 826: 200 CONTINUE 827: SCALE = SCALE*SCALOC 828: END IF 829: C( K1, L1 ) = X( 1, 1 ) 830: C( K2, L1 ) = X( 2, 1 ) 831: * 832: ELSE IF( L1.NE.L2 .AND. K1.EQ.K2 ) THEN 833: * 834: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 835: $ C( MIN( K1+1, M ), L1 ), 1 ) 836: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 837: $ B( L1, MIN( L2+1, N ) ), LDB ) 838: VEC( 1, 1 ) = SGN*( C( K1, L1 )-( SUML+SGN*SUMR ) ) 839: * 840: SUML = DDOT( M-K1, A( K1, MIN( K1+1, M ) ), LDA, 841: $ C( MIN( K1+1, M ), L2 ), 1 ) 842: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 843: $ B( L2, MIN( L2+1, N ) ), LDB ) 844: VEC( 2, 1 ) = SGN*( C( K1, L2 )-( SUML+SGN*SUMR ) ) 845: * 846: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, B( L1, L1 ), 847: $ LDB, ONE, ONE, VEC, 2, -SGN*A( K1, K1 ), 848: $ ZERO, X, 2, SCALOC, XNORM, IERR ) 849: IF( IERR.NE.0 ) 850: $ INFO = 1 851: * 852: IF( SCALOC.NE.ONE ) THEN 853: DO 210 J = 1, N 854: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 855: 210 CONTINUE 856: SCALE = SCALE*SCALOC 857: END IF 858: C( K1, L1 ) = X( 1, 1 ) 859: C( K1, L2 ) = X( 2, 1 ) 860: * 861: ELSE IF( L1.NE.L2 .AND. K1.NE.K2 ) THEN 862: * 863: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 864: $ C( MIN( K2+1, M ), L1 ), 1 ) 865: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 866: $ B( L1, MIN( L2+1, N ) ), LDB ) 867: VEC( 1, 1 ) = C( K1, L1 ) - ( SUML+SGN*SUMR ) 868: * 869: SUML = DDOT( M-K2, A( K1, MIN( K2+1, M ) ), LDA, 870: $ C( MIN( K2+1, M ), L2 ), 1 ) 871: SUMR = DDOT( N-L2, C( K1, MIN( L2+1, N ) ), LDC, 872: $ B( L2, MIN( L2+1, N ) ), LDB ) 873: VEC( 1, 2 ) = C( K1, L2 ) - ( SUML+SGN*SUMR ) 874: * 875: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 876: $ C( MIN( K2+1, M ), L1 ), 1 ) 877: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 878: $ B( L1, MIN( L2+1, N ) ), LDB ) 879: VEC( 2, 1 ) = C( K2, L1 ) - ( SUML+SGN*SUMR ) 880: * 881: SUML = DDOT( M-K2, A( K2, MIN( K2+1, M ) ), LDA, 882: $ C( MIN( K2+1, M ), L2 ), 1 ) 883: SUMR = DDOT( N-L2, C( K2, MIN( L2+1, N ) ), LDC, 884: $ B( L2, MIN( L2+1, N ) ), LDB ) 885: VEC( 2, 2 ) = C( K2, L2 ) - ( SUML+SGN*SUMR ) 886: * 887: CALL DLASY2( .FALSE., .TRUE., ISGN, 2, 2, A( K1, K1 ), 888: $ LDA, B( L1, L1 ), LDB, VEC, 2, SCALOC, X, 889: $ 2, XNORM, IERR ) 890: IF( IERR.NE.0 ) 891: $ INFO = 1 892: * 893: IF( SCALOC.NE.ONE ) THEN 894: DO 220 J = 1, N 895: CALL DSCAL( M, SCALOC, C( 1, J ), 1 ) 896: 220 CONTINUE 897: SCALE = SCALE*SCALOC 898: END IF 899: C( K1, L1 ) = X( 1, 1 ) 900: C( K1, L2 ) = X( 1, 2 ) 901: C( K2, L1 ) = X( 2, 1 ) 902: C( K2, L2 ) = X( 2, 2 ) 903: END IF 904: * 905: 230 CONTINUE 906: 240 CONTINUE 907: * 908: END IF 909: * 910: RETURN 911: * 912: * End of DTRSYL 913: * 914: END