subroutine daxpy0(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
! simple version of daxpy written by Greg
c
      double precision dx(n),dy(n),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return

      do 10 i = 1,n
        dy(i) = dy(i) + da*dx(i)
   10 continue
      return
      end

      subroutine daxpy2(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
! simple version of daxpy written by Greg, unrolled by 2
c
      double precision dx(n),dy(n),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return

      do 10 i = 1,n/2,2
        dy(i) = dy(i) + da*dx(i)
        dy(i+1) = dy(i+1) + da*dx(i+1)
   10 continue
      return
      end

      subroutine daxpy4(n,da,dx,incx,dy,incy)
c
c     constant times a vector plus a vector.
! simple version of daxpy written by Greg, unrolled by 4
c
      double precision dx(n),dy(n),da
      integer i,incx,incy,ix,iy,m,mp1,n
c
      if(n.le.0)return
      if (da .eq. 0.0d0) return

      do 10 i = 1,n/4,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)
   10 continue
      return
      end