C-----------------------------------------------------------------------
      program broyspec
C-----------------------------------------------------------------------
      implicit none
      integer M, M_MAX ,IRE,i,j,itmax
      include 'limits.h'
      include 'fcos2.h'
      double precision r(M_MAX+1),sum
      double precision pre(M_MAX+1,M_MAX+1),B2(M_MAX+1,M_MAX+1)
      double precision x(0:M_MAX), fftarry(M_MAX+2),pi
      double precision pim, rhs(M_MAX+1), ff, exact
      double precision u(M_MAX+1),scl(M_MAX+1),c(IRE)
      double precision wk(M_MAX+1),w(M_MAX+1,IRE)
      double precision tol,a,err,ei

      common /fun/ a

      M = 64     
      a=5.d-4 

      pi = acos(-1.d0)
      pim = pi/dble(M)
      call fftini(M)
   
      do i = 0,M
         x(i) = cos(pim*dble(i))
      end do
      call integ1(B2,pre,M,M_MAX)


*** Transform the right hand side **************
      do i = 1,M+1
         fftarry(i) = ff(x(i-1))
      end do 
      call cosft(M,fftarry,um,vm,wm,ibitm,1)

****** Add the boundary conditions and condition right hand side with
****** preconditioner      
      rhs(1) = exact(1.d0) 
      rhs(2) =  exact(-1.d0) 

      do i = 3,M+1
         sum = 0.d0
         do j = i,M+1
            sum = sum + pre(i,j)*fftarry(j-2)
         end do
         rhs(i) = sum
      end do   
*******Call the broyden solver********
      do i=1,M+1
        u(i)=0.d0
        if (i.le.2) then
          scl(i)=1.d0
        else
          scl(i)=1.d0
        end if 
      end do
      itmax=IRE      
      tol=5.d-14 
      call broyden(ire,M+1,M_MAX+1,scl,u,c,r,w,wk,itmax,tol,rhs)
********  Get an exact solution to check subroutine **************
      do i = 1,M+1
         fftarry(i) = u(i)
      end do
      call cosft(M,fftarry,um,vm,wm,ibitm,-1)
      err=0.d0
      do i = 1,M+1
         ei=abs(fftarry(i)-exact(x(i-1)))
         if (ei.gt.err) err=ei
      end do   
      write (6,*)err,M

      end

      double precision function ff(x)
      double precision ff,x,a
      common /fun/ a
       ff = (4.d0*((x-.5d0)/a)*((x-.5d0)/a)-2.d0/a+
     &       sin(x))*exp(-(x-.5d0)*(x-.5d0)/a)/sqrt(a)  
      end 

      double precision function exact(x)
      double precision exact,x,a
      common /fun/ a
      exact = exp(-(x-.5d0)*(x-.5d0)/a)/sqrt(a)
      end

      subroutine integ1(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 



