File:  [local] / rpl / lapack / blas / dsyr2.f
Revision 1.1.1.1 (vendor branch): download - view: text, annotated - select for diffs - revision graph
Tue Jan 26 15:22:45 2010 UTC (14 years, 3 months ago) by bertrand
Branches: JKB
CVS tags: start, rpl-4_0_14, rpl-4_0_13, rpl-4_0_12, rpl-4_0_11, rpl-4_0_10


Commit initial.

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

CVSweb interface <joel.bertrand@systella.fr>