File:  [local] / rpl / lapack / blas / ztrmm.f
Revision 1.5: download - view: text, annotated - select for diffs - revision graph
Fri Aug 13 21:03:42 2010 UTC (13 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_0_19, rpl-4_0_18, HEAD
Patches pour OS/2

    1:       SUBROUTINE ZTRMM(SIDE,UPLO,TRANSA,DIAG,M,N,ALPHA,A,LDA,B,LDB)
    2: *     .. Scalar Arguments ..
    3:       DOUBLE COMPLEX ALPHA
    4:       INTEGER LDA,LDB,M,N
    5:       CHARACTER DIAG,SIDE,TRANSA,UPLO
    6: *     ..
    7: *     .. Array Arguments ..
    8:       DOUBLE COMPLEX A(LDA,*),B(LDB,*)
    9: *     ..
   10: *
   11: *  Purpose
   12: *  =======
   13: *
   14: *  ZTRMM  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'   or   op( A ) = conjg( A' ).
   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'.
   53: *
   54: *              TRANSA = 'C' or 'c'   op( A ) = conjg( A' ).
   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  - COMPLEX*16      .
   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      - COMPLEX*16       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      - COMPLEX*16       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 DCONJG,MAX
  139: *     ..
  140: *     .. Local Scalars ..
  141:       DOUBLE COMPLEX TEMP
  142:       INTEGER I,INFO,J,K,NROWA
  143:       LOGICAL LSIDE,NOCONJ,NOUNIT,UPPER
  144: *     ..
  145: *     .. Parameters ..
  146:       DOUBLE COMPLEX ONE
  147:       PARAMETER (ONE= (1.0D+0,0.0D+0))
  148:       DOUBLE COMPLEX ZERO
  149:       PARAMETER (ZERO= (0.0D+0,0.0D+0))
  150: *     ..
  151: *
  152: *     Test the input parameters.
  153: *
  154:       LSIDE = LSAME(SIDE,'L')
  155:       IF (LSIDE) THEN
  156:           NROWA = M
  157:       ELSE
  158:           NROWA = N
  159:       END IF
  160:       NOCONJ = LSAME(TRANSA,'T')
  161:       NOUNIT = LSAME(DIAG,'N')
  162:       UPPER = LSAME(UPLO,'U')
  163: *
  164:       INFO = 0
  165:       IF ((.NOT.LSIDE) .AND. (.NOT.LSAME(SIDE,'R'))) THEN
  166:           INFO = 1
  167:       ELSE IF ((.NOT.UPPER) .AND. (.NOT.LSAME(UPLO,'L'))) THEN
  168:           INFO = 2
  169:       ELSE IF ((.NOT.LSAME(TRANSA,'N')) .AND.
  170:      +         (.NOT.LSAME(TRANSA,'T')) .AND.
  171:      +         (.NOT.LSAME(TRANSA,'C'))) THEN
  172:           INFO = 3
  173:       ELSE IF ((.NOT.LSAME(DIAG,'U')) .AND. (.NOT.LSAME(DIAG,'N'))) THEN
  174:           INFO = 4
  175:       ELSE IF (M.LT.0) THEN
  176:           INFO = 5
  177:       ELSE IF (N.LT.0) THEN
  178:           INFO = 6
  179:       ELSE IF (LDA.LT.MAX(1,NROWA)) THEN
  180:           INFO = 9
  181:       ELSE IF (LDB.LT.MAX(1,M)) THEN
  182:           INFO = 11
  183:       END IF
  184:       IF (INFO.NE.0) THEN
  185:           CALL XERBLA('ZTRMM ',INFO)
  186:           RETURN
  187:       END IF
  188: *
  189: *     Quick return if possible.
  190: *
  191:       IF (M.EQ.0 .OR. N.EQ.0) RETURN
  192: *
  193: *     And when  alpha.eq.zero.
  194: *
  195:       IF (ALPHA.EQ.ZERO) THEN
  196:           DO 20 J = 1,N
  197:               DO 10 I = 1,M
  198:                   B(I,J) = ZERO
  199:    10         CONTINUE
  200:    20     CONTINUE
  201:           RETURN
  202:       END IF
  203: *
  204: *     Start the operations.
  205: *
  206:       IF (LSIDE) THEN
  207:           IF (LSAME(TRANSA,'N')) THEN
  208: *
  209: *           Form  B := alpha*A*B.
  210: *
  211:               IF (UPPER) THEN
  212:                   DO 50 J = 1,N
  213:                       DO 40 K = 1,M
  214:                           IF (B(K,J).NE.ZERO) THEN
  215:                               TEMP = ALPHA*B(K,J)
  216:                               DO 30 I = 1,K - 1
  217:                                   B(I,J) = B(I,J) + TEMP*A(I,K)
  218:    30                         CONTINUE
  219:                               IF (NOUNIT) TEMP = TEMP*A(K,K)
  220:                               B(K,J) = TEMP
  221:                           END IF
  222:    40                 CONTINUE
  223:    50             CONTINUE
  224:               ELSE
  225:                   DO 80 J = 1,N
  226:                       DO 70 K = M,1,-1
  227:                           IF (B(K,J).NE.ZERO) THEN
  228:                               TEMP = ALPHA*B(K,J)
  229:                               B(K,J) = TEMP
  230:                               IF (NOUNIT) B(K,J) = B(K,J)*A(K,K)
  231:                               DO 60 I = K + 1,M
  232:                                   B(I,J) = B(I,J) + TEMP*A(I,K)
  233:    60                         CONTINUE
  234:                           END IF
  235:    70                 CONTINUE
  236:    80             CONTINUE
  237:               END IF
  238:           ELSE
  239: *
  240: *           Form  B := alpha*A'*B   or   B := alpha*conjg( A' )*B.
  241: *
  242:               IF (UPPER) THEN
  243:                   DO 120 J = 1,N
  244:                       DO 110 I = M,1,-1
  245:                           TEMP = B(I,J)
  246:                           IF (NOCONJ) THEN
  247:                               IF (NOUNIT) TEMP = TEMP*A(I,I)
  248:                               DO 90 K = 1,I - 1
  249:                                   TEMP = TEMP + A(K,I)*B(K,J)
  250:    90                         CONTINUE
  251:                           ELSE
  252:                               IF (NOUNIT) TEMP = TEMP*DCONJG(A(I,I))
  253:                               DO 100 K = 1,I - 1
  254:                                   TEMP = TEMP + DCONJG(A(K,I))*B(K,J)
  255:   100                         CONTINUE
  256:                           END IF
  257:                           B(I,J) = ALPHA*TEMP
  258:   110                 CONTINUE
  259:   120             CONTINUE
  260:               ELSE
  261:                   DO 160 J = 1,N
  262:                       DO 150 I = 1,M
  263:                           TEMP = B(I,J)
  264:                           IF (NOCONJ) THEN
  265:                               IF (NOUNIT) TEMP = TEMP*A(I,I)
  266:                               DO 130 K = I + 1,M
  267:                                   TEMP = TEMP + A(K,I)*B(K,J)
  268:   130                         CONTINUE
  269:                           ELSE
  270:                               IF (NOUNIT) TEMP = TEMP*DCONJG(A(I,I))
  271:                               DO 140 K = I + 1,M
  272:                                   TEMP = TEMP + DCONJG(A(K,I))*B(K,J)
  273:   140                         CONTINUE
  274:                           END IF
  275:                           B(I,J) = ALPHA*TEMP
  276:   150                 CONTINUE
  277:   160             CONTINUE
  278:               END IF
  279:           END IF
  280:       ELSE
  281:           IF (LSAME(TRANSA,'N')) THEN
  282: *
  283: *           Form  B := alpha*B*A.
  284: *
  285:               IF (UPPER) THEN
  286:                   DO 200 J = N,1,-1
  287:                       TEMP = ALPHA
  288:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  289:                       DO 170 I = 1,M
  290:                           B(I,J) = TEMP*B(I,J)
  291:   170                 CONTINUE
  292:                       DO 190 K = 1,J - 1
  293:                           IF (A(K,J).NE.ZERO) THEN
  294:                               TEMP = ALPHA*A(K,J)
  295:                               DO 180 I = 1,M
  296:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  297:   180                         CONTINUE
  298:                           END IF
  299:   190                 CONTINUE
  300:   200             CONTINUE
  301:               ELSE
  302:                   DO 240 J = 1,N
  303:                       TEMP = ALPHA
  304:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  305:                       DO 210 I = 1,M
  306:                           B(I,J) = TEMP*B(I,J)
  307:   210                 CONTINUE
  308:                       DO 230 K = J + 1,N
  309:                           IF (A(K,J).NE.ZERO) THEN
  310:                               TEMP = ALPHA*A(K,J)
  311:                               DO 220 I = 1,M
  312:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  313:   220                         CONTINUE
  314:                           END IF
  315:   230                 CONTINUE
  316:   240             CONTINUE
  317:               END IF
  318:           ELSE
  319: *
  320: *           Form  B := alpha*B*A'   or   B := alpha*B*conjg( A' ).
  321: *
  322:               IF (UPPER) THEN
  323:                   DO 280 K = 1,N
  324:                       DO 260 J = 1,K - 1
  325:                           IF (A(J,K).NE.ZERO) THEN
  326:                               IF (NOCONJ) THEN
  327:                                   TEMP = ALPHA*A(J,K)
  328:                               ELSE
  329:                                   TEMP = ALPHA*DCONJG(A(J,K))
  330:                               END IF
  331:                               DO 250 I = 1,M
  332:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  333:   250                         CONTINUE
  334:                           END IF
  335:   260                 CONTINUE
  336:                       TEMP = ALPHA
  337:                       IF (NOUNIT) THEN
  338:                           IF (NOCONJ) THEN
  339:                               TEMP = TEMP*A(K,K)
  340:                           ELSE
  341:                               TEMP = TEMP*DCONJG(A(K,K))
  342:                           END IF
  343:                       END IF
  344:                       IF (TEMP.NE.ONE) THEN
  345:                           DO 270 I = 1,M
  346:                               B(I,K) = TEMP*B(I,K)
  347:   270                     CONTINUE
  348:                       END IF
  349:   280             CONTINUE
  350:               ELSE
  351:                   DO 320 K = N,1,-1
  352:                       DO 300 J = K + 1,N
  353:                           IF (A(J,K).NE.ZERO) THEN
  354:                               IF (NOCONJ) THEN
  355:                                   TEMP = ALPHA*A(J,K)
  356:                               ELSE
  357:                                   TEMP = ALPHA*DCONJG(A(J,K))
  358:                               END IF
  359:                               DO 290 I = 1,M
  360:                                   B(I,J) = B(I,J) + TEMP*B(I,K)
  361:   290                         CONTINUE
  362:                           END IF
  363:   300                 CONTINUE
  364:                       TEMP = ALPHA
  365:                       IF (NOUNIT) THEN
  366:                           IF (NOCONJ) THEN
  367:                               TEMP = TEMP*A(K,K)
  368:                           ELSE
  369:                               TEMP = TEMP*DCONJG(A(K,K))
  370:                           END IF
  371:                       END IF
  372:                       IF (TEMP.NE.ONE) THEN
  373:                           DO 310 I = 1,M
  374:                               B(I,K) = TEMP*B(I,K)
  375:   310                     CONTINUE
  376:                       END IF
  377:   320             CONTINUE
  378:               END IF
  379:           END IF
  380:       END IF
  381: *
  382:       RETURN
  383: *
  384: *     End of ZTRMM .
  385: *
  386:       END

CVSweb interface <joel.bertrand@systella.fr>