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

CVSweb interface <joel.bertrand@systella.fr>