daxpy Subroutine

subroutine daxpy(N, Da, Dx, Incx, Dy, Incy)

NAME

daxpy(3f) - [M_odepack::matrix] Compute a constant times a vector plus a vector.

SYNOPSIS

subroutine daxpy(N,Da,Dx,Incx,Dy,Incy)
integer,intent(in)          :: N
real(kind=dp),intent(in)    :: Da
real(kind=dp),intent(in)    :: Dx(*)
integer,intent(in)          :: Incx
real(kind=dp),intent(inout) :: Dy(*)
integer,intent(in)          :: Incy

DESCRIPTION

daxpy(3f) computes a constant times a vector plus a vector. It uses unrolled loops for increments equal to one.

Overwrite double precision DY with double precision DADX + DY. For I = 0 to N-1, replace DY(LY+IINCY) with DADX(LX+IINCX) + DY(LY+I*INCY),

where LX = 1 if INCX .GE. 0, else LX = 1+(1-N)*INCX, and LY is defined in a similar way using INCY.

INPUT OPTIONS

N

number of elements in input vector(s)

DA

double precision scalar multiplier

DX

double precision vector with N elements

INCX

storage spacing between elements of DX

DY

double precision vector with N elements

INCY

storage spacing between elements of DY

RETURNS

DY

double precision result (unchanged if N .LE. 0)

REFERENCES

B L A S Subprogram

C. L. Lawson, R. J. Hanson, D. R. Kincaid and F. T. Krogh, Basic linear algebra subprograms for Fortran usage, Algorithm No. 539, Transactions on Mathematical Software 5, 3 (September 1979), pp. 308-323.

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: N
real(kind=dp), intent(in) :: Da
real(kind=dp), intent(in) :: Dx(*)
integer, intent(in) :: Incx
real(kind=dp), intent(inout) :: Dy(*)
integer, intent(in) :: Incy

Variables

Type Visibility Attributes Name Initial
integer, public :: i
integer, public :: ix
integer, public :: iy
integer, public :: m
integer, public :: mp1
integer, public :: ns

Source Code

subroutine daxpy(N,Da,Dx,Incx,Dy,Incy)
!
integer , intent(in)          :: N
real(kind=dp) , intent(in)    :: Da
real(kind=dp) , intent(in)    :: Dx(*)
integer , intent(in)          :: Incx
real(kind=dp) , intent(inout) :: Dy(*)
integer , intent(in)          :: Incy
!
integer :: i , ix , iy , m , mp1 , ns
!
if ( N<=0 .or. Da==0.0D0 ) return
if ( Incx==Incy ) then
   if ( Incx<1 ) then
   elseif ( Incx==1 ) then
!
!      Code for both increments equal to 1.
!
!      Clean-up loop so remaining vector length is a multiple of 4.
!
      m = mod(N,4)
      if ( m/=0 ) then
         do i = 1 , m
            Dy(i) = Dy(i) + Da*Dx(i)
         enddo
         if ( N<4 ) return
      endif
      mp1 = m + 1
      do i = mp1 , N , 4
         Dy(i) = Dy(i) + Da*Dx(i)
         Dy(i+1) = Dy(i+1) + Da*Dx(i+1)
         Dy(i+2) = Dy(i+2) + Da*Dx(i+2)
         Dy(i+3) = Dy(i+3) + Da*Dx(i+3)
      enddo
      return
   else
!
!      Code for equal, positive, non-unit increments.
!
      ns = N*Incx
      do i = 1 , ns , Incx
         Dy(i) = Da*Dx(i) + Dy(i)
      enddo
      return
   endif
endif
!
!      Code for unequal or nonpositive increments.
!
ix = 1
iy = 1
if ( Incx<0 ) ix = (-N+1)*Incx + 1
if ( Incy<0 ) iy = (-N+1)*Incy + 1
do i = 1 , N
   Dy(iy) = Dy(iy) + Da*Dx(ix)
   ix = ix + Incx
   iy = iy + Incy
enddo
end subroutine daxpy