File:  [local] / rpl / lapack / blas / dtpmv.f
Revision 1.16: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:44 2023 UTC (8 months, 3 weeks ago) by bertrand
Branches: MAIN
CVS tags: rpl-4_1_35, rpl-4_1_34, HEAD
Première mise à jour de lapack et blas.

    1: *> \brief \b DTPMV
    2: *
    3: *  =========== DOCUMENTATION ===========
    4: *
    5: * Online html documentation available at
    6: *            http://www.netlib.org/lapack/explore-html/
    7: *
    8: *  Definition:
    9: *  ===========
   10: *
   11: *       SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
   12: *
   13: *       .. Scalar Arguments ..
   14: *       INTEGER INCX,N
   15: *       CHARACTER DIAG,TRANS,UPLO
   16: *       ..
   17: *       .. Array Arguments ..
   18: *       DOUBLE PRECISION AP(*),X(*)
   19: *       ..
   20: *
   21: *
   22: *> \par Purpose:
   23: *  =============
   24: *>
   25: *> \verbatim
   26: *>
   27: *> DTPMV  performs one of the matrix-vector operations
   28: *>
   29: *>    x := A*x,   or   x := A**T*x,
   30: *>
   31: *> where x is an n element vector and  A is an n by n unit, or non-unit,
   32: *> upper or lower triangular matrix, supplied in packed form.
   33: *> \endverbatim
   34: *
   35: *  Arguments:
   36: *  ==========
   37: *
   38: *> \param[in] UPLO
   39: *> \verbatim
   40: *>          UPLO is CHARACTER*1
   41: *>           On entry, UPLO specifies whether the matrix is an upper or
   42: *>           lower triangular matrix as follows:
   43: *>
   44: *>              UPLO = 'U' or 'u'   A is an upper triangular matrix.
   45: *>
   46: *>              UPLO = 'L' or 'l'   A is a lower triangular matrix.
   47: *> \endverbatim
   48: *>
   49: *> \param[in] TRANS
   50: *> \verbatim
   51: *>          TRANS is CHARACTER*1
   52: *>           On entry, TRANS specifies the operation to be performed as
   53: *>           follows:
   54: *>
   55: *>              TRANS = 'N' or 'n'   x := A*x.
   56: *>
   57: *>              TRANS = 'T' or 't'   x := A**T*x.
   58: *>
   59: *>              TRANS = 'C' or 'c'   x := A**T*x.
   60: *> \endverbatim
   61: *>
   62: *> \param[in] DIAG
   63: *> \verbatim
   64: *>          DIAG is CHARACTER*1
   65: *>           On entry, DIAG specifies whether or not A is unit
   66: *>           triangular as follows:
   67: *>
   68: *>              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
   69: *>
   70: *>              DIAG = 'N' or 'n'   A is not assumed to be unit
   71: *>                                  triangular.
   72: *> \endverbatim
   73: *>
   74: *> \param[in] N
   75: *> \verbatim
   76: *>          N is INTEGER
   77: *>           On entry, N specifies the order of the matrix A.
   78: *>           N must be at least zero.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] AP
   82: *> \verbatim
   83: *>          AP is DOUBLE PRECISION array, dimension at least
   84: *>           ( ( n*( n + 1 ) )/2 ).
   85: *>           Before entry with  UPLO = 'U' or 'u', the array AP must
   86: *>           contain the upper triangular matrix packed sequentially,
   87: *>           column by column, so that AP( 1 ) contains a( 1, 1 ),
   88: *>           AP( 2 ) and AP( 3 ) contain a( 1, 2 ) and a( 2, 2 )
   89: *>           respectively, and so on.
   90: *>           Before entry with UPLO = 'L' or 'l', the array AP must
   91: *>           contain the lower triangular matrix packed sequentially,
   92: *>           column by column, so that AP( 1 ) contains a( 1, 1 ),
   93: *>           AP( 2 ) and AP( 3 ) contain a( 2, 1 ) and a( 3, 1 )
   94: *>           respectively, and so on.
   95: *>           Note that when  DIAG = 'U' or 'u', the diagonal elements of
   96: *>           A are not referenced, but are assumed to be unity.
   97: *> \endverbatim
   98: *>
   99: *> \param[in,out] X
  100: *> \verbatim
  101: *>          X is DOUBLE PRECISION array, dimension at least
  102: *>           ( 1 + ( n - 1 )*abs( INCX ) ).
  103: *>           Before entry, the incremented array X must contain the n
  104: *>           element vector x. On exit, X is overwritten with the
  105: *>           transformed vector x.
  106: *> \endverbatim
  107: *>
  108: *> \param[in] INCX
  109: *> \verbatim
  110: *>          INCX is INTEGER
  111: *>           On entry, INCX specifies the increment for the elements of
  112: *>           X. INCX must not be zero.
  113: *> \endverbatim
  114: *
  115: *  Authors:
  116: *  ========
  117: *
  118: *> \author Univ. of Tennessee
  119: *> \author Univ. of California Berkeley
  120: *> \author Univ. of Colorado Denver
  121: *> \author NAG Ltd.
  122: *
  123: *> \ingroup double_blas_level2
  124: *
  125: *> \par Further Details:
  126: *  =====================
  127: *>
  128: *> \verbatim
  129: *>
  130: *>  Level 2 Blas routine.
  131: *>  The vector and matrix arguments are not referenced when N = 0, or M = 0
  132: *>
  133: *>  -- Written on 22-October-1986.
  134: *>     Jack Dongarra, Argonne National Lab.
  135: *>     Jeremy Du Croz, Nag Central Office.
  136: *>     Sven Hammarling, Nag Central Office.
  137: *>     Richard Hanson, Sandia National Labs.
  138: *> \endverbatim
  139: *>
  140: *  =====================================================================
  141:       SUBROUTINE DTPMV(UPLO,TRANS,DIAG,N,AP,X,INCX)
  142: *
  143: *  -- Reference BLAS level2 routine --
  144: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
  145: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  146: *
  147: *     .. Scalar Arguments ..
  148:       INTEGER INCX,N
  149:       CHARACTER DIAG,TRANS,UPLO
  150: *     ..
  151: *     .. Array Arguments ..
  152:       DOUBLE PRECISION AP(*),X(*)
  153: *     ..
  154: *
  155: *  =====================================================================
  156: *
  157: *     .. Parameters ..
  158:       DOUBLE PRECISION ZERO
  159:       PARAMETER (ZERO=0.0D+0)
  160: *     ..
  161: *     .. Local Scalars ..
  162:       DOUBLE PRECISION TEMP
  163:       INTEGER I,INFO,IX,J,JX,K,KK,KX
  164:       LOGICAL NOUNIT
  165: *     ..
  166: *     .. External Functions ..
  167:       LOGICAL LSAME
  168:       EXTERNAL LSAME
  169: *     ..
  170: *     .. External Subroutines ..
  171:       EXTERNAL XERBLA
  172: *     ..
  173: *
  174: *     Test the input parameters.
  175: *
  176:       INFO = 0
  177:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  178:           INFO = 1
  179:       ELSE IF (.NOT.LSAME(TRANS,'N') .AND. .NOT.LSAME(TRANS,'T') .AND.
  180:      +         .NOT.LSAME(TRANS,'C')) THEN
  181:           INFO = 2
  182:       ELSE IF (.NOT.LSAME(DIAG,'U') .AND. .NOT.LSAME(DIAG,'N')) THEN
  183:           INFO = 3
  184:       ELSE IF (N.LT.0) THEN
  185:           INFO = 4
  186:       ELSE IF (INCX.EQ.0) THEN
  187:           INFO = 7
  188:       END IF
  189:       IF (INFO.NE.0) THEN
  190:           CALL XERBLA('DTPMV ',INFO)
  191:           RETURN
  192:       END IF
  193: *
  194: *     Quick return if possible.
  195: *
  196:       IF (N.EQ.0) RETURN
  197: *
  198:       NOUNIT = LSAME(DIAG,'N')
  199: *
  200: *     Set up the start point in X if the increment is not unity. This
  201: *     will be  ( N - 1 )*INCX  too small for descending loops.
  202: *
  203:       IF (INCX.LE.0) THEN
  204:           KX = 1 - (N-1)*INCX
  205:       ELSE IF (INCX.NE.1) THEN
  206:           KX = 1
  207:       END IF
  208: *
  209: *     Start the operations. In this version the elements of AP are
  210: *     accessed sequentially with one pass through AP.
  211: *
  212:       IF (LSAME(TRANS,'N')) THEN
  213: *
  214: *        Form  x:= A*x.
  215: *
  216:           IF (LSAME(UPLO,'U')) THEN
  217:               KK = 1
  218:               IF (INCX.EQ.1) THEN
  219:                   DO 20 J = 1,N
  220:                       IF (X(J).NE.ZERO) THEN
  221:                           TEMP = X(J)
  222:                           K = KK
  223:                           DO 10 I = 1,J - 1
  224:                               X(I) = X(I) + TEMP*AP(K)
  225:                               K = K + 1
  226:    10                     CONTINUE
  227:                           IF (NOUNIT) X(J) = X(J)*AP(KK+J-1)
  228:                       END IF
  229:                       KK = KK + J
  230:    20             CONTINUE
  231:               ELSE
  232:                   JX = KX
  233:                   DO 40 J = 1,N
  234:                       IF (X(JX).NE.ZERO) THEN
  235:                           TEMP = X(JX)
  236:                           IX = KX
  237:                           DO 30 K = KK,KK + J - 2
  238:                               X(IX) = X(IX) + TEMP*AP(K)
  239:                               IX = IX + INCX
  240:    30                     CONTINUE
  241:                           IF (NOUNIT) X(JX) = X(JX)*AP(KK+J-1)
  242:                       END IF
  243:                       JX = JX + INCX
  244:                       KK = KK + J
  245:    40             CONTINUE
  246:               END IF
  247:           ELSE
  248:               KK = (N* (N+1))/2
  249:               IF (INCX.EQ.1) THEN
  250:                   DO 60 J = N,1,-1
  251:                       IF (X(J).NE.ZERO) THEN
  252:                           TEMP = X(J)
  253:                           K = KK
  254:                           DO 50 I = N,J + 1,-1
  255:                               X(I) = X(I) + TEMP*AP(K)
  256:                               K = K - 1
  257:    50                     CONTINUE
  258:                           IF (NOUNIT) X(J) = X(J)*AP(KK-N+J)
  259:                       END IF
  260:                       KK = KK - (N-J+1)
  261:    60             CONTINUE
  262:               ELSE
  263:                   KX = KX + (N-1)*INCX
  264:                   JX = KX
  265:                   DO 80 J = N,1,-1
  266:                       IF (X(JX).NE.ZERO) THEN
  267:                           TEMP = X(JX)
  268:                           IX = KX
  269:                           DO 70 K = KK,KK - (N- (J+1)),-1
  270:                               X(IX) = X(IX) + TEMP*AP(K)
  271:                               IX = IX - INCX
  272:    70                     CONTINUE
  273:                           IF (NOUNIT) X(JX) = X(JX)*AP(KK-N+J)
  274:                       END IF
  275:                       JX = JX - INCX
  276:                       KK = KK - (N-J+1)
  277:    80             CONTINUE
  278:               END IF
  279:           END IF
  280:       ELSE
  281: *
  282: *        Form  x := A**T*x.
  283: *
  284:           IF (LSAME(UPLO,'U')) THEN
  285:               KK = (N* (N+1))/2
  286:               IF (INCX.EQ.1) THEN
  287:                   DO 100 J = N,1,-1
  288:                       TEMP = X(J)
  289:                       IF (NOUNIT) TEMP = TEMP*AP(KK)
  290:                       K = KK - 1
  291:                       DO 90 I = J - 1,1,-1
  292:                           TEMP = TEMP + AP(K)*X(I)
  293:                           K = K - 1
  294:    90                 CONTINUE
  295:                       X(J) = TEMP
  296:                       KK = KK - J
  297:   100             CONTINUE
  298:               ELSE
  299:                   JX = KX + (N-1)*INCX
  300:                   DO 120 J = N,1,-1
  301:                       TEMP = X(JX)
  302:                       IX = JX
  303:                       IF (NOUNIT) TEMP = TEMP*AP(KK)
  304:                       DO 110 K = KK - 1,KK - J + 1,-1
  305:                           IX = IX - INCX
  306:                           TEMP = TEMP + AP(K)*X(IX)
  307:   110                 CONTINUE
  308:                       X(JX) = TEMP
  309:                       JX = JX - INCX
  310:                       KK = KK - J
  311:   120             CONTINUE
  312:               END IF
  313:           ELSE
  314:               KK = 1
  315:               IF (INCX.EQ.1) THEN
  316:                   DO 140 J = 1,N
  317:                       TEMP = X(J)
  318:                       IF (NOUNIT) TEMP = TEMP*AP(KK)
  319:                       K = KK + 1
  320:                       DO 130 I = J + 1,N
  321:                           TEMP = TEMP + AP(K)*X(I)
  322:                           K = K + 1
  323:   130                 CONTINUE
  324:                       X(J) = TEMP
  325:                       KK = KK + (N-J+1)
  326:   140             CONTINUE
  327:               ELSE
  328:                   JX = KX
  329:                   DO 160 J = 1,N
  330:                       TEMP = X(JX)
  331:                       IX = JX
  332:                       IF (NOUNIT) TEMP = TEMP*AP(KK)
  333:                       DO 150 K = KK + 1,KK + N - J
  334:                           IX = IX + INCX
  335:                           TEMP = TEMP + AP(K)*X(IX)
  336:   150                 CONTINUE
  337:                       X(JX) = TEMP
  338:                       JX = JX + INCX
  339:                       KK = KK + (N-J+1)
  340:   160             CONTINUE
  341:               END IF
  342:           END IF
  343:       END IF
  344: *
  345:       RETURN
  346: *
  347: *     End of DTPMV
  348: *
  349:       END

CVSweb interface <joel.bertrand@systella.fr>