dcost Subroutine

      subroutine dcost(n, x, wsave)
         use fftpack_kind, only: rk
         implicit none
         integer, intent(in) :: n
         real(rk), intent(inout) :: wsave(*)
         real(rk), intent(inout) :: x(*)
         real(rk) :: c1, t1, t2, tx2, x1h, x1p3, &
                     xi, xim2
         integer :: i, k, kc, modn, nm1, np1, ns2
         nm1 = n - 1
         np1 = n + 1
         ns2 = n/2
         if (n < 2) return
         if (n == 2) then
            x1h = x(1) + x(2)
            x(2) = x(1) - x(2)
            x(1) = x1h
            return
         elseif (n > 3) then
            c1 = x(1) - x(n)
            x(1) = x(1) + x(n)
            do k = 2, ns2
               kc = np1 - k
               t1 = x(k) + x(kc)
               t2 = x(k) - x(kc)
               c1 = c1 + wsave(kc)*t2
               t2 = wsave(k)*t2
               x(k) = t1 - t2
               x(kc) = t1 + t2
            end do
            modn = mod(n, 2)
            if (modn /= 0) x(ns2 + 1) = x(ns2 + 1) + x(ns2 + 1)
            call dfftf(nm1, x, wsave(n + 1))
            xim2 = x(2)
            x(2) = c1
            do i = 4, n, 2
               xi = x(i)
               x(i) = x(i - 2) - x(i - 1)
               x(i - 1) = xim2
               xim2 = xi
            end do
            if (modn /= 0) x(n) = xim2
            return
         end if
         x1p3 = x(1) + x(3)
         tx2 = x(2) + x(2)
         x(2) = x(1) - x(3)
         x(1) = x1p3 + tx2
         x(3) = x1p3 - tx2
      end subroutine dcost