Annotation of rpl/lapack/lapack/dtrsyl.f, revision 1.1
1.1 ! bertrand 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
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtrsyl.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack