home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / blas / daxpy.f < prev    next >
Text File  |  1996-09-28  |  1KB  |  48 lines
  1.       SUBROUTINE DAXPY(N,DA,DX,INCX,DY,INCY)
  2. C
  3. C     CONSTANT TIMES A VECTOR PLUS A VECTOR.
  4. C     USES UNROLLED LOOPS FOR INCREMENTS EQUAL TO ONE.
  5. C     JACK DONGARRA, LINPACK, 3/11/78.
  6. C
  7.       DOUBLE PRECISION DX(1),DY(1),DA
  8.       INTEGER I,INCX,INCY,IXIY,M,MP1,N
  9. C
  10.       IF(N.LE.0)RETURN
  11.       IF (DA .EQ. 0.0D0) RETURN
  12.       IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
  13. C
  14. C        CODE FOR UNEQUAL INCREMENTS OR EQUAL INCREMENTS
  15. C          NOT EQUAL TO 1
  16. C
  17.       IX = 1
  18.       IY = 1
  19.       IF(INCX.LT.0)IX = (-N+1)*INCX + 1
  20.       IF(INCY.LT.0)IY = (-N+1)*INCY + 1
  21.       DO 10 I = 1,N
  22.         DY(IY) = DY(IY) + DA*DX(IX)
  23.         IX = IX + INCX
  24.         IY = IY + INCY
  25.    10 CONTINUE
  26.       RETURN
  27. C
  28. C        CODE FOR BOTH INCREMENTS EQUAL TO 1
  29. C
  30. C
  31. C        CLEAN-UP LOOP
  32. C
  33.    20 M = MOD(N,4)
  34.       IF( M .EQ. 0 ) GO TO 40
  35.       DO 30 I = 1,M
  36.         DY(I) = DY(I) + DA*DX(I)
  37.    30 CONTINUE
  38.       IF( N .LT. 4 ) RETURN
  39.    40 MP1 = M + 1
  40.       DO 50 I = MP1,N,4
  41.         DY(I) = DY(I) + DA*DX(I)
  42.         DY(I + 1) = DY(I + 1) + DA*DX(I + 1)
  43.         DY(I + 2) = DY(I + 2) + DA*DX(I + 2)
  44.         DY(I + 3) = DY(I + 3) + DA*DX(I + 3)
  45.    50 CONTINUE
  46.       RETURN
  47.       END
  48.