share/odepack/fortran/daxpy.f

   1 *DECK DAXPY
   2       SUBROUTINE DAXPY (N, DA, DX, INCX, DY, INCY)
   3 C***BEGIN PROLOGUE  DAXPY
   4 C***PURPOSE  Compute a constant times a vector plus a vector.
   5 C***CATEGORY  D1A7
   6 C***TYPE      DOUBLE PRECISION (SAXPY-S, DAXPY-D, CAXPY-C)
   7 C***KEYWORDS  BLAS, LINEAR ALGEBRA, TRIAD, VECTOR
   8 C***AUTHOR  Lawson, C. L., (JPL)
   9 C           Hanson, R. J., (SNLA)
  10 C           Kincaid, D. R., (U. of Texas)
  11 C           Krogh, F. T., (JPL)
  12 C***DESCRIPTION
  13 C
  14 C                B L A S  Subprogram
  15 C    Description of Parameters
  16 C
  17 C     --Input--
  18 C        N  number of elements in input vector(s)
  19 C       DA  double precision scalar multiplier
  20 C       DX  double precision vector with N elements
  21 C     INCX  storage spacing between elements of DX
  22 C       DY  double precision vector with N elements
  23 C     INCY  storage spacing between elements of DY
  24 C
  25 C     --Output--
  26 C       DY  double precision result (unchanged if N .LE. 0)
  27 C
  28 C     Overwrite double precision DY with double precision DA*DX + DY.
  29 C     For I = 0 to N-1, replace  DY(LY+I*INCY) with DA*DX(LX+I*INCX) +
  30 C       DY(LY+I*INCY),
  31 C     where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
  32 C     defined in a similar way using INCY.
  33 C
  34 C***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
  35 C                 Krogh, Basic linear algebra subprograms for Fortran
  36 C                 usage, Algorithm No. 539, Transactions on Mathematical
  37 C                 Software 5, 3 (September 1979), pp. 308-323.
  38 C***ROUTINES CALLED  (NONE)
  39 C***REVISION HISTORY  (YYMMDD)
  40 C   791001  DATE WRITTEN
  41 C   890831  Modified array declarations.  (WRB)
  42 C   890831  REVISION DATE from Version 3.2
  43 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  44 C   920310  Corrected definition of LX in DESCRIPTION.  (WRB)
  45 C   920501  Reformatted the REFERENCES section.  (WRB)
  46 C***END PROLOGUE  DAXPY
  47       DOUBLE PRECISION DX(*), DY(*), DA
  48 C***FIRST EXECUTABLE STATEMENT  DAXPY
  49       IF (N.LE.0 .OR. DA.EQ.0.0D0) RETURN
  50       IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
  51 C
  52 C     Code for unequal or nonpositive increments.
  53 C
  54     5 IX = 1
  55       IY = 1
  56       IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
  57       IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
  58       DO 10 I = 1,N
  59         DY(IY) = DY(IY) + DA*DX(IX)
  60         IX = IX + INCX
  61         IY = IY + INCY
  62    10 CONTINUE
  63       RETURN
  64 C
  65 C     Code for both increments equal to 1.
  66 C
  67 C     Clean-up loop so remaining vector length is a multiple of 4.
  68 C
  69    20 M = MOD(N,4)
  70       IF (M .EQ. 0) GO TO 40
  71       DO 30 I = 1,M
  72         DY(I) = DY(I) + DA*DX(I)
  73    30 CONTINUE
  74       IF (N .LT. 4) RETURN
  75    40 MP1 = M + 1
  76       DO 50 I = MP1,N,4
  77         DY(I) = DY(I) + DA*DX(I)
  78         DY(I+1) = DY(I+1) + DA*DX(I+1)
  79         DY(I+2) = DY(I+2) + DA*DX(I+2)
  80         DY(I+3) = DY(I+3) + DA*DX(I+3)
  81    50 CONTINUE
  82       RETURN
  83 C
  84 C     Code for equal, positive, non-unit increments.
  85 C
  86    60 NS = N*INCX
  87       DO 70 I = 1,NS,INCX
  88         DY(I) = DA*DX(I) + DY(I)
  89    70 CONTINUE
  90       RETURN
  91       END