File:  [local] / rpl / lapack / blas / dtrmv.f
Revision 1.3: download - view: text, annotated - select for diffs - revision graph
Fri Aug 6 15:32:20 2010 UTC (13 years, 9 months ago) by bertrand
Branches: MAIN
CVS tags: HEAD
Cohérence

    1:       SUBROUTINE DTRMV(UPLO,TRANS,DIAG,N,A,LDA,X,INCX)
    2: *     .. Scalar Arguments ..
    3:       INTEGER INCX,LDA,N
    4:       CHARACTER DIAG,TRANS,UPLO
    5: *     ..
    6: *     .. Array Arguments ..
    7:       DOUBLE PRECISION A(LDA,*),X(*)
    8: *     ..
    9: *
   10: *  Purpose
   11: *  =======
   12: *
   13: *  DTRMV  performs one of the matrix-vector operations
   14: *
   15: *     x := A*x,   or   x := A'*x,
   16: *
   17: *  where x is an n element vector and  A is an n by n unit, or non-unit,
   18: *  upper or lower triangular matrix.
   19: *
   20: *  Arguments
   21: *  ==========
   22: *
   23: *  UPLO   - CHARACTER*1.
   24: *           On entry, UPLO specifies whether the matrix is an upper or
   25: *           lower triangular matrix as follows:
   26: *
   27: *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
   28: *
   29: *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
   30: *
   31: *           Unchanged on exit.
   32: *
   33: *  TRANS  - CHARACTER*1.
   34: *           On entry, TRANS specifies the operation to be performed as
   35: *           follows:
   36: *
   37: *              TRANS = 'N' or 'n'   x := A*x.
   38: *
   39: *              TRANS = 'T' or 't'   x := A'*x.
   40: *
   41: *              TRANS = 'C' or 'c'   x := A'*x.
   42: *
   43: *           Unchanged on exit.
   44: *
   45: *  DIAG   - CHARACTER*1.
   46: *           On entry, DIAG specifies whether or not A is unit
   47: *           triangular as follows:
   48: *
   49: *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
   50: *
   51: *              DIAG = 'N' or 'n'   A is not assumed to be unit
   52: *                                  triangular.
   53: *
   54: *           Unchanged on exit.
   55: *
   56: *  N      - INTEGER.
   57: *           On entry, N specifies the order of the matrix A.
   58: *           N must be at least zero.
   59: *           Unchanged on exit.
   60: *
   61: *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
   62: *           Before entry with  UPLO = 'U' or 'u', the leading n by n
   63: *           upper triangular part of the array A must contain the upper
   64: *           triangular matrix and the strictly lower triangular part of
   65: *           A is not referenced.
   66: *           Before entry with UPLO = 'L' or 'l', the leading n by n
   67: *           lower triangular part of the array A must contain the lower
   68: *           triangular matrix and the strictly upper triangular part of
   69: *           A is not referenced.
   70: *           Note that when  DIAG = 'U' or 'u', the diagonal elements of
   71: *           A are not referenced either, but are assumed to be unity.
   72: *           Unchanged on exit.
   73: *
   74: *  LDA    - INTEGER.
   75: *           On entry, LDA specifies the first dimension of A as declared
   76: *           in the calling (sub) program. LDA must be at least
   77: *           max( 1, n ).
   78: *           Unchanged on exit.
   79: *
   80: *  X      - DOUBLE PRECISION array of dimension at least
   81: *           ( 1 + ( n - 1 )*abs( INCX ) ).
   82: *           Before entry, the incremented array X must contain the n
   83: *           element vector x. On exit, X is overwritten with the
   84: *           tranformed vector x.
   85: *
   86: *  INCX   - INTEGER.
   87: *           On entry, INCX specifies the increment for the elements of
   88: *           X. INCX must not be zero.
   89: *           Unchanged on exit.
   90: *
   91: *  Further Details
   92: *  ===============
   93: *
   94: *  Level 2 Blas routine.
   95: *
   96: *  -- Written on 22-October-1986.
   97: *     Jack Dongarra, Argonne National Lab.
   98: *     Jeremy Du Croz, Nag Central Office.
   99: *     Sven Hammarling, Nag Central Office.
  100: *     Richard Hanson, Sandia National Labs.
  101: *
  102: *  =====================================================================
  103: *
  104: *     .. Parameters ..
  105:       DOUBLE PRECISION ZERO
  106:       PARAMETER (ZERO=0.0D+0)
  107: *     ..
  108: *     .. Local Scalars ..
  109:       DOUBLE PRECISION TEMP
  110:       INTEGER I,INFO,IX,J,JX,KX
  111:       LOGICAL NOUNIT
  112: *     ..
  113: *     .. External Functions ..
  114:       LOGICAL LSAME
  115:       EXTERNAL LSAME
  116: *     ..
  117: *     .. External Subroutines ..
  118:       EXTERNAL XERBLA
  119: *     ..
  120: *     .. Intrinsic Functions ..
  121:       INTRINSIC MAX
  122: *     ..
  123: *
  124: *     Test the input parameters.
  125: *
  126:       INFO = 0
  127:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  128:           INFO = 1
  129:       ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
  130:      +         .NOT.LSAME(TRANS,'C')) THEN
  131:           INFO = 2
  132:       ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
  133:           INFO = 3
  134:       ELSE IF (N.LT.0) THEN
  135:           INFO = 4
  136:       ELSE IF (LDA.LT.MAX(1,N)) THEN
  137:           INFO = 6
  138:       ELSE IF (INCX.EQ.0) THEN
  139:           INFO = 8
  140:       END IF
  141:       IF (INFO.NE.0) THEN
  142:           CALL XERBLA('DTRMV ',INFO)
  143:           RETURN
  144:       END IF
  145: *
  146: *     Quick return if possible.
  147: *
  148:       IF (N.EQ.0) RETURN
  149: *
  150:       NOUNIT = LSAME(DIAG,'N')
  151: *
  152: *     Set up the start point in X if the increment is not unity. This
  153: *     will be  ( N - 1 )*INCX  too small for descending loops.
  154: *
  155:       IF (INCX.LE.0) THEN
  156:           KX = 1 - (N-1)*INCX
  157:       ELSE IF (INCX.NE.1) THEN
  158:           KX = 1
  159:       END IF
  160: *
  161: *     Start the operations. In this version the elements of A are
  162: *     accessed sequentially with one pass through A.
  163: *
  164:       IF (LSAME(TRANS,'N')) THEN
  165: *
  166: *        Form  x := A*x.
  167: *
  168:           IF (LSAME(UPLO,'U')) THEN
  169:               IF (INCX.EQ.1) THEN
  170:                   DO 20 J = 1,N
  171:                       IF (X(J).NE.ZERO) THEN
  172:                           TEMP = X(J)
  173:                           DO 10 I = 1,J - 1
  174:                               X(I) = X(I) + TEMP*A(I,J)
  175:    10                     CONTINUE
  176:                           IF (NOUNIT) X(J) = X(J)*A(J,J)
  177:                       END IF
  178:    20             CONTINUE
  179:               ELSE
  180:                   JX = KX
  181:                   DO 40 J = 1,N
  182:                       IF (X(JX).NE.ZERO) THEN
  183:                           TEMP = X(JX)
  184:                           IX = KX
  185:                           DO 30 I = 1,J - 1
  186:                               X(IX) = X(IX) + TEMP*A(I,J)
  187:                               IX = IX + INCX
  188:    30                     CONTINUE
  189:                           IF (NOUNIT) X(JX) = X(JX)*A(J,J)
  190:                       END IF
  191:                       JX = JX + INCX
  192:    40             CONTINUE
  193:               END IF
  194:           ELSE
  195:               IF (INCX.EQ.1) THEN
  196:                   DO 60 J = N,1,-1
  197:                       IF (X(J).NE.ZERO) THEN
  198:                           TEMP = X(J)
  199:                           DO 50 I = N,J + 1,-1
  200:                               X(I) = X(I) + TEMP*A(I,J)
  201:    50                     CONTINUE
  202:                           IF (NOUNIT) X(J) = X(J)*A(J,J)
  203:                       END IF
  204:    60             CONTINUE
  205:               ELSE
  206:                   KX = KX + (N-1)*INCX
  207:                   JX = KX
  208:                   DO 80 J = N,1,-1
  209:                       IF (X(JX).NE.ZERO) THEN
  210:                           TEMP = X(JX)
  211:                           IX = KX
  212:                           DO 70 I = N,J + 1,-1
  213:                               X(IX) = X(IX) + TEMP*A(I,J)
  214:                               IX = IX - INCX
  215:    70                     CONTINUE
  216:                           IF (NOUNIT) X(JX) = X(JX)*A(J,J)
  217:                       END IF
  218:                       JX = JX - INCX
  219:    80             CONTINUE
  220:               END IF
  221:           END IF
  222:       ELSE
  223: *
  224: *        Form  x := A'*x.
  225: *
  226:           IF (LSAME(UPLO,'U')) THEN
  227:               IF (INCX.EQ.1) THEN
  228:                   DO 100 J = N,1,-1
  229:                       TEMP = X(J)
  230:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  231:                       DO 90 I = J - 1,1,-1
  232:                           TEMP = TEMP + A(I,J)*X(I)
  233:    90                 CONTINUE
  234:                       X(J) = TEMP
  235:   100             CONTINUE
  236:               ELSE
  237:                   JX = KX + (N-1)*INCX
  238:                   DO 120 J = N,1,-1
  239:                       TEMP = X(JX)
  240:                       IX = JX
  241:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  242:                       DO 110 I = J - 1,1,-1
  243:                           IX = IX - INCX
  244:                           TEMP = TEMP + A(I,J)*X(IX)
  245:   110                 CONTINUE
  246:                       X(JX) = TEMP
  247:                       JX = JX - INCX
  248:   120             CONTINUE
  249:               END IF
  250:           ELSE
  251:               IF (INCX.EQ.1) THEN
  252:                   DO 140 J = 1,N
  253:                       TEMP = X(J)
  254:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  255:                       DO 130 I = J + 1,N
  256:                           TEMP = TEMP + A(I,J)*X(I)
  257:   130                 CONTINUE
  258:                       X(J) = TEMP
  259:   140             CONTINUE
  260:               ELSE
  261:                   JX = KX
  262:                   DO 160 J = 1,N
  263:                       TEMP = X(JX)
  264:                       IX = JX
  265:                       IF (NOUNIT) TEMP = TEMP*A(J,J)
  266:                       DO 150 I = J + 1,N
  267:                           IX = IX + INCX
  268:                           TEMP = TEMP + A(I,J)*X(IX)
  269:   150                 CONTINUE
  270:                       X(JX) = TEMP
  271:                       JX = JX + INCX
  272:   160             CONTINUE
  273:               END IF
  274:           END IF
  275:       END IF
  276: *
  277:       RETURN
  278: *
  279: *     End of DTRMV .
  280: *
  281:       END

CVSweb interface <joel.bertrand@systella.fr>