ddot.f (3542B)
1 !> 2 ! @license Apache-2.0 3 ! 4 ! Copyright (c) 2019 The Stdlib Authors. 5 ! 6 ! Licensed under the Apache License, Version 2.0 (the "License"); 7 ! you may not use this file except in compliance with the License. 8 ! You may obtain a copy of the License at 9 ! 10 ! http://www.apache.org/licenses/LICENSE-2.0 11 ! 12 ! Unless required by applicable law or agreed to in writing, software 13 ! distributed under the License is distributed on an "AS IS" BASIS, 14 ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 15 ! See the License for the specific language governing permissions and 16 ! limitations under the License. 17 !< 18 19 !> Computes the dot product of two double-precision floating-point vectors. 20 ! 21 ! ## Notes 22 ! 23 ! * Modified version of reference BLAS level1 routine (version 3.7.0). Updated to "free form" Fortran 95. 24 ! 25 ! ## Authors 26 ! 27 ! * Univ. of Tennessee 28 ! * Univ. of California Berkeley 29 ! * Univ. of Colorado Denver 30 ! * NAG Ltd. 31 ! 32 ! ## History 33 ! 34 ! * Jack Dongarra, linpack, 3/11/78. 35 ! 36 ! - modified 12/3/93, array(1) declarations changed to array(*) 37 ! 38 ! ## License 39 ! 40 ! From <http://netlib.org/blas/faq.html>: 41 ! 42 ! > The reference BLAS is a freely-available software package. It is available from netlib via anonymous ftp and the World Wide Web. Thus, it can be included in commercial software packages (and has been). We only ask that proper credit be given to the authors. 43 ! > 44 ! > Like all software, it is copyrighted. It is not trademarked, but we do ask the following: 45 ! > 46 ! > * If you modify the source for these routines we ask that you change the name of the routine and comment the changes made to the original. 47 ! > 48 ! > * We will gladly answer any questions regarding the software. If a modification is done, however, it is the responsibility of the person who modified the routine to provide support. 49 ! 50 ! @param {integer} N - number of values over which to compute the dot product 51 ! @param {Array<double>} dx - first array 52 ! @param {integer} strideX - `dx` stride length 53 ! @param {Array<double>} dy - second array 54 ! @param {integer} strideY - `dy` stride length 55 ! @returns {double} the dot product of `dx` and `dy` 56 !< 57 double precision function ddot( N, dx, strideX, dy, strideY ) 58 implicit none 59 ! .. 60 ! Scalar arguments: 61 integer :: strideX, strideY, N 62 ! .. 63 ! Array arguments: 64 double precision, intent(in) :: dx(*), dy(*) 65 ! .. 66 ! Local scalars: 67 double precision :: dtemp 68 integer :: mp1, ix, iy, i, m 69 ! .. 70 ! Intrinsic functions: 71 intrinsic mod 72 ! .. 73 dtemp = 0.0d0 74 ddot = 0.0d0 75 ! .. 76 if ( N <= 0 ) then 77 return 78 end if 79 ! .. 80 ! If both strides are equal to `1`, use unrolled loops... 81 if ( strideX == 1 .AND. strideY == 1 ) then 82 m = mod( N, 5 ) 83 ! .. 84 ! If we have a remainder, do a clean-up loop... 85 if ( m /= 0 ) then 86 do i = 1, m 87 dtemp = dtemp + ( dx( i ) * dy( i ) ) 88 end do 89 end if 90 if ( N < M ) then 91 ddot = dtemp 92 return 93 end if 94 mp1 = m + 1 95 do i = mp1, N, 5 96 dtemp = dtemp + & 97 ( dx( i ) * dy( i ) ) + & 98 ( dx( i+1 ) * dy( i+1 ) ) + & 99 ( dx( i+2 ) * dy( i+2 ) ) + & 100 ( dx( i+3 ) * dy( i+3 ) ) + & 101 ( dx( i+4 ) * dy( i+4 ) ) 102 end do 103 else 104 if ( strideX < 0 ) then 105 ix = ((1-N)*strideX) + 1 106 else 107 ix = 1 108 endif 109 if ( strideY < 0 ) then 110 iy = ((1-N)*strideY) + 1 111 else 112 iy = 1 113 endif 114 do i = 1, N 115 dtemp = dtemp + ( dx( ix ) * dy( iy ) ) 116 ix = ix + strideX 117 iy = iy + strideY 118 end do 119 endif 120 ddot = dtemp 121 return 122 end function ddot