share/odepack/fortran/ddot.f

   1 *DECK DDOT
   2       DOUBLE PRECISION FUNCTION DDOT (N, DX, INCX, DY, INCY)
   3 C***BEGIN PROLOGUE  DDOT
   4 C***PURPOSE  Compute the inner product of two vectors.
   5 C***CATEGORY  D1A4
   6 C***TYPE      DOUBLE PRECISION (SDOT-S, DDOT-D, CDOTU-C)
   7 C***KEYWORDS  BLAS, INNER PRODUCT, LINEAR ALGEBRA, 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       DX  double precision vector with N elements
  20 C     INCX  storage spacing between elements of DX
  21 C       DY  double precision vector with N elements
  22 C     INCY  storage spacing between elements of DY
  23 C
  24 C     --Output--
  25 C     DDOT  double precision dot product (zero if N .LE. 0)
  26 C
  27 C     Returns the dot product of double precision DX and DY.
  28 C     DDOT = sum for I = 0 to N-1 of  DX(LX+I*INCX) * DY(LY+I*INCY),
  29 C     where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is
  30 C     defined in a similar way using INCY.
  31 C
  32 C***REFERENCES  C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T.
  33 C                 Krogh, Basic linear algebra subprograms for Fortran
  34 C                 usage, Algorithm No. 539, Transactions on Mathematical
  35 C                 Software 5, 3 (September 1979), pp. 308-323.
  36 C***ROUTINES CALLED  (NONE)
  37 C***REVISION HISTORY  (YYMMDD)
  38 C   791001  DATE WRITTEN
  39 C   890831  Modified array declarations.  (WRB)
  40 C   890831  REVISION DATE from Version 3.2
  41 C   891214  Prologue converted to Version 4.0 format.  (BAB)
  42 C   920310  Corrected definition of LX in DESCRIPTION.  (WRB)
  43 C   920501  Reformatted the REFERENCES section.  (WRB)
  44 C***END PROLOGUE  DDOT
  45       DOUBLE PRECISION DX(*), DY(*)
  46 C***FIRST EXECUTABLE STATEMENT  DDOT
  47       DDOT = 0.0D0
  48       IF (N .LE. 0) RETURN
  49       IF (INCX .EQ. INCY) IF (INCX-1) 5,20,60
  50 C
  51 C     Code for unequal or nonpositive increments.
  52 C
  53     5 IX = 1
  54       IY = 1
  55       IF (INCX .LT. 0) IX = (-N+1)*INCX + 1
  56       IF (INCY .LT. 0) IY = (-N+1)*INCY + 1
  57       DO 10 I = 1,N
  58         DDOT = DDOT + DX(IX)*DY(IY)
  59         IX = IX + INCX
  60         IY = IY + INCY
  61    10 CONTINUE
  62       RETURN
  63 C
  64 C     Code for both increments equal to 1.
  65 C
  66 C     Clean-up loop so remaining vector length is a multiple of 5.
  67 C
  68    20 M = MOD(N,5)
  69       IF (M .EQ. 0) GO TO 40
  70       DO 30 I = 1,M
  71          DDOT = DDOT + DX(I)*DY(I)
  72    30 CONTINUE
  73       IF (N .LT. 5) RETURN
  74    40 MP1 = M + 1
  75       DO 50 I = MP1,N,5
  76       DDOT = DDOT + DX(I)*DY(I) + DX(I+1)*DY(I+1) + DX(I+2)*DY(I+2) +
  77      1              DX(I+3)*DY(I+3) + DX(I+4)*DY(I+4)
  78    50 CONTINUE
  79       RETURN
  80 C
  81 C     Code for equal, positive, non-unit increments.
  82 C
  83    60 NS = N*INCX
  84       DO 70 I = 1,NS,INCX
  85         DDOT = DDOT + DX(I)*DY(I)
  86    70 CONTINUE
  87       RETURN
  88       END