reference/dsyr2f.f

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