SUBROUTINE saxpy (n, sa, sx, incx, sy, incy)  
!     constant times a vector plus a vector.
!     uses unrolled loop for increments equal to one.
!     jack dongarra, linpack, 3/11/78.
!     modified 12/3/93, array(1) declarations changed to array(*)
IMPLICIT NONE

INTEGER, INTENT(IN) :: n, incx, incy
REAL   , INTENT(IN) :: sx(*), sa
REAL   , INTENT(OUT):: sy(*)
integer :: i, ix, iy, m, mp1

if (n.le.0) return  
if (sa.eq.0.0) return  
if (incx.eq.1.and.incy.eq.1) goto 20  
!        code for unequal increments or equal increments
!          not equal to 1
ix = 1  
iy = 1  
if (incx.lt.0) ix = ( - n + 1) * incx + 1  
if (incy.lt.0) iy = ( - n + 1) * incy + 1  
do 10 i = 1, n  
   sy (iy) = sy (iy) + sa * sx (ix)  
   ix = ix + incx  
   iy = iy + incy  
10 end do  
return  
!        code for both increments equal to 1
!        clean-up loop
20 m = mod (n, 4)  
if (m.eq.0) goto 40  
do 30 i = 1, m  
   sy (i) = sy (i) + sa * sx (i)  
30 end do  
if (n.lt.4) return  
40 mp1 = m + 1  
do 50 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)  
50 end do  

return  
END SUBROUTINE saxpy
