File:  [local] / rpl / lapack / blas / dspr2.f
Revision 1.16: download - view: text, annotated - select for diffs - revision graph
Mon Aug 7 08:38:43 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 DSPR2
    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 DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
   12: *
   13: *       .. Scalar Arguments ..
   14: *       DOUBLE PRECISION ALPHA
   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: *> DSPR2  performs the symmetric rank 2 operation
   29: *>
   30: *>    A := alpha*x*y**T + alpha*y*x**T + A,
   31: *>
   32: *> where alpha is a scalar, x and y are n element vectors and A is an
   33: *> 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] X
   67: *> \verbatim
   68: *>          X is DOUBLE PRECISION array, dimension at least
   69: *>           ( 1 + ( n - 1 )*abs( INCX ) ).
   70: *>           Before entry, the incremented array X must contain the n
   71: *>           element vector x.
   72: *> \endverbatim
   73: *>
   74: *> \param[in] INCX
   75: *> \verbatim
   76: *>          INCX is INTEGER
   77: *>           On entry, INCX specifies the increment for the elements of
   78: *>           X. INCX must not be zero.
   79: *> \endverbatim
   80: *>
   81: *> \param[in] Y
   82: *> \verbatim
   83: *>          Y is DOUBLE PRECISION array, dimension at least
   84: *>           ( 1 + ( n - 1 )*abs( INCY ) ).
   85: *>           Before entry, the incremented array Y must contain the n
   86: *>           element vector y.
   87: *> \endverbatim
   88: *>
   89: *> \param[in] INCY
   90: *> \verbatim
   91: *>          INCY is INTEGER
   92: *>           On entry, INCY specifies the increment for the elements of
   93: *>           Y. INCY must not be zero.
   94: *> \endverbatim
   95: *>
   96: *> \param[in,out] AP
   97: *> \verbatim
   98: *>          AP is DOUBLE PRECISION array, dimension at least
   99: *>           ( ( n*( n + 1 ) )/2 ).
  100: *>           Before entry with  UPLO = 'U' or 'u', the array AP must
  101: *>           contain the upper triangular part of the symmetric matrix
  102: *>           packed sequentially, column by column, so that AP( 1 )
  103: *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 1, 2 )
  104: *>           and a( 2, 2 ) respectively, and so on. On exit, the array
  105: *>           AP is overwritten by the upper triangular part of the
  106: *>           updated matrix.
  107: *>           Before entry with UPLO = 'L' or 'l', the array AP must
  108: *>           contain the lower triangular part of the symmetric matrix
  109: *>           packed sequentially, column by column, so that AP( 1 )
  110: *>           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  111: *>           and a( 3, 1 ) respectively, and so on. On exit, the array
  112: *>           AP is overwritten by the lower triangular part of the
  113: *>           updated matrix.
  114: *> \endverbatim
  115: *
  116: *  Authors:
  117: *  ========
  118: *
  119: *> \author Univ. of Tennessee
  120: *> \author Univ. of California Berkeley
  121: *> \author Univ. of Colorado Denver
  122: *> \author NAG Ltd.
  123: *
  124: *> \ingroup double_blas_level2
  125: *
  126: *> \par Further Details:
  127: *  =====================
  128: *>
  129: *> \verbatim
  130: *>
  131: *>  Level 2 Blas routine.
  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 DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
  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:       DOUBLE PRECISION ALPHA
  149:       INTEGER INCX,INCY,N
  150:       CHARACTER UPLO
  151: *     ..
  152: *     .. Array Arguments ..
  153:       DOUBLE PRECISION AP(*),X(*),Y(*)
  154: *     ..
  155: *
  156: *  =====================================================================
  157: *
  158: *     .. Parameters ..
  159:       DOUBLE PRECISION ZERO
  160:       PARAMETER (ZERO=0.0D+0)
  161: *     ..
  162: *     .. Local Scalars ..
  163:       DOUBLE PRECISION TEMP1,TEMP2
  164:       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
  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 (N.LT.0) THEN
  180:           INFO = 2
  181:       ELSE IF (INCX.EQ.0) THEN
  182:           INFO = 5
  183:       ELSE IF (INCY.EQ.0) THEN
  184:           INFO = 7
  185:       END IF
  186:       IF (INFO.NE.0) THEN
  187:           CALL XERBLA('DSPR2 ',INFO)
  188:           RETURN
  189:       END IF
  190: *
  191: *     Quick return if possible.
  192: *
  193:       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
  194: *
  195: *     Set up the start points in X and Y if the increments are not both
  196: *     unity.
  197: *
  198:       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
  199:           IF (INCX.GT.0) THEN
  200:               KX = 1
  201:           ELSE
  202:               KX = 1 - (N-1)*INCX
  203:           END IF
  204:           IF (INCY.GT.0) THEN
  205:               KY = 1
  206:           ELSE
  207:               KY = 1 - (N-1)*INCY
  208:           END IF
  209:           JX = KX
  210:           JY = KY
  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:       KK = 1
  217:       IF (LSAME(UPLO,'U')) THEN
  218: *
  219: *        Form  A  when upper triangle is stored in AP.
  220: *
  221:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  222:               DO 20 J = 1,N
  223:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  224:                       TEMP1 = ALPHA*Y(J)
  225:                       TEMP2 = ALPHA*X(J)
  226:                       K = KK
  227:                       DO 10 I = 1,J
  228:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  229:                           K = K + 1
  230:    10                 CONTINUE
  231:                   END IF
  232:                   KK = KK + J
  233:    20         CONTINUE
  234:           ELSE
  235:               DO 40 J = 1,N
  236:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  237:                       TEMP1 = ALPHA*Y(JY)
  238:                       TEMP2 = ALPHA*X(JX)
  239:                       IX = KX
  240:                       IY = KY
  241:                       DO 30 K = KK,KK + J - 1
  242:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  243:                           IX = IX + INCX
  244:                           IY = IY + INCY
  245:    30                 CONTINUE
  246:                   END IF
  247:                   JX = JX + INCX
  248:                   JY = JY + INCY
  249:                   KK = KK + J
  250:    40         CONTINUE
  251:           END IF
  252:       ELSE
  253: *
  254: *        Form  A  when lower triangle is stored in AP.
  255: *
  256:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  257:               DO 60 J = 1,N
  258:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  259:                       TEMP1 = ALPHA*Y(J)
  260:                       TEMP2 = ALPHA*X(J)
  261:                       K = KK
  262:                       DO 50 I = J,N
  263:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  264:                           K = K + 1
  265:    50                 CONTINUE
  266:                   END IF
  267:                   KK = KK + N - J + 1
  268:    60         CONTINUE
  269:           ELSE
  270:               DO 80 J = 1,N
  271:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  272:                       TEMP1 = ALPHA*Y(JY)
  273:                       TEMP2 = ALPHA*X(JX)
  274:                       IX = JX
  275:                       IY = JY
  276:                       DO 70 K = KK,KK + N - J
  277:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  278:                           IX = IX + INCX
  279:                           IY = IY + INCY
  280:    70                 CONTINUE
  281:                   END IF
  282:                   JX = JX + INCX
  283:                   JY = JY + INCY
  284:                   KK = KK + N - J + 1
  285:    80         CONTINUE
  286:           END IF
  287:       END IF
  288: *
  289:       RETURN
  290: *
  291: *     End of DSPR2
  292: *
  293:       END

CVSweb interface <joel.bertrand@systella.fr>