File:  [local] / rpl / lapack / blas / dtrmm.f
Revision 1.7: download - view: text, annotated - select for diffs - revision graph
Fri Jul 22 07:38:02 2011 UTC (12 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_3, rpl-4_1_2, rpl-4_1_1, HEAD
En route vers la 4.4.1.

    1:       SUBROUTINE DTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
    2: *     .. Scalar Arguments ..
    3:       DOUBLE PRECISION ALPHA
    4:       INTEGER LDA,LDB,M,N
    5:       CHARACTER DIAG,SIDE,TRANSA,UPLO
    6: *     ..
    7: *     .. Array Arguments ..
    8:       DOUBLE PRECISION A(LDA,*),B(LDB,*)
    9: *     ..
   10: *
   11: *  Purpose
   12: *  =======
   13: *
   14: *  DTRMM  performs one of the matrix-matrix operations
   15: *
   16: *     B := alpha*op( A )*B,   or   B := alpha*B*op( A ),
   17: *
   18: *  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
   19: *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
   20: *
   21: *     op( A ) = A   or   op( A ) = A**T.
   22: *
   23: *  Arguments
   24: *  ==========
   25: *
   26: *  SIDE   - CHARACTER*1.
   27: *           On entry,  SIDE specifies whether  op( A ) multiplies B from
   28: *           the left or right as follows:
   29: *
   30: *              SIDE = 'L' or 'l'   B := alpha*op( A )*B.
   31: *
   32: *              SIDE = 'R' or 'r'   B := alpha*B*op( A ).
   33: *
   34: *           Unchanged on exit.
   35: *
   36: *  UPLO   - CHARACTER*1.
   37: *           On entry, UPLO specifies whether the matrix A is an upper or
   38: *           lower triangular matrix as follows:
   39: *
   40: *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
   41: *
   42: *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
   43: *
   44: *           Unchanged on exit.
   45: *
   46: *  TRANSA - CHARACTER*1.
   47: *           On entry, TRANSA specifies the form of op( A ) to be used in
   48: *           the matrix multiplication as follows:
   49: *
   50: *              TRANSA = 'N' or 'n'   op( A ) = A.
   51: *
   52: *              TRANSA = 'T' or 't'   op( A ) = A**T.
   53: *
   54: *              TRANSA = 'C' or 'c'   op( A ) = A**T.
   55: *
   56: *           Unchanged on exit.
   57: *
   58: *  DIAG   - CHARACTER*1.
   59: *           On entry, DIAG specifies whether or not A is unit triangular
   60: *           as follows:
   61: *
   62: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
   63: *
   64: *              DIAG = 'N' or 'n'   A is not assumed to be unit
   65: *                                  triangular.
   66: *
   67: *           Unchanged on exit.
   68: *
   69: *  M      - INTEGER.
   70: *           On entry, M specifies the number of rows of B. M must be at
   71: *           least zero.
   72: *           Unchanged on exit.
   73: *
   74: *  N      - INTEGER.
   75: *           On entry, N specifies the number of columns of B.  N must be
   76: *           at least zero.
   77: *           Unchanged on exit.
   78: *
   79: *  ALPHA  - DOUBLE PRECISION.
   80: *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
   81: *           zero then  A is not referenced and  B need not be set before
   82: *           entry.
   83: *           Unchanged on exit.
   84: *
   85: *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
   86: *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
   87: *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
   88: *           upper triangular part of the array  A must contain the upper
   89: *           triangular matrix  and the strictly lower triangular part of
   90: *           A is not referenced.
   91: *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
   92: *           lower triangular part of the array  A must contain the lower
   93: *           triangular matrix  and the strictly upper triangular part of
   94: *           A is not referenced.
   95: *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
   96: *           A  are not referenced either,  but are assumed to be  unity.
   97: *           Unchanged on exit.
   98: *
   99: *  LDA    - INTEGER.
  100: *           On entry, LDA specifies the first dimension of A as declared
  101: *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
  102: *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
  103: *           then LDA must be at least max( 1, n ).
  104: *           Unchanged on exit.
  105: *
  106: *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
  107: *           Before entry,  the leading  m by n part of the array  B must
  108: *           contain the matrix  B,  and  on exit  is overwritten  by the
  109: *           transformed matrix.
  110: *
  111: *  LDB    - INTEGER.
  112: *           On entry, LDB specifies the first dimension of B as declared
  113: *           in  the  calling  (sub)  program.   LDB  must  be  at  least
  114: *           max( 1, m ).
  115: *           Unchanged on exit.
  116: *
  117: *  Further Details
  118: *  ===============
  119: *
  120: *  Level 3 Blas routine.
  121: *
  122: *  -- Written on 8-February-1989.
  123: *     Jack Dongarra, Argonne National Laboratory.
  124: *     Iain Duff, AERE Harwell.
  125: *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  126: *     Sven Hammarling, Numerical Algorithms Group Ltd.
  127: *
  128: *  =====================================================================
  129: *
  130: *     .. External Functions ..
  131:       LOGICAL LSAME
  132:       EXTERNAL LSAME
  133: *     ..
  134: *     .. External Subroutines ..
  135:       EXTERNAL XERBLA
  136: *     ..
  137: *     .. Intrinsic Functions ..
  138:       INTRINSIC MAX
  139: *     ..
  140: *     .. Local Scalars ..
  141:       DOUBLE PRECISION TEMP
  142:       INTEGER I,INFO,J,K,NROWA
  143:       LOGICAL LSIDE,NOUNIT,UPPER
  144: *     ..
  145: *     .. Parameters ..
  146:       DOUBLE PRECISION ONE,ZERO
  147:       PARAMETER (ONE=1.0D+0,ZERO=0.0D+0)
  148: *     ..
  149: *
  150: *     Test the input parameters.
  151: *
  152:       LSIDE = LSAME(SIDE,'L')
  153:       IF (LSIDE) THEN
  154:           NROWA = M
  155:       ELSE
  156:           NROWA = N
  157:       END IF
  158:       NOUNIT = LSAME(DIAG,'N')
  159:       UPPER = LSAME(UPLO,'U')
  160: *
  161:       INFO = 0
  162:       IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
  163:           INFO = 1
  164:       ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
  165:           INFO = 2
  166:       ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
  167:      +         (.NOT.LSAME(TRANSA,'T')) .AND.
  168:      +         (.NOT.LSAME(TRANSA,'C'))) THEN
  169:           INFO = 3
  170:       ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
  171:           INFO = 4
  172:       ELSE IF (M.LT.0) THEN
  173:           INFO = 5
  174:       ELSE IF (N.LT.0) THEN
  175:           INFO = 6
  176:       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
  177:           INFO = 9
  178:       ELSE IF (LDB.LT.MAX(1,M)) THEN
  179:           INFO = 11
  180:       END IF
  181:       IF (INFO.NE.0) THEN
  182:           CALL XERBLA('DTRMM ',INFO)
  183:           RETURN
  184:       END IF
  185: *
  186: *     Quick return if possible.
  187: *
  188:       IF (M.EQ.0 .OR. N.EQ.0) RETURN
  189: *
  190: *     And when  alpha.eq.zero.
  191: *
  192:       IF (ALPHA.EQ.ZERO) THEN
  193:           DO 20 J = 1,N
  194:               DO 10 I = 1,M
  195:                   B(I,J) = ZERO
  196:    10         CONTINUE
  197:    20     CONTINUE
  198:           RETURN
  199:       END IF
  200: *
  201: *     Start the operations.
  202: *
  203:       IF (LSIDE) THEN
  204:           IF (LSAME(TRANSA,'N')) THEN
  205: *
  206: *           Form  B := alpha*A*B.
  207: *
  208:               IF (UPPER) THEN
  209:                   DO 50 J = 1,N
  210:                       DO 40 K = 1,M
  211:                           IF (B(K,J).NE.ZERO) THEN
  212:                               TEMP = ALPHA*B(K,J)
  213:                               DO 30 I = 1,K - 1
  214:                                   B(I,J) = B(I,J) + TEMP*A(I,K)
  215:    30                         CONTINUE
  216:                               IF (NOUNIT) TEMP = TEMP*A(K,K)
  217:                               B(K,J) = TEMP
  218:                           END IF
  219:    40                 CONTINUE
  220:    50             CONTINUE
  221:               ELSE
  222:                   DO 80 J = 1,N
  223:                       DO 70 K = M,1,-1
  224:                           IF (B(K,J).NE.ZERO) THEN
  225:                               TEMP = ALPHA*B(K,J)
  226:                               B(K,J) = TEMP
  227:                               IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
  228:                               DO 60 I = K + 1,M
  229:                                   B(I,J) = B(I,J) + TEMP*A(I,K)
  230:    60                         CONTINUE
  231:                           END IF
  232:    70                 CONTINUE
  233:    80             CONTINUE
  234:               END IF
  235:           ELSE
  236: *
  237: *           Form  B := alpha*A**T*B.
  238: *
  239:               IF (UPPER) THEN
  240:                   DO 110 J = 1,N
  241:                       DO 100 I = M,1,-1
  242:                           TEMP = B(I,J)
  243:                           IF (NOUNIT) TEMP = TEMP*A(I,I)
  244:                           DO 90 K = 1,I - 1
  245:                               TEMP = TEMP + A(K,I)*B(K,J)
  246:    90                     CONTINUE
  247:                           B(I,J) = ALPHA*TEMP
  248:   100                 CONTINUE
  249:   110             CONTINUE
  250:               ELSE
  251:                   DO 140 J = 1,N
  252:                       DO 130 I = 1,M
  253:                           TEMP = B(I,J)
  254:                           IF (NOUNIT) TEMP = TEMP*A(I,I)
  255:                           DO 120 K = I + 1,M
  256:                               TEMP = TEMP + A(K,I)*B(K,J)
  257:   120                     CONTINUE
  258:                           B(I,J) = ALPHA*TEMP
  259:   130                 CONTINUE
  260:   140             CONTINUE
  261:               END IF
  262:           END IF
  263:       ELSE
  264:           IF (LSAME(TRANSA,'N')) THEN
  265: *
  266: *           Form  B := alpha*B*A.
  267: *
  268:               IF (UPPER) THEN
  269:                   DO 180 J = N,1,-1
  270:                       TEMP = ALPHA
  271:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  272:                       DO 150 I = 1,M
  273:                           B(I,J) = TEMP*B(I,J)
  274:   150                 CONTINUE
  275:                       DO 170 K = 1,J - 1
  276:                           IF (A(K,J).NE.ZERO) THEN
  277:                               TEMP = ALPHA*A(K,J)
  278:                               DO 160 I = 1,M
  279:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  280:   160                         CONTINUE
  281:                           END IF
  282:   170                 CONTINUE
  283:   180             CONTINUE
  284:               ELSE
  285:                   DO 220 J = 1,N
  286:                       TEMP = ALPHA
  287:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  288:                       DO 190 I = 1,M
  289:                           B(I,J) = TEMP*B(I,J)
  290:   190                 CONTINUE
  291:                       DO 210 K = J + 1,N
  292:                           IF (A(K,J).NE.ZERO) THEN
  293:                               TEMP = ALPHA*A(K,J)
  294:                               DO 200 I = 1,M
  295:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  296:   200                         CONTINUE
  297:                           END IF
  298:   210                 CONTINUE
  299:   220             CONTINUE
  300:               END IF
  301:           ELSE
  302: *
  303: *           Form  B := alpha*B*A**T.
  304: *
  305:               IF (UPPER) THEN
  306:                   DO 260 K = 1,N
  307:                       DO 240 J = 1,K - 1
  308:                           IF (A(J,K).NE.ZERO) THEN
  309:                               TEMP = ALPHA*A(J,K)
  310:                               DO 230 I = 1,M
  311:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  312:   230                         CONTINUE
  313:                           END IF
  314:   240                 CONTINUE
  315:                       TEMP = ALPHA
  316:                       IF (NOUNIT) TEMP = TEMP*A(K,K)
  317:                       IF (TEMP.NE.ONE) THEN
  318:                           DO 250 I = 1,M
  319:                               B(I,K) = TEMP*B(I,K)
  320:   250                     CONTINUE
  321:                       END IF
  322:   260             CONTINUE
  323:               ELSE
  324:                   DO 300 K = N,1,-1
  325:                       DO 280 J = K + 1,N
  326:                           IF (A(J,K).NE.ZERO) THEN
  327:                               TEMP = ALPHA*A(J,K)
  328:                               DO 270 I = 1,M
  329:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  330:   270                         CONTINUE
  331:                           END IF
  332:   280                 CONTINUE
  333:                       TEMP = ALPHA
  334:                       IF (NOUNIT) TEMP = TEMP*A(K,K)
  335:                       IF (TEMP.NE.ONE) THEN
  336:                           DO 290 I = 1,M
  337:                               B(I,K) = TEMP*B(I,K)
  338:   290                     CONTINUE
  339:                       END IF
  340:   300             CONTINUE
  341:               END IF
  342:           END IF
  343:       END IF
  344: *
  345:       RETURN
  346: *
  347: *     End of DTRMM .
  348: *
  349:       END

CVSweb interface <joel.bertrand@systella.fr>