c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      subroutine broyden(rhot,rhin,xvalws,
     <                   iter,alpha,beta,komp,nbasis,
     <                       mspn,mtasa,jtop,rms)
      implicit real*8 (a-h,o-z)
c************************************************************
c Written by D.D. Johnson (from Phys. Rev. B38, 12807 (1988))
c note: A few typo's exist in that paper involving weights!
c       Code below is correct.
c************************************************************
c      double precision b,d on machines with low accuracy
c      such as old VAXes, etc.
c************************************************************
c iter  = iteration number
c alpha = empirical mix parameter, density (about .1 -> .3)
c         sometimes reduce alpha with large number of atoms.
c beta  = emperical mix parameter, moments (about .5 -> .8)
c komp  = number of components in alloy
c nbasis= number of sublattices in alloy
c mspn  = number of electron spin components (1 or 2)
c mtasa = Is this MT or ASA calculation?     (0=MT or 1=ASA)
c jtop  = size of vector to be used in mixing
c rms   = chg. density rms deviation used for weighting
c
c For LDA use, charge density mixing with total rho and
c magnetization (i.e., difference rho) is found optimal.
c As such, simple  mix of magnetization and 
c          Broyden mix of densities.
c
c************************************************************
c*  the vectors ui(maxsiz) and vti(maxsiz) are johnson's    *
c*  u(of i)  and df(transpose), respectively. these are     *
c*  continually updated. all iterations are stored on tape  *
c*  32 . this is done to prevent the prohibitive storage    *
c*  costs associated with holding onto the entire jacobian. *
c*  vector ul is the vt of earlier iterations. vector f is: *
c*  vector(output) - vector(in). vector df is:  f(m+1)-f(m) *
c*  finally,vector dumvi(maxsiz) is the predicted vector.   *
c************************************************************
c
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      include 'mkkrcpa.h'
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer    maxsiz,imatsz
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      parameter (maxsiz=ipsublat*ipcomp*iprpts)
      parameter (imatsz=40)
c++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
c***** arrays strictly for Broyden construction: do not alter *****
      dimension a(imatsz,imatsz),b(imatsz,imatsz),cm(imatsz)
      dimension d(imatsz,imatsz),w(imatsz)
      dimension f(maxsiz),ui(maxsiz),vti(maxsiz),t1(maxsiz)
      dimension vector(maxsiz,2),dumvi(maxsiz),df(maxsiz)
c
c==============  program dependent arrays for load/unload ==========
cab      character*30 obroy1,obroy2
      character*50 obroy1,obroy2
      common/broydn/obroy1,obroy2
c
      dimension   rhin(iprpts,ipcomp,ipsublat,ipspin)
      dimension   rhot(iprpts,ipcomp,ipsublat,ipspin)
      dimension   xvalws(ipcomp,ipsublat,ipspin)
c
      integer     komp(ipsublat)
c===================================================================
c
c++++++ files for BROYDEN mixing ,                            ++++++
c      files can be deleted any time for improved convergence ++++++
      open(61,file=obroy1,status='unknown',form='unformatted')
      open(62,file=obroy2,status='unknown',form='unformatted')
c
c++++++ set up the vector of the current iteration for mixing ++++++
c
c  For this method we have only saved input/output chg. densities.
c  For MT/ASA also include XVALWS in mix in jtop+1 array location.
      jtop1=jtop+1
c  Keep track of basis indexing
       indx=0
      do ib=1,nbasis
       do nc=1,komp(ib)
           indx=indx+1
           is=(indx-1)*jtop1
c
         do k=1,jtop1
           vector(k+is,1)=rhin(k,nc,ib,1) + (mspn-1)*rhin(k,nc,ib,mspn)
           vector(k+is,2)=rhot(k,nc,ib,1) + (mspn-1)*rhot(k,nc,ib,mspn)
         enddo
       enddo
      enddo
c  ivsiz is the length of the vector
      ivsiz=is + jtop1
      if(iter.lt.3)write( 6,1001)ivsiz
      if(ivsiz.gt.maxsiz)then
         write(6,'('' in mixing: exceeded max. vector length'')')
         stop
      endif
c++++++ end of program specific loading of vector from main ++++++++
c
c
c*******************  begin broyden's method  **********************
c
c  weighting factor for the zeroth iteration
      w0=.01       ! w0**2 must be larger than allowed precision
