share/lapack/blas/fortran/dsyr.f

   1       SUBROUTINE DSYR  ( UPLO, N, ALPHA, X, INCX, A, LDA )
   2 *     .. Scalar Arguments ..
   3       DOUBLE PRECISION   ALPHA
   4       INTEGER            INCX, LDA, N
   5       CHARACTER*1        UPLO
   6 *     .. Array Arguments ..
   7       DOUBLE PRECISION   A( LDA, * ), X( * )
   8 *     ..
   9 *
  10 *  Purpose
  11 *  =======
  12 *
  13 *  DSYR   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.
  19 *
  20 *  Parameters
  21 *  ==========
  22 *
  23 *  UPLO   - CHARACTER*1.
  24 *           On entry, UPLO specifies whether the upper or lower
  25 *           triangular part of the array A is to be referenced as
  26 *           follows:
  27 *
  28 *              UPLO = 'U' or 'u'   Only the upper triangular part of A
  29 *                                  is to be referenced.
  30 *
  31 *              UPLO = 'L' or 'l'   Only the lower triangular part of A
  32 *                                  is to be referenced.
  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 *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, n ).
  57 *           Before entry with  UPLO = 'U' or 'u', the leading n by n
  58 *           upper triangular part of the array A must contain the upper
  59 *           triangular part of the symmetric matrix and the strictly
  60 *           lower triangular part of A is not referenced. On exit, the
  61 *           upper triangular part of the array A is overwritten by the
  62 *           upper triangular part of the updated matrix.
  63 *           Before entry with UPLO = 'L' or 'l', the leading n by n
  64 *           lower triangular part of the array A must contain the lower
  65 *           triangular part of the symmetric matrix and the strictly
  66 *           upper triangular part of A is not referenced. On exit, the
  67 *           lower triangular part of the array A is overwritten by the
  68 *           lower triangular part of the updated matrix.
  69 *
  70 *  LDA    - INTEGER.
  71 *           On entry, LDA specifies the first dimension of A as declared
  72 *           in the calling (sub) program. LDA must be at least
  73 *           max( 1, n ).
  74 *           Unchanged on exit.
  75 *
  76 *
  77 *  Level 2 Blas routine.
  78 *
  79 *  -- Written on 22-October-1986.
  80 *     Jack Dongarra, Argonne National Lab.
  81 *     Jeremy Du Croz, Nag Central Office.
  82 *     Sven Hammarling, Nag Central Office.
  83 *     Richard Hanson, Sandia National Labs.
  84 *
  85 *
  86 *     .. Parameters ..
  87       DOUBLE PRECISION   ZERO
  88       PARAMETER        ( ZERO = 0.0D+0 )
  89 *     .. Local Scalars ..
  90       DOUBLE PRECISION   TEMP
  91       INTEGER            I, INFO, IX, J, JX, KX
  92 *     .. External Functions ..
  93       LOGICAL            LSAME
  94       EXTERNAL           LSAME
  95 *     .. External Subroutines ..
  96       EXTERNAL           XERBLA
  97 *     .. Intrinsic Functions ..
  98       INTRINSIC          MAX
  99 *     ..
 100 *     .. Executable Statements ..
 101 *
 102 *     Test the input parameters.
 103 *
 104       INFO = 0
 105       IF     ( .NOT.LSAME( UPLO, 'U' ).AND.
 106      $         .NOT.LSAME( UPLO, 'L' )      )THEN
 107          INFO = 1
 108       ELSE IF( N.LT.0 )THEN
 109          INFO = 2
 110       ELSE IF( INCX.EQ.0 )THEN
 111          INFO = 5
 112       ELSE IF( LDA.LT.MAX( 1, N ) )THEN
 113          INFO = 7
 114       END IF
 115       IF( INFO.NE.0 )THEN
 116          CALL XERBLA( 'DSYR  ', INFO )
 117          RETURN
 118       END IF
 119 *
 120 *     Quick return if possible.
 121 *
 122       IF( ( N.EQ.0 ).OR.( ALPHA.EQ.ZERO ) )
 123      $   RETURN
 124 *
 125 *     Set the start point in X if the increment is not unity.
 126 *
 127       IF( INCX.LE.0 )THEN
 128          KX = 1 - ( N - 1 )*INCX
 129       ELSE IF( INCX.NE.1 )THEN
 130          KX = 1
 131       END IF
 132 *
 133 *     Start the operations. In this version the elements of A are
 134 *     accessed sequentially with one pass through the triangular part
 135 *     of A.
 136 *
 137       IF( LSAME( UPLO, 'U' ) )THEN
 138 *
 139 *        Form  A  when A is stored in upper triangle.
 140 *
 141          IF( INCX.EQ.1 )THEN
 142             DO 20, J = 1, N
 143                IF( X( J ).NE.ZERO )THEN
 144                   TEMP = ALPHA*X( J )
 145                   DO 10, I = 1, J
 146                      A( I, J ) = A( I, J ) + X( I )*TEMP
 147    10             CONTINUE
 148                END IF
 149    20       CONTINUE
 150          ELSE
 151             JX = KX
 152             DO 40, J = 1, N
 153                IF( X( JX ).NE.ZERO )THEN
 154                   TEMP = ALPHA*X( JX )
 155                   IX   = KX
 156                   DO 30, I = 1, J
 157                      A( I, J ) = A( I, J ) + X( IX )*TEMP
 158                      IX        = IX        + INCX
 159    30             CONTINUE
 160                END IF
 161                JX = JX + INCX
 162    40       CONTINUE
 163          END IF
 164       ELSE
 165 *
 166 *        Form  A  when A is stored in lower triangle.
 167 *
 168          IF( INCX.EQ.1 )THEN
 169             DO 60, J = 1, N
 170                IF( X( J ).NE.ZERO )THEN
 171                   TEMP = ALPHA*X( J )
 172                   DO 50, I = J, N
 173                      A( I, J ) = A( I, J ) + X( I )*TEMP
 174    50             CONTINUE
 175                END IF
 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, I = J, N
 184                      A( I, J ) = A( I, J ) + X( IX )*TEMP
 185                      IX        = IX        + INCX
 186    70             CONTINUE
 187                END IF
 188                JX = JX + INCX
 189    80       CONTINUE
 190          END IF
 191       END IF
 192 *
 193       RETURN
 194 *
 195 *     End of DSYR  .
 196 *
 197       END