share/lapack/blas/fortran/dspr.f

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