Return to dspr2.f CVS log

Up to [local] / rpl / lapack / blas

Cohérence.

    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 of 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 of 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 of 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: *> \date November 2011
  125: *
  126: *> \ingroup double_blas_level2
  127: *
  128: *> \par Further Details:
  129: *  =====================
  130: *>
  131: *> \verbatim
  132: *>
  133: *>  Level 2 Blas routine.
  134: *>
  135: *>  -- Written on 22-October-1986.
  136: *>     Jack Dongarra, Argonne National Lab.
  137: *>     Jeremy Du Croz, Nag Central Office.
  138: *>     Sven Hammarling, Nag Central Office.
  139: *>     Richard Hanson, Sandia National Labs.
  140: *> \endverbatim
  141: *>
  142: *  =====================================================================
  143:       SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
  144: *
  145: *  -- Reference BLAS level2 routine (version 3.4.0) --
  146: *  -- Reference BLAS is a software package provided by Univ. of Tennessee,    --
  147: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  148: *     November 2011
  149: *
  150: *     .. Scalar Arguments ..
  151:       DOUBLE PRECISION ALPHA
  152:       INTEGER INCX,INCY,N
  153:       CHARACTER UPLO
  154: *     ..
  155: *     .. Array Arguments ..
  156:       DOUBLE PRECISION AP(*),X(*),Y(*)
  157: *     ..
  158: *
  159: *  =====================================================================
  160: *
  161: *     .. Parameters ..
  162:       DOUBLE PRECISION ZERO
  163:       PARAMETER (ZERO=0.0D+0)
  164: *     ..
  165: *     .. Local Scalars ..
  166:       DOUBLE PRECISION TEMP1,TEMP2
  167:       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
  168: *     ..
  169: *     .. External Functions ..
  170:       LOGICAL LSAME
  171:       EXTERNAL LSAME
  172: *     ..
  173: *     .. External Subroutines ..
  174:       EXTERNAL XERBLA
  175: *     ..
  176: *
  177: *     Test the input parameters.
  178: *
  179:       INFO = 0
  180:       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
  181:           INFO = 1
  182:       ELSE IF (N.LT.0) THEN
  183:           INFO = 2
  184:       ELSE IF (INCX.EQ.0) THEN
  185:           INFO = 5
  186:       ELSE IF (INCY.EQ.0) THEN
  187:           INFO = 7
  188:       END IF
  189:       IF (INFO.NE.0) THEN
  190:           CALL XERBLA('DSPR2 ',INFO)
  191:           RETURN
  192:       END IF
  193: *
  194: *     Quick return if possible.
  195: *
  196:       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
  197: *
  198: *     Set up the start points in X and Y if the increments are not both
  199: *     unity.
  200: *
  201:       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
  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:           JX = KX
  213:           JY = KY
  214:       END IF
  215: *
  216: *     Start the operations. In this version the elements of the array AP
  217: *     are accessed sequentially with one pass through AP.
  218: *
  219:       KK = 1
  220:       IF (LSAME(UPLO,'U')) THEN
  221: *
  222: *        Form  A  when upper triangle is stored in AP.
  223: *
  224:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  225:               DO 20 J = 1,N
  226:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  227:                       TEMP1 = ALPHA*Y(J)
  228:                       TEMP2 = ALPHA*X(J)
  229:                       K = KK
  230:                       DO 10 I = 1,J
  231:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  232:                           K = K + 1
  233:    10                 CONTINUE
  234:                   END IF
  235:                   KK = KK + J
  236:    20         CONTINUE
  237:           ELSE
  238:               DO 40 J = 1,N
  239:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  240:                       TEMP1 = ALPHA*Y(JY)
  241:                       TEMP2 = ALPHA*X(JX)
  242:                       IX = KX
  243:                       IY = KY
  244:                       DO 30 K = KK,KK + J - 1
  245:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  246:                           IX = IX + INCX
  247:                           IY = IY + INCY
  248:    30                 CONTINUE
  249:                   END IF
  250:                   JX = JX + INCX
  251:                   JY = JY + INCY
  252:                   KK = KK + J
  253:    40         CONTINUE
  254:           END IF
  255:       ELSE
  256: *
  257: *        Form  A  when lower triangle is stored in AP.
  258: *
  259:           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
  260:               DO 60 J = 1,N
  261:                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
  262:                       TEMP1 = ALPHA*Y(J)
  263:                       TEMP2 = ALPHA*X(J)
  264:                       K = KK
  265:                       DO 50 I = J,N
  266:                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
  267:                           K = K + 1
  268:    50                 CONTINUE
  269:                   END IF
  270:                   KK = KK + N - J + 1
  271:    60         CONTINUE
  272:           ELSE
  273:               DO 80 J = 1,N
  274:                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
  275:                       TEMP1 = ALPHA*Y(JY)
  276:                       TEMP2 = ALPHA*X(JX)
  277:                       IX = JX
  278:                       IY = JY
  279:                       DO 70 K = KK,KK + N - J
  280:                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
  281:                           IX = IX + INCX
  282:                           IY = IY + INCY
  283:    70                 CONTINUE
  284:                   END IF
  285:                   JX = JX + INCX
  286:                   JY = JY + INCY
  287:                   KK = KK + N - J + 1
  288:    80         CONTINUE
  289:           END IF
  290:       END IF
  291: *
  292:       RETURN
  293: *
  294: *     End of DSPR2 .
  295: *
  296:       END