C-----------------------------------------------------------------------
      subroutine resisub(u,rhs,residue,M,tcall)
C-----------------------------------------------------------------------
c
c the residual for the linear problem
c
********** INPUTS 
****       u: Approximate solution 
****       rhs: Right hand side of differential equation
****!!!!!  M: A power of 2 - 
****       tcall: Let tcall = 1 the first time resisub is called
****              Let tcall = 2 thereafter
*********  OUTPUTS
****       residue: Lu - f , where L is the operator we spoke about 
****       u,rhs,tcall,M: All remain unchanged

*********  REMARKS
****       M_MAX should be changed in the file  limits.h so that
****       the dimension of  u  and rhs does not exceed M_MAX+1   
****       Anytime M_MAX is changed, the file residue.f should be
****       recompiled

****       The file  limits.h should be included in the directory
****       one is working in

****       The active elements of u and rhs are 1 through M+1 (not M)
****       Example:  Suppose N elements of u and rhs N are active in the
****       main program - The call to resisub should be 
****       resisub(r,rhs,residue,N-1,tcall)
****       Be sure N-1 is a power of 2

****       f is not just the right hand transformed to Chebyshev space
****       First, f must be conditioned by the preconditioner (B2) 
****       Then, two boundary conditions are inserted in the 
****       first two empty slots


      implicit none
      integer M,M_MAX,i,j,tcall
      include 'limits.h'
      include 'fcos2.h'
      double precision residue(M_MAX+1),sum
      double precision B2(M_MAX+1,M_MAX+1)
      double precision pre(M_MAX+1,M_MAX+1)
      double precision x(0:M_MAX), fftarry(M_MAX+2),pi
      double precision pim, rhs(M_MAX+1)
      double precision u(M_MAX+1),Q,sumbc1,sumbc2
      save

      if (tcall .le. 1) then
      pi = acos(-1.d0)
      pim = pi/dble(M)
      call fftini(M)
      
      do i = 0,M
         x(i) = cos(pim*dble(i))
      end do
      call integ(B2,pre,M,M_MAX)
      do i = 1,M+1
         B2(1,i) = 1.d0
         B2(2,i) = (-1.d0)**(i-1)
      end do

      end if


      do i = 1,M+1
         residue(i) = 0.d0
      end do   
   
      do i = 1,M+1
         fftarry(i) = u(i)
      end do   
      call cosft(M,fftarry,um,vm,wm,ibitm,-1)   
      do i = 1,M+1
         fftarry(i) = fftarry(i)*Q(x(i-1))
      end do   
      call cosft(M,fftarry,um,vm,wm,ibitm,1)     


      do i = 3,M+1
         sum = 0.d0
         do j = i-2,min(i+2,M+1),2
            sum = sum + B2(i,j)*fftarry(j)
         end do
         residue(i) = sum
      end do

      sumbc1 = 0.d0
      sumbc2 = 0.d0
      do i = 1,M+1
         sumbc1 = sumbc1 + B2(1,i)*u(i)
         sumbc2 = sumbc2 + B2(2,i)*u(i)
      end do   

      residue(1) = sumbc1
      residue(2) = sumbc2

      do i = 3,M+1
         residue(i) = residue(i)+u(i)
      end do

      do i = 1,M+1
         residue(i) = residue(i) - rhs(i)
      end do

      end

      double precision function Q(x)
      double precision Q,x
      Q = sin(x)
      end

      subroutine integ(B2,pre,M,M_MAX)
      implicit none
      integer M,M_MAX,i,j
      double precision B2(M_MAX+1,M_MAX+1)
      double precision pre(M_MAX+1,M_MAX+1)
      do i = 3,M+1
         pre(i,i) = 1.d0/(dble(2*i-4)*dble(2*i-2))
         pre(i,i+2) = -.5d0/(dble(i)*dble(i-2))
         pre(i,i+4) = 1.d0/(dble(2*i)*dble(2*i-2))
      end do   
      pre(3,3) = 1.d0/4.d0

      do i = 3,M+1
         do j = 1,M-1
            B2(i,j) = pre(i,j+2)
         end do
      end do   

      end 

c *************************************************

      subroutine cosft(n,y,u,v,w,ibit,isign)

c....   this subroutine calculates the cosine transform of
c....   an array of n+1 real numbers, f_j (j = 0,1,...,n) :
c....
c....                 n            pi
c....        f_j = sum    a_k*cos( -- *j*k)
c....                 k=0          n
c....
c....
c....   variables: (i=input, o=output)
c....   ==========
c....   n     integer          (i)   n+1=number of real numbers
c....   y     real             (i,o) dimension: nmax+2
c....                                array containing the n+1 real numbers
c....   u     real array       (i)   dimension: nmax-2
c....                                trigonometric factors used in cosft
c....   v     real array       (i)   dimension: (nmax/2)+2
c....                                trigonometric factors used in realft
c....   w     real array       (i)   dimension: nmax/2
c....                                trigonometric factors used in cfft
c....   ibit  integer array    (i)   dimension: nmax
c....                                bit-reversal array used in cfft
c....   isign integer          (i)   specifies direction of transformation (see
c....                                below)
c....
c....   isign=1:  forward transformation
c....   ========
c....      input:  y(i)=f_j    ... n+1 real function values
c....      output: y(i)=a_k    ... n+1 cosine coefficients,
c....
c....   isign=-1: invers transformation
c....   =========
c....      input:  y(i)=a_k    ... n+1 cosine coefficients
c....      output: y(i)=f_j    ... n+1 real function values
c....
c....               -------------------------------------------
c....
c....   note : a_0 and a_(n) do not contain any weight factor (see
c....            the above def. of f_j)
c....
c.........................................................................

      implicit none

      integer n, ibit(*), m, j, isign
      double precision wr,wi,sum,y1,y2,c1
      double precision y(*), u(*), v(*), w(*)

      if (isign.eq.-1) then
        y(1)=y(1)*2.d0
        y(n+1)=y(n+1)*2.d0
      end if
      sum =0.5d0*(y(1)-y(n+1))
      y(1)=0.5d0*(y(1)+y(n+1))
      m=n/2
      do 11 j=1,m-1
        wr=u(2*(j-1)+1)
        wi=u(2*(j-1)+2)
        y1=0.5d0*(y(j+1)+y(n-j+1))
        y2=(y(j+1)-y(n-j+1))
        y(j+1)  =y1-wi*y2
        y(n-j+1)=y1+wi*y2
        sum=sum+wr*y2
11    continue
      call realft(n,y,v,w,ibit,1)
      sum=sum*2.d0/dble(n)
      y(2)=sum
      do 12 j=4,n,2
        sum=sum+y(j)
        y(j)=sum
12    continue
      if (isign.eq.-1) then
        y(1)=dble(n)*y(1)
        y(n+1)=dble(n)*y(n+1)
        c1=n/2.d0
        do j=2,n
          y(j)=c1*y(j)
        end do
      end if
      return
      end

c ******************************************************

      subroutine realft(n,data,v,w,ibit,isign)

c....   this subroutine calculates the fourier transform of
c....   an array of n real numbers, f_j (j = 0,1,...,n-1) :
c....
c....        f_j = a_0 + a_(n/2)*cos(pi*j) +
c....
c....                 (n/2)-1           2*pi                   2*pi
c....              sum       { a_k*cos( ----- *j*k) + b_k*sin( ---- *j*k) }
c....                  k=1                n                      n
c....
c....
c....   variables: (i=input, o=output)
c....   ==========
c....   n     integer          (i)   number of real numbers
c....   data  real             (i,o) dimension: nmax+2
c....                                array containing the n real numbers
c....   v     real array       (i)   dimension: (nmax/2)+2
c....                                trigonometric factors used in realft
c....   w     real array       (i)   dimension: nmax/2
c....                                trigonometric factors used in cfft
c....   ibit  integer array    (i)   dimension: nmax
c....                                bit-reversal array used in cfft
c....   isign integer          (i)   specifies direction of transformation (see
c....                                below)
c....
c....   isign=1:  forward transformation
c....   ========
c....      input:  data(i)=f_j           ... n real function values
c....      output: data(i)=a_k and b_k   ... n real fourier coefficients,
c....                                        organized as:
c....                                         data(1)    = a_0
c....                                         data(2)    = 0
c....                                         data(2k+1) = a_k, k=1 ... (n/2)-1
c....                                         data(2k+2) = b_k, k=1 ... (n/2)-1
c....                                         data(n+1)  = a_(n/2)
c....                                         data(n+2)  = 0
c....
c....   isign=-1: invers transformation
c....   =========
c....      input:  data(i)=a_k and b_k    ... n real fourier coefficients
c....      output: data(i)=f_j            ... n real function values
c....
c....               -------------------------------------------
c....
c....   note : - definition of "forwards" and "inverse"
c....          - the full value of n is passed in the call to realft
c....          - the array data must be dimensioned to nmax+2
c....          - cosine and sine formulation rather than exponential (which
c....              means a change in sign of the sin terms)
c....          - a_0 and a_(n/2) do not contain any weight factor (see
c....              the above def. of f_j)
c....
c.........................................................................

      implicit none

      integer n, isign, ibit(*), np3, i, i1, i2, i3, i4
      double precision data(*), v(*), w(*)
      double precision c1, c2, h1r, h1i, h2r, h2i, wr, wi

      c1=0.5d0
      if (isign.eq.1) then
        c2=-0.5d0
        call cfft(n/2,data,w,ibit,+1)
        data(n+1)=data(1)
        data(n+2)=data(2)
      else
        c2=0.5d0
        data(1)  =data(1)*2.d0
        data(n+1)=data(n+1)*2.d0
      endif
      np3=n+3
      do 11 i=1,n/4+1
        i1=2*i-1
        i2=i1+1
        i3=np3-i2
        i4=i3+1
        h1r= c1*(data(i1)+data(i3))
        h1i= c1*(data(i2)-data(i4))*isign
        h2r=-c2*(data(i2)+data(i4))*isign
        h2i= c2*(data(i1)-data(i3))
        wr=v(2*(i-1)+1)
        wi=v(2*(i-1)+2)*isign
        data(i1)= h1r+wr*h2r+wi*h2i
        data(i2)=isign*(-h1i-wr*h2i+wi*h2r)
        data(i3)= h1r-wr*h2r-wi*h2i
        data(i4)=isign*(h1i-wr*h2i+wi*h2r)
