Annotation of rpl/lapack/lapack/dtfsm.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DTFSM( TRANSR, SIDE, UPLO, TRANS, DIAG, M, N, ALPHA, A,
! 2: + B, LDB )
! 3: *
! 4: * -- LAPACK routine (version 3.2.2) --
! 5: *
! 6: * -- Contributed by Fred Gustavson of the IBM Watson Research Center --
! 7: * -- June 2010 --
! 8: *
! 9: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 10: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 11: *
! 12: * ..
! 13: * .. Scalar Arguments ..
! 14: CHARACTER TRANSR, DIAG, SIDE, TRANS, UPLO
! 15: INTEGER LDB, M, N
! 16: DOUBLE PRECISION ALPHA
! 17: * ..
! 18: * .. Array Arguments ..
! 19: DOUBLE PRECISION A( 0: * ), B( 0: LDB-1, 0: * )
! 20: * ..
! 21: *
! 22: * Purpose
! 23: * =======
! 24: *
! 25: * Level 3 BLAS like routine for A in RFP Format.
! 26: *
! 27: * DTFSM solves the matrix equation
! 28: *
! 29: * op( A )*X = alpha*B or X*op( A ) = alpha*B
! 30: *
! 31: * where alpha is a scalar, X and B are m by n matrices, A is a unit, or
! 32: * non-unit, upper or lower triangular matrix and op( A ) is one of
! 33: *
! 34: * op( A ) = A or op( A ) = A'.
! 35: *
! 36: * A is in Rectangular Full Packed (RFP) Format.
! 37: *
! 38: * The matrix X is overwritten on B.
! 39: *
! 40: * Arguments
! 41: * ==========
! 42: *
! 43: * TRANSR (input) CHARACTER
! 44: * = 'N': The Normal Form of RFP A is stored;
! 45: * = 'T': The Transpose Form of RFP A is stored.
! 46: *
! 47: * SIDE (input) CHARACTER
! 48: * On entry, SIDE specifies whether op( A ) appears on the left
! 49: * or right of X as follows:
! 50: *
! 51: * SIDE = 'L' or 'l' op( A )*X = alpha*B.
! 52: *
! 53: * SIDE = 'R' or 'r' X*op( A ) = alpha*B.
! 54: *
! 55: * Unchanged on exit.
! 56: *
! 57: * UPLO (input) CHARACTER
! 58: * On entry, UPLO specifies whether the RFP matrix A came from
! 59: * an upper or lower triangular matrix as follows:
! 60: * UPLO = 'U' or 'u' RFP A came from an upper triangular matrix
! 61: * UPLO = 'L' or 'l' RFP A came from a lower triangular matrix
! 62: *
! 63: * Unchanged on exit.
! 64: *
! 65: * TRANS (input) CHARACTER
! 66: * On entry, TRANS specifies the form of op( A ) to be used
! 67: * in the matrix multiplication as follows:
! 68: *
! 69: * TRANS = 'N' or 'n' op( A ) = A.
! 70: *
! 71: * TRANS = 'T' or 't' op( A ) = A'.
! 72: *
! 73: * Unchanged on exit.
! 74: *
! 75: * DIAG (input) CHARACTER
! 76: * On entry, DIAG specifies whether or not RFP A is unit
! 77: * triangular as follows:
! 78: *
! 79: * DIAG = 'U' or 'u' A is assumed to be unit triangular.
! 80: *
! 81: * DIAG = 'N' or 'n' A is not assumed to be unit
! 82: * triangular.
! 83: *
! 84: * Unchanged on exit.
! 85: *
! 86: * M (input) INTEGER
! 87: * On entry, M specifies the number of rows of B. M must be at
! 88: * least zero.
! 89: * Unchanged on exit.
! 90: *
! 91: * N (input) INTEGER
! 92: * On entry, N specifies the number of columns of B. N must be
! 93: * at least zero.
! 94: * Unchanged on exit.
! 95: *
! 96: * ALPHA (input) DOUBLE PRECISION
! 97: * On entry, ALPHA specifies the scalar alpha. When alpha is
! 98: * zero then A is not referenced and B need not be set before
! 99: * entry.
! 100: * Unchanged on exit.
! 101: *
! 102: * A (input) DOUBLE PRECISION array, dimension (NT)
! 103: * NT = N*(N+1)/2. On entry, the matrix A in RFP Format.
! 104: * RFP Format is described by TRANSR, UPLO and N as follows:
! 105: * If TRANSR='N' then RFP A is (0:N,0:K-1) when N is even;
! 106: * K=N/2. RFP A is (0:N-1,0:K) when N is odd; K=N/2. If
! 107: * TRANSR = 'T' then RFP is the transpose of RFP A as
! 108: * defined when TRANSR = 'N'. The contents of RFP A are defined
! 109: * by UPLO as follows: If UPLO = 'U' the RFP A contains the NT
! 110: * elements of upper packed A either in normal or
! 111: * transpose Format. If UPLO = 'L' the RFP A contains
! 112: * the NT elements of lower packed A either in normal or
! 113: * transpose Format. The LDA of RFP A is (N+1)/2 when
! 114: * TRANSR = 'T'. When TRANSR is 'N' the LDA is N+1 when N is
! 115: * even and is N when is odd.
! 116: * See the Note below for more details. Unchanged on exit.
! 117: *
! 118: * B (input/output) DOUBLE PRECISION array, dimension (LDB,N)
! 119: * Before entry, the leading m by n part of the array B must
! 120: * contain the right-hand side matrix B, and on exit is
! 121: * overwritten by the solution matrix X.
! 122: *
! 123: * LDB (input) INTEGER
! 124: * On entry, LDB specifies the first dimension of B as declared
! 125: * in the calling (sub) program. LDB must be at least
! 126: * max( 1, m ).
! 127: * Unchanged on exit.
! 128: *
! 129: * Further Details
! 130: * ===============
! 131: *
! 132: * We first consider Rectangular Full Packed (RFP) Format when N is
! 133: * even. We give an example where N = 6.
! 134: *
! 135: * AP is Upper AP is Lower
! 136: *
! 137: * 00 01 02 03 04 05 00
! 138: * 11 12 13 14 15 10 11
! 139: * 22 23 24 25 20 21 22
! 140: * 33 34 35 30 31 32 33
! 141: * 44 45 40 41 42 43 44
! 142: * 55 50 51 52 53 54 55
! 143: *
! 144: *
! 145: * Let TRANSR = 'N'. RFP holds AP as follows:
! 146: * For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
! 147: * three columns of AP upper. The lower triangle A(4:6,0:2) consists of
! 148: * the transpose of the first three columns of AP upper.
! 149: * For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
! 150: * three columns of AP lower. The upper triangle A(0:2,0:2) consists of
! 151: * the transpose of the last three columns of AP lower.
! 152: * This covers the case N even and TRANSR = 'N'.
! 153: *
! 154: * RFP A RFP A
! 155: *
! 156: * 03 04 05 33 43 53
! 157: * 13 14 15 00 44 54
! 158: * 23 24 25 10 11 55
! 159: * 33 34 35 20 21 22
! 160: * 00 44 45 30 31 32
! 161: * 01 11 55 40 41 42
! 162: * 02 12 22 50 51 52
! 163: *
! 164: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
! 165: * transpose of RFP A above. One therefore gets:
! 166: *
! 167: *
! 168: * RFP A RFP A
! 169: *
! 170: * 03 13 23 33 00 01 02 33 00 10 20 30 40 50
! 171: * 04 14 24 34 44 11 12 43 44 11 21 31 41 51
! 172: * 05 15 25 35 45 55 22 53 54 55 22 32 42 52
! 173: *
! 174: *
! 175: * We then consider Rectangular Full Packed (RFP) Format when N is
! 176: * odd. We give an example where N = 5.
! 177: *
! 178: * AP is Upper AP is Lower
! 179: *
! 180: * 00 01 02 03 04 00
! 181: * 11 12 13 14 10 11
! 182: * 22 23 24 20 21 22
! 183: * 33 34 30 31 32 33
! 184: * 44 40 41 42 43 44
! 185: *
! 186: *
! 187: * Let TRANSR = 'N'. RFP holds AP as follows:
! 188: * For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
! 189: * three columns of AP upper. The lower triangle A(3:4,0:1) consists of
! 190: * the transpose of the first two columns of AP upper.
! 191: * For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
! 192: * three columns of AP lower. The upper triangle A(0:1,1:2) consists of
! 193: * the transpose of the last two columns of AP lower.
! 194: * This covers the case N odd and TRANSR = 'N'.
! 195: *
! 196: * RFP A RFP A
! 197: *
! 198: * 02 03 04 00 33 43
! 199: * 12 13 14 10 11 44
! 200: * 22 23 24 20 21 22
! 201: * 00 33 34 30 31 32
! 202: * 01 11 44 40 41 42
! 203: *
! 204: * Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
! 205: * transpose of RFP A above. One therefore gets:
! 206: *
! 207: * RFP A RFP A
! 208: *
! 209: * 02 12 22 00 01 00 10 20 30 40 50
! 210: * 03 13 23 33 11 33 11 21 31 41 51
! 211: * 04 14 24 34 44 43 44 22 32 42 52
! 212: *
! 213: * Reference
! 214: * =========
! 215: *
! 216: * =====================================================================
! 217: *
! 218: * ..
! 219: * .. Parameters ..
! 220: DOUBLE PRECISION ONE, ZERO
! 221: PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
! 222: * ..
! 223: * .. Local Scalars ..
! 224: LOGICAL LOWER, LSIDE, MISODD, NISODD, NORMALTRANSR,
! 225: + NOTRANS
! 226: INTEGER M1, M2, N1, N2, K, INFO, I, J
! 227: * ..
! 228: * .. External Functions ..
! 229: LOGICAL LSAME
! 230: EXTERNAL LSAME
! 231: * ..
! 232: * .. External Subroutines ..
! 233: EXTERNAL XERBLA, DGEMM, DTRSM
! 234: * ..
! 235: * .. Intrinsic Functions ..
! 236: INTRINSIC MAX, MOD
! 237: * ..
! 238: * .. Executable Statements ..
! 239: *
! 240: * Test the input parameters.
! 241: *
! 242: INFO = 0
! 243: NORMALTRANSR = LSAME( TRANSR, 'N' )
! 244: LSIDE = LSAME( SIDE, 'L' )
! 245: LOWER = LSAME( UPLO, 'L' )
! 246: NOTRANS = LSAME( TRANS, 'N' )
! 247: IF( .NOT.NORMALTRANSR .AND. .NOT.LSAME( TRANSR, 'T' ) ) THEN
! 248: INFO = -1
! 249: ELSE IF( .NOT.LSIDE .AND. .NOT.LSAME( SIDE, 'R' ) ) THEN
! 250: INFO = -2
! 251: ELSE IF( .NOT.LOWER .AND. .NOT.LSAME( UPLO, 'U' ) ) THEN
! 252: INFO = -3
! 253: ELSE IF( .NOT.NOTRANS .AND. .NOT.LSAME( TRANS, 'T' ) ) THEN
! 254: INFO = -4
! 255: ELSE IF( .NOT.LSAME( DIAG, 'N' ) .AND. .NOT.LSAME( DIAG, 'U' ) )
! 256: + THEN
! 257: INFO = -5
! 258: ELSE IF( M.LT.0 ) THEN
! 259: INFO = -6
! 260: ELSE IF( N.LT.0 ) THEN
! 261: INFO = -7
! 262: ELSE IF( LDB.LT.MAX( 1, M ) ) THEN
! 263: INFO = -11
! 264: END IF
! 265: IF( INFO.NE.0 ) THEN
! 266: CALL XERBLA( 'DTFSM ', -INFO )
! 267: RETURN
! 268: END IF
! 269: *
! 270: * Quick return when ( (N.EQ.0).OR.(M.EQ.0) )
! 271: *
! 272: IF( ( M.EQ.0 ) .OR. ( N.EQ.0 ) )
! 273: + RETURN
! 274: *
! 275: * Quick return when ALPHA.EQ.(0D+0)
! 276: *
! 277: IF( ALPHA.EQ.ZERO ) THEN
! 278: DO 20 J = 0, N - 1
! 279: DO 10 I = 0, M - 1
! 280: B( I, J ) = ZERO
! 281: 10 CONTINUE
! 282: 20 CONTINUE
! 283: RETURN
! 284: END IF
! 285: *
! 286: IF( LSIDE ) THEN
! 287: *
! 288: * SIDE = 'L'
! 289: *
! 290: * A is M-by-M.
! 291: * If M is odd, set NISODD = .TRUE., and M1 and M2.
! 292: * If M is even, NISODD = .FALSE., and M.
! 293: *
! 294: IF( MOD( M, 2 ).EQ.0 ) THEN
! 295: MISODD = .FALSE.
! 296: K = M / 2
! 297: ELSE
! 298: MISODD = .TRUE.
! 299: IF( LOWER ) THEN
! 300: M2 = M / 2
! 301: M1 = M - M2
! 302: ELSE
! 303: M1 = M / 2
! 304: M2 = M - M1
! 305: END IF
! 306: END IF
! 307: *
! 308: *
! 309: IF( MISODD ) THEN
! 310: *
! 311: * SIDE = 'L' and N is odd
! 312: *
! 313: IF( NORMALTRANSR ) THEN
! 314: *
! 315: * SIDE = 'L', N is odd, and TRANSR = 'N'
! 316: *
! 317: IF( LOWER ) THEN
! 318: *
! 319: * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'L'
! 320: *
! 321: IF( NOTRANS ) THEN
! 322: *
! 323: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
! 324: * TRANS = 'N'
! 325: *
! 326: IF( M.EQ.1 ) THEN
! 327: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
! 328: + A, M, B, LDB )
! 329: ELSE
! 330: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
! 331: + A( 0 ), M, B, LDB )
! 332: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( M1 ),
! 333: + M, B, LDB, ALPHA, B( M1, 0 ), LDB )
! 334: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
! 335: + A( M ), M, B( M1, 0 ), LDB )
! 336: END IF
! 337: *
! 338: ELSE
! 339: *
! 340: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'L', and
! 341: * TRANS = 'T'
! 342: *
! 343: IF( M.EQ.1 ) THEN
! 344: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ALPHA,
! 345: + A( 0 ), M, B, LDB )
! 346: ELSE
! 347: CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
! 348: + A( M ), M, B( M1, 0 ), LDB )
! 349: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( M1 ),
! 350: + M, B( M1, 0 ), LDB, ALPHA, B, LDB )
! 351: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
! 352: + A( 0 ), M, B, LDB )
! 353: END IF
! 354: *
! 355: END IF
! 356: *
! 357: ELSE
! 358: *
! 359: * SIDE ='L', N is odd, TRANSR = 'N', and UPLO = 'U'
! 360: *
! 361: IF( .NOT.NOTRANS ) THEN
! 362: *
! 363: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
! 364: * TRANS = 'N'
! 365: *
! 366: CALL DTRSM( 'L', 'L', 'N', DIAG, M1, N, ALPHA,
! 367: + A( M2 ), M, B, LDB )
! 368: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE, A( 0 ), M,
! 369: + B, LDB, ALPHA, B( M1, 0 ), LDB )
! 370: CALL DTRSM( 'L', 'U', 'T', DIAG, M2, N, ONE,
! 371: + A( M1 ), M, B( M1, 0 ), LDB )
! 372: *
! 373: ELSE
! 374: *
! 375: * SIDE ='L', N is odd, TRANSR = 'N', UPLO = 'U', and
! 376: * TRANS = 'T'
! 377: *
! 378: CALL DTRSM( 'L', 'U', 'N', DIAG, M2, N, ALPHA,
! 379: + A( M1 ), M, B( M1, 0 ), LDB )
! 380: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE, A( 0 ), M,
! 381: + B( M1, 0 ), LDB, ALPHA, B, LDB )
! 382: CALL DTRSM( 'L', 'L', 'T', DIAG, M1, N, ONE,
! 383: + A( M2 ), M, B, LDB )
! 384: *
! 385: END IF
! 386: *
! 387: END IF
! 388: *
! 389: ELSE
! 390: *
! 391: * SIDE = 'L', N is odd, and TRANSR = 'T'
! 392: *
! 393: IF( LOWER ) THEN
! 394: *
! 395: * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'L'
! 396: *
! 397: IF( NOTRANS ) THEN
! 398: *
! 399: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
! 400: * TRANS = 'N'
! 401: *
! 402: IF( M.EQ.1 ) THEN
! 403: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
! 404: + A( 0 ), M1, B, LDB )
! 405: ELSE
! 406: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
! 407: + A( 0 ), M1, B, LDB )
! 408: CALL DGEMM( 'T', 'N', M2, N, M1, -ONE,
! 409: + A( M1*M1 ), M1, B, LDB, ALPHA,
! 410: + B( M1, 0 ), LDB )
! 411: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
! 412: + A( 1 ), M1, B( M1, 0 ), LDB )
! 413: END IF
! 414: *
! 415: ELSE
! 416: *
! 417: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'L', and
! 418: * TRANS = 'T'
! 419: *
! 420: IF( M.EQ.1 ) THEN
! 421: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ALPHA,
! 422: + A( 0 ), M1, B, LDB )
! 423: ELSE
! 424: CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
! 425: + A( 1 ), M1, B( M1, 0 ), LDB )
! 426: CALL DGEMM( 'N', 'N', M1, N, M2, -ONE,
! 427: + A( M1*M1 ), M1, B( M1, 0 ), LDB,
! 428: + ALPHA, B, LDB )
! 429: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
! 430: + A( 0 ), M1, B, LDB )
! 431: END IF
! 432: *
! 433: END IF
! 434: *
! 435: ELSE
! 436: *
! 437: * SIDE ='L', N is odd, TRANSR = 'T', and UPLO = 'U'
! 438: *
! 439: IF( .NOT.NOTRANS ) THEN
! 440: *
! 441: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
! 442: * TRANS = 'N'
! 443: *
! 444: CALL DTRSM( 'L', 'U', 'T', DIAG, M1, N, ALPHA,
! 445: + A( M2*M2 ), M2, B, LDB )
! 446: CALL DGEMM( 'N', 'N', M2, N, M1, -ONE, A( 0 ), M2,
! 447: + B, LDB, ALPHA, B( M1, 0 ), LDB )
! 448: CALL DTRSM( 'L', 'L', 'N', DIAG, M2, N, ONE,
! 449: + A( M1*M2 ), M2, B( M1, 0 ), LDB )
! 450: *
! 451: ELSE
! 452: *
! 453: * SIDE ='L', N is odd, TRANSR = 'T', UPLO = 'U', and
! 454: * TRANS = 'T'
! 455: *
! 456: CALL DTRSM( 'L', 'L', 'T', DIAG, M2, N, ALPHA,
! 457: + A( M1*M2 ), M2, B( M1, 0 ), LDB )
! 458: CALL DGEMM( 'T', 'N', M1, N, M2, -ONE, A( 0 ), M2,
! 459: + B( M1, 0 ), LDB, ALPHA, B, LDB )
! 460: CALL DTRSM( 'L', 'U', 'N', DIAG, M1, N, ONE,
! 461: + A( M2*M2 ), M2, B, LDB )
! 462: *
! 463: END IF
! 464: *
! 465: END IF
! 466: *
! 467: END IF
! 468: *
! 469: ELSE
! 470: *
! 471: * SIDE = 'L' and N is even
! 472: *
! 473: IF( NORMALTRANSR ) THEN
! 474: *
! 475: * SIDE = 'L', N is even, and TRANSR = 'N'
! 476: *
! 477: IF( LOWER ) THEN
! 478: *
! 479: * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'L'
! 480: *
! 481: IF( NOTRANS ) THEN
! 482: *
! 483: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
! 484: * and TRANS = 'N'
! 485: *
! 486: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
! 487: + A( 1 ), M+1, B, LDB )
! 488: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( K+1 ),
! 489: + M+1, B, LDB, ALPHA, B( K, 0 ), LDB )
! 490: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
! 491: + A( 0 ), M+1, B( K, 0 ), LDB )
! 492: *
! 493: ELSE
! 494: *
! 495: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'L',
! 496: * and TRANS = 'T'
! 497: *
! 498: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
! 499: + A( 0 ), M+1, B( K, 0 ), LDB )
! 500: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( K+1 ),
! 501: + M+1, B( K, 0 ), LDB, ALPHA, B, LDB )
! 502: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
! 503: + A( 1 ), M+1, B, LDB )
! 504: *
! 505: END IF
! 506: *
! 507: ELSE
! 508: *
! 509: * SIDE ='L', N is even, TRANSR = 'N', and UPLO = 'U'
! 510: *
! 511: IF( .NOT.NOTRANS ) THEN
! 512: *
! 513: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
! 514: * and TRANS = 'N'
! 515: *
! 516: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ALPHA,
! 517: + A( K+1 ), M+1, B, LDB )
! 518: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), M+1,
! 519: + B, LDB, ALPHA, B( K, 0 ), LDB )
! 520: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ONE,
! 521: + A( K ), M+1, B( K, 0 ), LDB )
! 522: *
! 523: ELSE
! 524: *
! 525: * SIDE ='L', N is even, TRANSR = 'N', UPLO = 'U',
! 526: * and TRANS = 'T'
! 527: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ALPHA,
! 528: + A( K ), M+1, B( K, 0 ), LDB )
! 529: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), M+1,
! 530: + B( K, 0 ), LDB, ALPHA, B, LDB )
! 531: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ONE,
! 532: + A( K+1 ), M+1, B, LDB )
! 533: *
! 534: END IF
! 535: *
! 536: END IF
! 537: *
! 538: ELSE
! 539: *
! 540: * SIDE = 'L', N is even, and TRANSR = 'T'
! 541: *
! 542: IF( LOWER ) THEN
! 543: *
! 544: * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'L'
! 545: *
! 546: IF( NOTRANS ) THEN
! 547: *
! 548: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
! 549: * and TRANS = 'N'
! 550: *
! 551: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
! 552: + A( K ), K, B, LDB )
! 553: CALL DGEMM( 'T', 'N', K, N, K, -ONE,
! 554: + A( K*( K+1 ) ), K, B, LDB, ALPHA,
! 555: + B( K, 0 ), LDB )
! 556: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
! 557: + A( 0 ), K, B( K, 0 ), LDB )
! 558: *
! 559: ELSE
! 560: *
! 561: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'L',
! 562: * and TRANS = 'T'
! 563: *
! 564: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
! 565: + A( 0 ), K, B( K, 0 ), LDB )
! 566: CALL DGEMM( 'N', 'N', K, N, K, -ONE,
! 567: + A( K*( K+1 ) ), K, B( K, 0 ), LDB,
! 568: + ALPHA, B, LDB )
! 569: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
! 570: + A( K ), K, B, LDB )
! 571: *
! 572: END IF
! 573: *
! 574: ELSE
! 575: *
! 576: * SIDE ='L', N is even, TRANSR = 'T', and UPLO = 'U'
! 577: *
! 578: IF( .NOT.NOTRANS ) THEN
! 579: *
! 580: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
! 581: * and TRANS = 'N'
! 582: *
! 583: CALL DTRSM( 'L', 'U', 'T', DIAG, K, N, ALPHA,
! 584: + A( K*( K+1 ) ), K, B, LDB )
! 585: CALL DGEMM( 'N', 'N', K, N, K, -ONE, A( 0 ), K, B,
! 586: + LDB, ALPHA, B( K, 0 ), LDB )
! 587: CALL DTRSM( 'L', 'L', 'N', DIAG, K, N, ONE,
! 588: + A( K*K ), K, B( K, 0 ), LDB )
! 589: *
! 590: ELSE
! 591: *
! 592: * SIDE ='L', N is even, TRANSR = 'T', UPLO = 'U',
! 593: * and TRANS = 'T'
! 594: *
! 595: CALL DTRSM( 'L', 'L', 'T', DIAG, K, N, ALPHA,
! 596: + A( K*K ), K, B( K, 0 ), LDB )
! 597: CALL DGEMM( 'T', 'N', K, N, K, -ONE, A( 0 ), K,
! 598: + B( K, 0 ), LDB, ALPHA, B, LDB )
! 599: CALL DTRSM( 'L', 'U', 'N', DIAG, K, N, ONE,
! 600: + A( K*( K+1 ) ), K, B, LDB )
! 601: *
! 602: END IF
! 603: *
! 604: END IF
! 605: *
! 606: END IF
! 607: *
! 608: END IF
! 609: *
! 610: ELSE
! 611: *
! 612: * SIDE = 'R'
! 613: *
! 614: * A is N-by-N.
! 615: * If N is odd, set NISODD = .TRUE., and N1 and N2.
! 616: * If N is even, NISODD = .FALSE., and K.
! 617: *
! 618: IF( MOD( N, 2 ).EQ.0 ) THEN
! 619: NISODD = .FALSE.
! 620: K = N / 2
! 621: ELSE
! 622: NISODD = .TRUE.
! 623: IF( LOWER ) THEN
! 624: N2 = N / 2
! 625: N1 = N - N2
! 626: ELSE
! 627: N1 = N / 2
! 628: N2 = N - N1
! 629: END IF
! 630: END IF
! 631: *
! 632: IF( NISODD ) THEN
! 633: *
! 634: * SIDE = 'R' and N is odd
! 635: *
! 636: IF( NORMALTRANSR ) THEN
! 637: *
! 638: * SIDE = 'R', N is odd, and TRANSR = 'N'
! 639: *
! 640: IF( LOWER ) THEN
! 641: *
! 642: * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'L'
! 643: *
! 644: IF( NOTRANS ) THEN
! 645: *
! 646: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
! 647: * TRANS = 'N'
! 648: *
! 649: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
! 650: + A( N ), N, B( 0, N1 ), LDB )
! 651: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
! 652: + LDB, A( N1 ), N, ALPHA, B( 0, 0 ),
! 653: + LDB )
! 654: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
! 655: + A( 0 ), N, B( 0, 0 ), LDB )
! 656: *
! 657: ELSE
! 658: *
! 659: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'L', and
! 660: * TRANS = 'T'
! 661: *
! 662: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
! 663: + A( 0 ), N, B( 0, 0 ), LDB )
! 664: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
! 665: + LDB, A( N1 ), N, ALPHA, B( 0, N1 ),
! 666: + LDB )
! 667: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
! 668: + A( N ), N, B( 0, N1 ), LDB )
! 669: *
! 670: END IF
! 671: *
! 672: ELSE
! 673: *
! 674: * SIDE ='R', N is odd, TRANSR = 'N', and UPLO = 'U'
! 675: *
! 676: IF( NOTRANS ) THEN
! 677: *
! 678: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
! 679: * TRANS = 'N'
! 680: *
! 681: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N1, ALPHA,
! 682: + A( N2 ), N, B( 0, 0 ), LDB )
! 683: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
! 684: + LDB, A( 0 ), N, ALPHA, B( 0, N1 ),
! 685: + LDB )
! 686: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N2, ONE,
! 687: + A( N1 ), N, B( 0, N1 ), LDB )
! 688: *
! 689: ELSE
! 690: *
! 691: * SIDE ='R', N is odd, TRANSR = 'N', UPLO = 'U', and
! 692: * TRANS = 'T'
! 693: *
! 694: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N2, ALPHA,
! 695: + A( N1 ), N, B( 0, N1 ), LDB )
! 696: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
! 697: + LDB, A( 0 ), N, ALPHA, B( 0, 0 ), LDB )
! 698: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N1, ONE,
! 699: + A( N2 ), N, B( 0, 0 ), LDB )
! 700: *
! 701: END IF
! 702: *
! 703: END IF
! 704: *
! 705: ELSE
! 706: *
! 707: * SIDE = 'R', N is odd, and TRANSR = 'T'
! 708: *
! 709: IF( LOWER ) THEN
! 710: *
! 711: * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'L'
! 712: *
! 713: IF( NOTRANS ) THEN
! 714: *
! 715: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
! 716: * TRANS = 'N'
! 717: *
! 718: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
! 719: + A( 1 ), N1, B( 0, N1 ), LDB )
! 720: CALL DGEMM( 'N', 'T', M, N1, N2, -ONE, B( 0, N1 ),
! 721: + LDB, A( N1*N1 ), N1, ALPHA, B( 0, 0 ),
! 722: + LDB )
! 723: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
! 724: + A( 0 ), N1, B( 0, 0 ), LDB )
! 725: *
! 726: ELSE
! 727: *
! 728: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'L', and
! 729: * TRANS = 'T'
! 730: *
! 731: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
! 732: + A( 0 ), N1, B( 0, 0 ), LDB )
! 733: CALL DGEMM( 'N', 'N', M, N2, N1, -ONE, B( 0, 0 ),
! 734: + LDB, A( N1*N1 ), N1, ALPHA, B( 0, N1 ),
! 735: + LDB )
! 736: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
! 737: + A( 1 ), N1, B( 0, N1 ), LDB )
! 738: *
! 739: END IF
! 740: *
! 741: ELSE
! 742: *
! 743: * SIDE ='R', N is odd, TRANSR = 'T', and UPLO = 'U'
! 744: *
! 745: IF( NOTRANS ) THEN
! 746: *
! 747: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
! 748: * TRANS = 'N'
! 749: *
! 750: CALL DTRSM( 'R', 'U', 'N', DIAG, M, N1, ALPHA,
! 751: + A( N2*N2 ), N2, B( 0, 0 ), LDB )
! 752: CALL DGEMM( 'N', 'T', M, N2, N1, -ONE, B( 0, 0 ),
! 753: + LDB, A( 0 ), N2, ALPHA, B( 0, N1 ),
! 754: + LDB )
! 755: CALL DTRSM( 'R', 'L', 'T', DIAG, M, N2, ONE,
! 756: + A( N1*N2 ), N2, B( 0, N1 ), LDB )
! 757: *
! 758: ELSE
! 759: *
! 760: * SIDE ='R', N is odd, TRANSR = 'T', UPLO = 'U', and
! 761: * TRANS = 'T'
! 762: *
! 763: CALL DTRSM( 'R', 'L', 'N', DIAG, M, N2, ALPHA,
! 764: + A( N1*N2 ), N2, B( 0, N1 ), LDB )
! 765: CALL DGEMM( 'N', 'N', M, N1, N2, -ONE, B( 0, N1 ),
! 766: + LDB, A( 0 ), N2, ALPHA, B( 0, 0 ),
! 767: + LDB )
! 768: CALL DTRSM( 'R', 'U', 'T', DIAG, M, N1, ONE,
! 769: + A( N2*N2 ), N2, B( 0, 0 ), LDB )
! 770: *
! 771: END IF
! 772: *
! 773: END IF
! 774: *
! 775: END IF
! 776: *
! 777: ELSE
! 778: *
! 779: * SIDE = 'R' and N is even
! 780: *
! 781: IF( NORMALTRANSR ) THEN
! 782: *
! 783: * SIDE = 'R', N is even, and TRANSR = 'N'
! 784: *
! 785: IF( LOWER ) THEN
! 786: *
! 787: * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'L'
! 788: *
! 789: IF( NOTRANS ) THEN
! 790: *
! 791: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
! 792: * and TRANS = 'N'
! 793: *
! 794: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
! 795: + A( 0 ), N+1, B( 0, K ), LDB )
! 796: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
! 797: + LDB, A( K+1 ), N+1, ALPHA, B( 0, 0 ),
! 798: + LDB )
! 799: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
! 800: + A( 1 ), N+1, B( 0, 0 ), LDB )
! 801: *
! 802: ELSE
! 803: *
! 804: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'L',
! 805: * and TRANS = 'T'
! 806: *
! 807: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
! 808: + A( 1 ), N+1, B( 0, 0 ), LDB )
! 809: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
! 810: + LDB, A( K+1 ), N+1, ALPHA, B( 0, K ),
! 811: + LDB )
! 812: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
! 813: + A( 0 ), N+1, B( 0, K ), LDB )
! 814: *
! 815: END IF
! 816: *
! 817: ELSE
! 818: *
! 819: * SIDE ='R', N is even, TRANSR = 'N', and UPLO = 'U'
! 820: *
! 821: IF( NOTRANS ) THEN
! 822: *
! 823: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
! 824: * and TRANS = 'N'
! 825: *
! 826: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ALPHA,
! 827: + A( K+1 ), N+1, B( 0, 0 ), LDB )
! 828: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
! 829: + LDB, A( 0 ), N+1, ALPHA, B( 0, K ),
! 830: + LDB )
! 831: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ONE,
! 832: + A( K ), N+1, B( 0, K ), LDB )
! 833: *
! 834: ELSE
! 835: *
! 836: * SIDE ='R', N is even, TRANSR = 'N', UPLO = 'U',
! 837: * and TRANS = 'T'
! 838: *
! 839: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ALPHA,
! 840: + A( K ), N+1, B( 0, K ), LDB )
! 841: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
! 842: + LDB, A( 0 ), N+1, ALPHA, B( 0, 0 ),
! 843: + LDB )
! 844: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ONE,
! 845: + A( K+1 ), N+1, B( 0, 0 ), LDB )
! 846: *
! 847: END IF
! 848: *
! 849: END IF
! 850: *
! 851: ELSE
! 852: *
! 853: * SIDE = 'R', N is even, and TRANSR = 'T'
! 854: *
! 855: IF( LOWER ) THEN
! 856: *
! 857: * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'L'
! 858: *
! 859: IF( NOTRANS ) THEN
! 860: *
! 861: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
! 862: * and TRANS = 'N'
! 863: *
! 864: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
! 865: + A( 0 ), K, B( 0, K ), LDB )
! 866: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, K ),
! 867: + LDB, A( ( K+1 )*K ), K, ALPHA,
! 868: + B( 0, 0 ), LDB )
! 869: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
! 870: + A( K ), K, B( 0, 0 ), LDB )
! 871: *
! 872: ELSE
! 873: *
! 874: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'L',
! 875: * and TRANS = 'T'
! 876: *
! 877: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
! 878: + A( K ), K, B( 0, 0 ), LDB )
! 879: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, 0 ),
! 880: + LDB, A( ( K+1 )*K ), K, ALPHA,
! 881: + B( 0, K ), LDB )
! 882: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
! 883: + A( 0 ), K, B( 0, K ), LDB )
! 884: *
! 885: END IF
! 886: *
! 887: ELSE
! 888: *
! 889: * SIDE ='R', N is even, TRANSR = 'T', and UPLO = 'U'
! 890: *
! 891: IF( NOTRANS ) THEN
! 892: *
! 893: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
! 894: * and TRANS = 'N'
! 895: *
! 896: CALL DTRSM( 'R', 'U', 'N', DIAG, M, K, ALPHA,
! 897: + A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
! 898: CALL DGEMM( 'N', 'T', M, K, K, -ONE, B( 0, 0 ),
! 899: + LDB, A( 0 ), K, ALPHA, B( 0, K ), LDB )
! 900: CALL DTRSM( 'R', 'L', 'T', DIAG, M, K, ONE,
! 901: + A( K*K ), K, B( 0, K ), LDB )
! 902: *
! 903: ELSE
! 904: *
! 905: * SIDE ='R', N is even, TRANSR = 'T', UPLO = 'U',
! 906: * and TRANS = 'T'
! 907: *
! 908: CALL DTRSM( 'R', 'L', 'N', DIAG, M, K, ALPHA,
! 909: + A( K*K ), K, B( 0, K ), LDB )
! 910: CALL DGEMM( 'N', 'N', M, K, K, -ONE, B( 0, K ),
! 911: + LDB, A( 0 ), K, ALPHA, B( 0, 0 ), LDB )
! 912: CALL DTRSM( 'R', 'U', 'T', DIAG, M, K, ONE,
! 913: + A( ( K+1 )*K ), K, B( 0, 0 ), LDB )
! 914: *
! 915: END IF
! 916: *
! 917: END IF
! 918: *
! 919: END IF
! 920: *
! 921: END IF
! 922: END IF
! 923: *
! 924: RETURN
! 925: *
! 926: * End of DTFSM
! 927: *
! 928: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtfsm.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack