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

CVSweb interface <joel.bertrand@systella.fr>