DAXPY constant times a vector plus a vector. uses unrolled loops for increments equal to one.
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
integer(kind=ip) | :: | n | ||||
real(kind=wp) | :: | da | ||||
real(kind=wp) | :: | dx(*) | ||||
integer(kind=ip) | :: | incx | ||||
real(kind=wp) | :: | dy(*) | ||||
integer(kind=ip) | :: | incy |
subroutine daxpy(n, da, dx, incx, dy, incy) !! DAXPY constant times a vector plus a vector. !! uses unrolled loops for increments equal to one. real(wp) :: da integer(ip) :: incx, incy, n real(wp) :: dx(*), dy(*) integer(ip) :: i, ix, iy, m, mp1 if (n <= 0_ip) return if (da == 0.0_wp) return if (incx == 1_ip .and. incy == 1_ip) then ! code for both increments equal to 1 ! clean-up loop m = mod(n, 4_ip) if (m /= 0_ip) then do i = 1_ip, m dy(i) = dy(i) + da*dx(i) end do end if if (n < 4_ip) return mp1 = m + 1_ip do i = mp1, n, 4_ip dy(i) = dy(i) + da*dx(i) dy(i + 1_ip) = dy(i + 1_ip) + da*dx(i + 1_ip) dy(i + 2_ip) = dy(i + 2_ip) + da*dx(i + 2_ip) dy(i + 3_ip) = dy(i + 3_ip) + da*dx(i + 3_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 dy(iy) = dy(iy) + da*dx(ix) ix = ix + incx iy = iy + incy end do end if end subroutine daxpy