saxpy Subroutine

  subroutine saxpy(n,sa,sx,incx,sy,incy)

    !! Routine to scale a vector by a constant, then add another vector
    !! author: Jack Dongarra, Linpack
    !! author: P J Knight, CCFE, Culham Science Centre
    !! n        : input integer : order of the matrices sx, sy
    !! sa       : input real array : constant multiplier
    !! sx(n*incx) : input real array : matrix to be scaled
    !! incx     : input integer : interval in storage between sx array elements
    !! sy(n*incy) : input/output real array : On entry, matrix being added;
    !! On exit, the final result
    !! incy     : input integer : interval in storage between sy array elements
    !! This routine calculates <CODE>sa*sx(:) + sy(:)</CODE>,
    !! using unrolled loops for increments equal to 1.
    !!     !
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    implicit none

    !  Arguments

    integer, intent(in) :: n,incx,incy
    real(dp), intent(in) :: sa
    real(dp), dimension(n*incx), intent(in) :: sx
    real(dp), dimension(n*incy), intent(inout) :: sy

    !  Local variables

    integer :: i,ix,iy,m,mp1

    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    if ((n <= 0).or.(sa == 0.0D0)) return

    if ((incx /= 1).or.(incy /= 1)) then

       ix = 1 ; iy = 1
       if (incx < 0) ix = (-n+1)*incx + 1
       if (incy < 0) iy = (-n+1)*incy + 1
       do i = 1,n
          sy(iy) = sy(iy) + sa*sx(ix)
          ix = ix + incx
          iy = iy + incy
       end do

    else

       m = mod(n,4)
       if (m /= 0) then
          do i = 1,m
             sy(i) = sy(i) + sa*sx(i)
          end do
          if (n < 4) return
       end if

       mp1 = m + 1
       do i = mp1,n,4
          sy(i)     = sy(i)     + sa*sx(i)
          sy(i + 1) = sy(i + 1) + sa*sx(i + 1)
          sy(i + 2) = sy(i + 2) + sa*sx(i + 2)
          sy(i + 3) = sy(i + 3) + sa*sx(i + 3)
       end do

    end if

  end subroutine saxpy