BLAS/SRC/dspr2.f

   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