C-------------------------------------------------------------------
      subroutine broyden(ire,neqs,ndim,scl,x,c,r,w,wk,itmax,tol,rhs)
C-------------------------------------------------------------------
      implicit double precision (a-h,o-z)
      dimension x(ndim),scl(ndim),c(ire),r(ndim),w(ndim,ire)
      dimension wk(ndim),rhs(ndim)
c
c     this program solves r(x)=0 using Broyden's method 
c     we follow the incremental storage approach of Kelley
c     here x is the solution, scl is a scaling which defines
c     the inner product, ire is a restart parameter, tol is
c     the tolerance and itmax the maximum number of iterations
c
c     first we initialize
c
      itct=0
      ir=0
      small=1.d-2
c
c  the main loop
c
 100  continue
      itct=itct+1
      ir=ir+1
      icase=2
      if (itct.eq.1) icase=1
      call resisub(x,rhs,r,neqs-1,icase)
      romo=sqrt(dsp(neqs,r,r,scl))
      write (6,*)itct,romo
      if (romo.le.tol) then
        go to 400
      end if 
      if (itct.gt.itmax) then
        write (6,*)' no convergence '
        stop
      end if 
c
c  apply H_(k-1)
c
      if (ir.gt.2) then
        do 120 ip=1,ir-2
          do 105 i=1,neqs
            wk(i)=w(i,ip)
 105      continue
          str=dsp(neqs,r,wk,scl)
          fact=str/c(ip)
          do 110 i=1,neqs
            r(i)=r(i)+fact*w(i,ip+1)
 110      continue
 120    continue
      end if
      if (ir.gt.1) then 
        do 125 i=1,neqs
          wk(i)=w(i,ir-1)
 125    continue
        stz=dsp(neqs,r,wk,scl)
        sntst=abs(1.d0-(stz/c(ir-1)))
        if (sntst.lt.small) then
          write (6,*)' stability correction ',itct
          sfact=.5d0
        else
          sfact=1.d0
        end if
        fact=c(ir-1)/(c(ir-1)-sfact*stz) 
        do 130 i=1,neqs
          r(i)=fact*r(i)
 130    continue
      end if
      do 140 i=1,neqs
        x(i)=x(i)+r(i)
 140  continue
      if (ir.eq.1) sfact=1.d0 
      c(ir)=dsp(neqs,r,r,scl)/sfact
      do 150 i=1,neqs
        w(i,ir)=r(i)
 150  continue
c
c   this completes the main iteration loop
c
      if (ir.gt.ire) ir=0
      go to 100
 400  continue
      do 420 i=1,neqs
        x(i)=x(i)+r(i)
 420  continue
      return
      end
c
      double precision function dsp(neqs,v,w,scl)
      implicit double precision (a-h,o-z)
      dimension v(*),w(*),scl(*)
c
c  this routine computes the scaled dot product of v and w
c
      dsp=0.d0
      do 100 i=1,neqs
        dsp=dsp+v(i)*w(i)*scl(i)
 100  continue
      return
      end 
  