c
c      f:  the difference of previous output and input vectors
c  dumvi:  a dummy vector, here it is the previous input vector
      rewind 61
      read(61,err=119,end=119)amix,lastit
      read(61)(f(k),k=1,ivsiz)
      read(61)(dumvi(k),k=1,ivsiz)
      if(iter.eq.1 .and. lastit.gt.1)then
      read(61)ltmp,((a(i,j),i=1,ltmp),j=1,ltmp)
      read(61)(w(i),i=1,ltmp)
      endif
c
c  alpha (or amix, when saved) and beta are previously chosen
c  simple mixing parameters for rho and mag.
      write(6,1002)amix,beta,iter+1
c
      do 104 k=1,ivsiz
      dumvi(k)=vector(k,1)-dumvi(k)
  104 df(k)=vector(k,2)-vector(k,1)-f(k)
      do 114 k=1,ivsiz
  114 f(k)=vector(k,2)-vector(k,1)
c
c  for i-th iter.,dfnorm is ( f(i) minus f(i-1) ), used for normalization.
      dfnorm=zero
      fnorm=zero
      do 113 k=1,ivsiz
      dfnorm=dfnorm + df(k)*df(k)
  113 fnorm=fnorm + f(k)*f(k)
      dfnorm=sqrt(dfnorm)
      fnorm=sqrt(fnorm)
      write(6,'(''  dfnorm '',e12.6,'' fnorm '',e12.6)')dfnorm,fnorm
c
      fac2=one/dfnorm
      fac1=amix*fac2
      do 105 k=1,ivsiz
      ui(k) = fac1*df(k) + fac2*dumvi(k)
 105  vti(k)= fac2*df(k)
c
c*********** calculation of coefficient matrices *************
c***********    and the sum for corrections      *************
c
c recall: a(i,j) is a symmetric matrix
c       : b(i,j) is the inverse of [ w0**2 I + a ]
c
         lastit=lastit+1
         lastm1=lastit-1
         lastm2=lastit-2
c
c dumvi is the u(of i) and t1 is the vt(of i)
c from the previous iterations
      rewind 62
      write(6,1003)lastit,lastm1
      if(lastit.gt.2)then
      do 500 j=1,lastm2
      read(62)(dumvi(k),k=1,ivsiz)
      read(62)(t1(k),k=1,ivsiz)
c
      aij=zero
      cmj=zero
      do 501 k=1,ivsiz
      cmj=cmj + t1(k)*f(k)
  501 aij=aij + t1(k)*vti(k)
      a(lastm1,j)=aij
      a(j,lastm1)=aij
            cm(j)=cmj
  500 continue
      endif
c
      aij=zero
      cmj=zero
      do 106 k=1,ivsiz
      cmj= cmj + vti(k)*f(k)
  106 aij= aij + vti(k)*vti(k)
      a(lastm1,lastm1)=aij
            cm(lastm1)=cmj
c
      write(62)(ui(k),k=1,ivsiz)
      write(62)(vti(k),k=1,ivsiz)
      rewind 62
c
c the weighting factors for each iteration have been chosen
c equal to one over the r.m.s. error. This need not be the case.
c       wtmp=0.010D0/fnorm
       wtmp=0.0d0
       if(rms.gt. 1.0d-09)wtmp=two*sqrt(0.010D0/rms)
       if(wtmp.lt. one)wtmp=one
       w(lastm1)=wtmp
       write(6,'(''  WEIGHTING SET =  '',e12.6)')wtmp
c
c
c with the current iterations f and vector calculated,
c write them to unit 61 for use later.
      rewind 61
      write(61)amix,lastit
      write(61)(f(k),k=1,ivsiz)
      write(61)(vector(k,1),k=1,ivsiz)
      write(61)lastm1,((a(i,j),i=1,lastm1),j=1,lastm1)
      write(61)(w(i),i=1,lastm1)
c
c set up and calculate beta matrix
      do 506 lm=1,lastm1
      do 507 ln=1,lastm1
         d(ln,lm)= a(ln,lm)*w(ln)*w(lm)
 507     b(ln,lm)= zero
         b(lm,lm)= one
 506     d(lm,lm)= w0**2 + a(lm,lm)*w(lm)*w(lm)
c
      call invert(d,b,lastm1)
c
c  calculate the vector for the new iteration
      do 505 k=1,ivsiz
  505 dumvi(k)= vector(k,1) + amix*f(k)
c
      do 504 i=1,lastm1
      read(62)(ui(k),k=1,ivsiz)
      read(62)(vti(k),k=1,ivsiz)
      gmi=zero
      do 503 ip=1,lastm1
  503 gmi=gmi + cm(ip)*b(ip,i)*w(ip)
      do 504 k=1,ivsiz
  504 dumvi(k)=dumvi(k)-gmi*ui(k)*w(i)