11    continue
      if (isign.eq.1) then
        data(2)=  0.d0
        data(n+2)=0.d0
        data(1)  =data(1)/2.d0
        data(n+1)=data(n+1)/2.d0
      else
        call cfft(n/2,data,w,ibit,-1)
        data(n+1)=0.d0
        data(n+2)=0.d0
      endif
      return
      end

 
c ******************************************************

	subroutine cfft(n,d,w,ibit,isign)

c....   this subroutine calculates the complex fourier transform of
c....   an array of n complex numbers, f_j (j = 0,1,...,n-1) :
c....
c....                 n-1               2*pi
c....        f_j = sum     c_k * exp(i* ----- * j*k)
c....                 k=0                 n
c....
c....
c....   variables: (i=input, o=output)
c....   ==========
c....   n     integer          (i)   number of complex numbers
c....   d     complex or real  (i,o) dimension: nmax(=number of real numbers)
c....                                array containing n complex numbers,
c....                                or equivalently 2n real numbers
c....                                designating alternately the
c....                                real and imaginary parts
c....   w     real array       (i)   dimension: nmax/2
c....                                trigonometric factors
c....   ibit  integer array    (i)   dimension: nmax
c....                                bit-reversal array
c....   isign integer          (i)   specifies direction of transformation (see
c....                                below)
c....
c....   isign=1:  forward transformation
c....   ========
c....      input:  d(i)=f_j   ... n complex function values
c....      output: d(i)=c_k   ... n complex fourier coefficients
c....
c....   isign=-1: invers transformation
c....   =========
c....      input:  d(i)=c_k   ... n complex fourier coefficients
c....      output: d(i)=f_j   ... n complex function values
c....
c....               -------------------------------------------
c....
c....   note !!!: the definition of forward and invers transformations
c....               is opposite of the definition in "numerical recipies" !!
c....             the values of c_k have been normalized by division with n,
c....               as they should.
c....
c.........................................................................

        implicit none

	integer n,m,i,j,mmax,k,kd,istep,isign,ibit(*),nn
	double precision d(*),tempr,tempi,tr,ti,w(*)

        nn=2*n
	do 10 i=1,nn,2
                j=ibit(i)
		if(j.gt.i) then
			tempr = d(j)
			tempi = d(j+1)
			d(j)  = d(i)
			d(j+1)= d(i+1)
			d(i)  = tempr
			d(i+1)= tempi
		endif
   10	continue
	mmax = 2
	kd= n/2
  11	if(nn.gt.mmax) then
		istep = 2*mmax
		do 13 m=1,mmax,2
			k = 1+(m-1)*kd
			do 12 i=m,nn,istep
                                j=i+mmax
				tr=d(j)*w(k)+d(j+1)*w(k+1)*isign
				ti=-d(j)*w(k+1)*isign+d(j+1)*w(k)
				d(j)=d(i)-tr
				d(j+1)=d(i+1)-ti
				d(i)=d(i)+tr
				d(i+1)=d(i+1)+ti
   12			continue
   13		continue
		kd = kd/2
		mmax = istep
		goto 11
	endif

        if (isign.eq.1) then
          do 20 i=1,nn
            d(i) = d(i)/dble(n)
20        continue
        end if

	return
	end

 
c *************************************************

      subroutine fftini(m)

c....   this subroutine initializes the bit-reversal
c....   and trigonometric arrays used in
c....   cosft, realft, and cfft (and therefore also in cosfou).
c....   the index m refers to cosft, and n refers to realft, which
c....   are the two different transformations in cosfou.

      implicit none

      integer m

      integer i, j, mm
      double precision pi, factor
      parameter(pi=3.1415926535897932384626d0)

      include 'limits.h'
      include 'fcos2.h'

      j=1
      do i=1,m,2
        ibitm(i)=j
        mm=m/2
 1      if (mm.ge.2 .and. j.gt.mm) then
          j=j-mm
          mm=mm/2
          go to 1
        end if
        j=j+mm
      end do


      factor = pi/dble(m)
      do i=1, m/2-1
        um(2*(i-1)+1) = cos(factor*i)
        um(2*(i-1)+2) = sin(factor*i)
      end do

      factor = 2*pi/dble(m)
      do i=1, m/4+1
        vm(2*(i-1)+1) = cos(factor*(i-1))
        vm(2*(i-1)+2) = sin(factor*(i-1))
      end do

      factor = 2*pi/(m/2)
      do i=1, m/4
        wm(2*(i-1)+1) = cos(factor*(i-1))
        wm(2*(i-1)+2) = sin(factor*(i-1))
      end do

      end





