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