ddot forms the dot product of two vectors. uses unrolled loops for increments equal to one.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer(kind=ip) | :: | n | ||||
real(kind=wp) | :: | dx(*) | ||||
integer(kind=ip) | :: | incx | ||||
real(kind=wp) | :: | dy(*) | ||||
integer(kind=ip) | :: | incy |
real(wp) function ddot(n, dx, incx, dy, incy) !! ddot forms the dot product of two vectors. !! uses unrolled loops for increments equal to one. integer(ip) :: incx, incy, n real(wp) :: dx(*), dy(*) real(wp) :: dtemp integer(ip) :: i, ix, iy, m, mp1 ddot = 0.0_wp dtemp = 0.0_wp if (n <= 0_ip) return if (incx == 1_ip .and. incy == 1_ip) then ! code for both increments equal to 1 ! clean-up loop m = mod(n, 5_ip) if (m /= 0_ip) then do i = 1_ip, m dtemp = dtemp + dx(i)*dy(i) end do if (n < 5_ip) then ddot = dtemp return end if end if mp1 = m + 1_ip do i = mp1, n, 5_ip dtemp = dtemp + dx(i)*dy(i) + & dx(i + 1_ip)*dy(i + 1_ip) + dx(i + 2_ip)*dy(i + 2_ip) + & dx(i + 3_ip)*dy(i + 3_ip) + dx(i + 4_ip)*dy(i + 4_ip) end do else ! code for unequal increments or equal increments ! not equal to 1 ix = 1_ip iy = 1_ip if (incx < 0_ip) ix = (-n + 1_ip)*incx + 1_ip if (incy < 0_ip) iy = (-n + 1_ip)*incy + 1_ip do i = 1_ip, n dtemp = dtemp + dx(ix)*dy(iy) ix = ix + incx iy = iy + incy end do end if ddot = dtemp end function ddot