c  end of the calculation of dumvi, the new vector
c
      rewind 61
      rewind 62
c
      goto 120
c if this is the first iteration, then load
c    f=vector(out)-vector(in) and vector(in)
  119 rewind 61
      lastit=1
c  alpha and beta are simple mixing parameters for rho and mag.
      write(6,1002)alpha,beta,lastit
      amix=alpha
      write(61)amix,lastit
      do 101 k=1,ivsiz
  101 f(k)=vector(k,2)-vector(k,1)
      write(61)(f(k),k=1,ivsiz)
      write(61)(vector(k,1),k=1,ivsiz)
c
c since we are on the first iteration, simple mix the vector.
      do 102 k=1,ivsiz
  102 dumvi(k)= vector(k,1) + amix*f(k)
      write( 6,1000)
  120 continue
c
      close(61,status='keep')
      close(62,status='keep')
c
c*************  the end of the broyden method **************
c
c+++++++ program specific code of reloading arrays +++++++++
c
c need to unload the new vector into the appropriate arrays.
      xx=one-beta
      nspm1=mspn-1
      sp=one/mspn
c
c ************************** simple mix moments, then unload vector!!!!
         indx=0
      do ib=1,nbasis
        do nc=1,komp(ib)
         indx=indx+1
         is=(indx-1)*jtop1
c **************** unload new XVALWS
             a1= dumvi(jtop1+is)
             b1= xx*( rhin(jtop1,nc,ib,1)-rhin(jtop1,nc,ib,mspn) ) + 
     <         beta*( rhot(jtop1,nc,ib,1)-rhot(jtop1,nc,ib,mspn) )
              xvalws(nc,ib,mspn)=sp*nspm1*( a1 - b1 )
              xvalws(nc,ib,   1)=sp*( a1 + b1*nspm1 )
c **************** finished unload of XVALWS
          do k=1,jtop
             a1= dumvi(k+is)
             b1= xx*( rhin(k,nc,ib,1)-rhin(k,nc,ib,mspn) ) + 
     <         beta*( rhot(k,nc,ib,1)-rhot(k,nc,ib,mspn) )
              rhot(k,nc,ib,mspn)=sp*nspm1*( a1 - b1 )
              rhot(k,nc,ib,   1)=sp*( a1 + b1*nspm1 )
          enddo
        enddo
      enddo
c ************************** end simple mix of moments, unloading!!!!
c
c+++++++++ end of program specific reloading of arrays +++++++++
c
      write(6,1004)iter+1
      return
c
 1000 format(' ---->  straight mixing on this iteration')
 1001 format(' in mixing:   ivsiz =',i7,/)
 1002 format(' in mixing: simple mix parameters',2(f10.6,',')
     >      ,'  for iter=',i5)
 1003 format(' current iter= ',i5,' includes values from iter=',i5)
 1004 format(10x,'density for iteration',i4,' prepared')
c
      end
c     ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      subroutine invert(a,b,m)
c     =============================================================
      implicit real*8  (a-h,o-z)
c      complex a,b,td,ad,bd
c
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      parameter (imatsz=40)
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c
      dimension a(imatsz,imatsz),b(imatsz,imatsz)
      dimension td(imatsz),ad(imatsz),bd(imatsz)
c
c subroutine to preform gaussian elimination
c            no zeros along the diagonal
      save
c
      n=m
      if(n.gt.imatsz)then
       write(6,'('' invert: matrix a too large'')')
       stop
      endif
c
      do 14 i=1,n
      atmp=a(i,i)
      if(abs(atmp).lt. 1.0d-08)then
        write(6,'('' invert: matrix has zero diagonal'',
     <            '' element in the '',i4,'' row'')')i
        stop
      endif
  14  continue
c
      if(n.eq.1) go to 605
c
      do 23 i=1,n
c
      do 35 j=1,n
 35      td(j)=a(j,i)/a(i,i)
c
      td(i)=(0.0e+00,0.0e+00)
c
      do 71 k=1,n
         bd(k)=b(i,k)
 71      ad(k)=a(i,k)
c
      do 601 k=1,n
      do 601 j=1,n
         b(j,k)=b(j,k)-(td(j)*bd(k))
 601     a(j,k)=a(j,k)-(td(j)*ad(k))
c
 23   continue
c
      do 603 i=1,n
      do 603 j=1,n
 603     b(j,i)=b(j,i)/a(j,j)
c
      return
c
 605  b(1,1)=1.0e+00/a(1,1)
      return
      end
c     ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc