src/BLAS/dspr2.f

   1       SUBROUTINE DSPR2(UPLO,N,ALPHA,X,INCX,Y,INCY,AP)
   2 *     .. Scalar Arguments ..
   3       DOUBLE PRECISION ALPHA
   4       INTEGER INCX,INCY,N
   5       CHARACTER UPLO
   6 *     ..
   7 *     .. Array Arguments ..
   8       DOUBLE PRECISION AP(*),X(*),Y(*)
   9 *     ..
  10 *
  11 *  Purpose
  12 *  =======
  13 *
  14 *  DSPR2  performs the symmetric rank 2 operation
  15 *
  16 *     A := alpha*x*y' + alpha*y*x' + A,
  17 *
  18 *  where alpha is a scalar, x and y are n element vectors and A is an
  19 *  n by n symmetric matrix, supplied in packed form.
  20 *
  21 *  Arguments
  22 *  ==========
  23 *
  24 *  UPLO   - CHARACTER*1.
  25 *           On entry, UPLO specifies whether the upper or lower
  26 *           triangular part of the matrix A is supplied in the packed
  27 *           array AP as follows:
  28 *
  29 *              UPLO = 'U' or 'u'   The upper triangular part of A is
  30 *                                  supplied in AP.
  31 *
  32 *              UPLO = 'L' or 'l'   The lower triangular part of A is
  33 *                                  supplied in AP.
  34 *
  35 *           Unchanged on exit.
  36 *
  37 *  N      - INTEGER.
  38 *           On entry, N specifies the order of the matrix A.
  39 *           N must be at least zero.
  40 *           Unchanged on exit.
  41 *
  42 *  ALPHA  - DOUBLE PRECISION.
  43 *           On entry, ALPHA specifies the scalar alpha.
  44 *           Unchanged on exit.
  45 *
  46 *  X      - DOUBLE PRECISION array of dimension at least
  47 *           ( 1 + ( n - 1 )*abs( INCX ) ).
  48 *           Before entry, the incremented array X must contain the n
  49 *           element vector x.
  50 *           Unchanged on exit.
  51 *
  52 *  INCX   - INTEGER.
  53 *           On entry, INCX specifies the increment for the elements of
  54 *           X. INCX must not be zero.
  55 *           Unchanged on exit.
  56 *
  57 *  Y      - DOUBLE PRECISION array of dimension at least
  58 *           ( 1 + ( n - 1 )*abs( INCY ) ).
  59 *           Before entry, the incremented array Y must contain the n
  60 *           element vector y.
  61 *           Unchanged on exit.
  62 *
  63 *  INCY   - INTEGER.
  64 *           On entry, INCY specifies the increment for the elements of
  65 *           Y. INCY must not be zero.
  66 *           Unchanged on exit.
  67 *
  68 *  AP     - DOUBLE PRECISION array of 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. On exit, the array
  75 *           AP is overwritten by the upper triangular part of the
  76 *           updated matrix.
  77 *           Before entry with UPLO = 'L' or 'l', the array AP must
  78 *           contain the lower triangular part of the symmetric matrix
  79 *           packed sequentially, column by column, so that AP( 1 )
  80 *           contains a( 1, 1 ), AP( 2 ) and AP( 3 ) contain a( 2, 1 )
  81 *           and a( 3, 1 ) respectively, and so on. On exit, the array
  82 *           AP is overwritten by the lower triangular part of the
  83 *           updated matrix.
  84 *
  85 *
  86 *  Level 2 Blas routine.
  87 *
  88 *  -- Written on 22-October-1986.
  89 *     Jack Dongarra, Argonne National Lab.
  90 *     Jeremy Du Croz, Nag Central Office.
  91 *     Sven Hammarling, Nag Central Office.
  92 *     Richard Hanson, Sandia National Labs.
  93 *
  94 *
  95 *     .. Parameters ..
  96       DOUBLE PRECISION ZERO
  97       PARAMETER (ZERO=0.0D+0)
  98 *     ..
  99 *     .. Local Scalars ..
 100       DOUBLE PRECISION TEMP1,TEMP2
 101       INTEGER I,INFO,IX,IY,J,JX,JY,K,KK,KX,KY
 102 *     ..
 103 *     .. External Functions ..
 104       LOGICAL LSAME
 105       EXTERNAL LSAME
 106 *     ..
 107 *     .. External Subroutines ..
 108       EXTERNAL XERBLA
 109 *     ..
 110 *
 111 *     Test the input parameters.
 112 *
 113       INFO = 0
 114       IF (.NOT.LSAME(UPLO,'U') .AND. .NOT.LSAME(UPLO,'L')) THEN
 115           INFO = 1
 116       ELSE IF (N.LT.0) THEN
 117           INFO = 2
 118       ELSE IF (INCX.EQ.0) THEN
 119           INFO = 5
 120       ELSE IF (INCY.EQ.0) THEN
 121           INFO = 7
 122       END IF
 123       IF (INFO.NE.0) THEN
 124           CALL XERBLA('DSPR2 ',INFO)
 125           RETURN
 126       END IF
 127 *
 128 *     Quick return if possible.
 129 *
 130       IF ((N.EQ.0) .OR. (ALPHA.EQ.ZERO)) RETURN
 131 *
 132 *     Set up the start points in X and Y if the increments are not both
 133 *     unity.
 134 *
 135       IF ((INCX.NE.1) .OR. (INCY.NE.1)) THEN
 136           IF (INCX.GT.0) THEN
 137               KX = 1
 138           ELSE
 139               KX = 1 - (N-1)*INCX
 140           END IF
 141           IF (INCY.GT.0) THEN
 142               KY = 1
 143           ELSE
 144               KY = 1 - (N-1)*INCY
 145           END IF
 146           JX = KX
 147           JY = KY
 148       END IF
 149 *
 150 *     Start the operations. In this version the elements of the array AP
 151 *     are accessed sequentially with one pass through AP.
 152 *
 153       KK = 1
 154       IF (LSAME(UPLO,'U')) THEN
 155 *
 156 *        Form  A  when upper triangle is stored in AP.
 157 *
 158           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
 159               DO 20 J = 1,N
 160                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
 161                       TEMP1 = ALPHA*Y(J)
 162                       TEMP2 = ALPHA*X(J)
 163                       K = KK
 164                       DO 10 I = 1,J
 165                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
 166                           K = K + 1
 167    10                 CONTINUE
 168                   END IF
 169                   KK = KK + J
 170    20         CONTINUE
 171           ELSE
 172               DO 40 J = 1,N
 173                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
 174                       TEMP1 = ALPHA*Y(JY)
 175                       TEMP2 = ALPHA*X(JX)
 176                       IX = KX
 177                       IY = KY
 178                       DO 30 K = KK,KK + J - 1
 179                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
 180                           IX = IX + INCX
 181                           IY = IY + INCY
 182    30                 CONTINUE
 183                   END IF
 184                   JX = JX + INCX
 185                   JY = JY + INCY
 186                   KK = KK + J
 187    40         CONTINUE
 188           END IF
 189       ELSE
 190 *
 191 *        Form  A  when lower triangle is stored in AP.
 192 *
 193           IF ((INCX.EQ.1) .AND. (INCY.EQ.1)) THEN
 194               DO 60 J = 1,N
 195                   IF ((X(J).NE.ZERO) .OR. (Y(J).NE.ZERO)) THEN
 196                       TEMP1 = ALPHA*Y(J)
 197                       TEMP2 = ALPHA*X(J)
 198                       K = KK
 199                       DO 50 I = J,N
 200                           AP(K) = AP(K) + X(I)*TEMP1 + Y(I)*TEMP2
 201                           K = K + 1
 202    50                 CONTINUE
 203                   END IF
 204                   KK = KK + N - J + 1
 205    60         CONTINUE
 206           ELSE
 207               DO 80 J = 1,N
 208                   IF ((X(JX).NE.ZERO) .OR. (Y(JY).NE.ZERO)) THEN
 209                       TEMP1 = ALPHA*Y(JY)
 210                       TEMP2 = ALPHA*X(JX)
 211                       IX = JX
 212                       IY = JY
 213                       DO 70 K = KK,KK + N - J
 214                           AP(K) = AP(K) + X(IX)*TEMP1 + Y(IY)*TEMP2
 215                           IX = IX + INCX
 216                           IY = IY + INCY
 217    70                 CONTINUE
 218                   END IF
 219                   JX = JX + INCX
 220                   JY = JY + INCY
 221                   KK = KK + N - J + 1
 222    80         CONTINUE
 223           END IF
 224       END IF
 225 *
 226       RETURN
 227 *
 228 *     End of DSPR2 .
 229 *
 230       END