#undef USEMODULE

#ifdef USEMODULE
      module m_mixing_simple
      implicit none
      interface mixing_simple
      module procedure mixing_simple_r, mixing_simple_c
      end interface
      contains
#endif

      subroutine mixing_simple_r(ite,n,alpha,rhoin,fin,pout)
      implicit none
      integer::ite,n
      real(8):: alpha
      real(8),intent(in)::fin(n),rhoin(n)
      real(8),intent(out):: pout(n)

      integer::i

      do i=1,n
         pout(i)=rhoin(i)*(1.0d0-alpha)+alpha*fin(i)
      enddo

      end subroutine mixing_simple_r

      subroutine mixing_simple_c(ite,n,alpha,rhoin,fin,pout)
      implicit none
      integer::ite,n
      real(8):: alpha
      complex(8),intent(in)::fin(n),rhoin(n)
      complex(8),intent(out):: pout(n)

      integer::i

      do i=1,n
         pout(i)=rhoin(i)*(1.0d0-alpha)+alpha*fin(i)
      enddo

      end subroutine mixing_simple_c


#ifdef USEMODULE
      end module m_mixing_simple
#endif

#ifdef USEMODULE 
      module m_mixing_diis
      implicit none
      interface mixing_diis3
      module procedure mixing_diis3_r,mixing_diis3_c
      end interface 
      contains
#endif

c------------------

c     Anderson mixing or DIIS mixing level 3
c     input ite,n,pin,a
c     output  Fout
c     work Fwork
      subroutine mixing_diis3_r(ite,n,alpha, fin, pout, Fwork,nstepFwork)
      implicit none
      integer,intent(in)::ite,n,nstepFwork
      real(8),intent(in):: Fin(n)
      real(8),intent(out):: pout(n)
      real(8),intent(out):: Fwork(nstepFwork,3) ! 3=2*level-1
      real(8),intent(in) :: alpha !mixing
c     Fwork(n,1) : F(i-1)
c     Fwork(n,2) : p(i-1)
c     Fwork(n,3) : p(i)

      integer::i
      real(8) :: deno,nume,varFn,varFn1,dFn,theta


      if (ite==1) then
c     Fout = pin
         do i=1,n
            pout(i) = fin(i)
         enddo
         do i=1,n
            Fwork(i,3)=pout(i)
         enddo
      else if (ite==2) then
c     Fout(2) = Fout(1)*a+pin*(1-a)
         do i=1,n
            pout(i)= Fwork(i,3)*(1-alpha)+fin(i)*alpha
         enddo
         do i=1,n
            Fwork(i,1)=fin(i)
         enddo
         do i=1,n
            Fwork(i,2)=Fwork(i,3)
         enddo
         do i=1,n
            Fwork(i,3)=pout(i)
         enddo 
      else
#if 0
         write(*,*) ite,'p(n-1)=',Fwork(:,2)
         write(*,*) ite,'p(n)  =',Fwork(:,3)
         write(*,*) ite,'F(n-1)=',Fwork(:,1)
         write(*,*) ite,'F(n)  =',fin(:)
#endif
         nume=0.0d0
         deno=0.0d0
         do i=1,n
c     varF(n) =  F(n)-p(n)
c     dF(n) =  varF(n)-varF(n-1)
            varFn =  Fin(i)-Fwork(i,3)
            varFn1 = Fwork(i,1)-Fwork(i,2)
            dFn = varFn -varFn1
c     numerator = <dFn | varF(n-1) >
c     demoninator=<dF(n)|dF(n)>
            nume = nume + dFn*varFn1
            deno = deno + dFn*dFn
         enddo
         theta = -nume/deno
c     if (theta<0.0) theta=0.0
c     if (theta>1.0) theta=1.0
#if 0
         write(*,'(a,i,f10.5)') 'ite,theta=',ite,theta
#endif
         do i=1,n
c     p(n+1) = (1-alpha)*varrho(n) + alpha*varF(n)
c     varrho(n) = (1-theta)*rho(n-1)+theta*rho(n)
c     varF(n)   = (1-theta)*F(n-1)+theta*F(n)
            pout(i) = (1.0d0-alpha)*( (1.0d0-theta)*Fwork(i,2) + theta*Fwork(i,3) )
     .           +  alpha*( (1.0d0-theta)*Fwork(i,1)+theta*Fin(i) )  
         enddo
         do i=1,n
            Fwork(i,1) = Fin(i)
         enddo
         do i=1,n
            Fwork(i,2) = Fwork(i,3)
         enddo
         do i=1,n
            Fwork(i,3) = pout(i)
         enddo
      endif

      end subroutine mixing_diis3_r

c-----------------------------------------

c     Anderson mixing or DIIS mixing level 3
c     input ite,n,pin,a
c     output  Fout
c     work Fwork
      subroutine mixing_diis3_c(ite,n,alpha, fin, pout, Fwork,nstepFwork)
      implicit none
      integer,intent(in)::ite,n,nstepFwork
      complex(8),intent(in):: Fin(n)
      complex(8),intent(out):: pout(n)
      complex(8),intent(out):: Fwork(nstepFwork,3) ! 3=2*level-1
      real(8),intent(in) :: alpha !mixing
c     Fwork(n,1) : F(i-1)
c     Fwork(n,2) : p(i-1)
c     Fwork(n,3) : p(i)

      integer::i
      real(8) :: deno,nume,theta
      complex(8)::varFn,varFn1,dFn

      if (n>nstepFwork) then
         write(*,*) 'n>nstepFwork n=',n,'nstepFwork=',nstepFwork
         stop 'in mixing_diis3_c'
      endif

c        write(*,*)'ite,n,alpha=',ite,n,alpha


      if (ite==1) then
c         write(*,*) 'diis3 mixing, only set workarea'
c     Fout = pin
         do i=1,n
            pout(i) = fin(i)
         enddo
         do i=1,n
            Fwork(i,3)=pout(i)
         enddo
      else if (ite==2) then
c         write(*,*)'diis3 mixing, adopt simple mixing'
c     Fout(2) = Fout(1)*a+pin*(1-a)
         do i=1,n
            pout(i)= Fwork(i,3)*(1-alpha)+fin(i)*alpha
         enddo
         do i=1,n
            Fwork(i,1)=fin(i)
         enddo
         do i=1,n
            Fwork(i,2)=Fwork(i,3)
         enddo
         do i=1,n
            Fwork(i,3)=pout(i)
         enddo 
      else
#if 0
         write(*,*) ite,'p(n-1)=',Fwork(:,2)
         write(*,*) ite,'p(n)  =',Fwork(:,3)
         write(*,*) ite,'F(n-1)=',Fwork(:,1)
         write(*,*) ite,'F(n)  =',fin(:)
#endif
         nume=0.0d0
         deno=0.0d0
         do i=1,n
c     varF(n) =  F(n)-p(n)
c     dF(n) =  varF(n)-varF(n-1)
            varFn =  Fin(i)-Fwork(i,3)
            varFn1 = Fwork(i,1)-Fwork(i,2)
            dFn = varFn -varFn1
c     numerator = <dFn | varF(n-1) >
c     demoninator=<dF(n)|dF(n)>
            nume = nume + dreal(dFn)*dreal(varFn1)+dimag(dFn)*dimag(varFn1)
            deno = deno + dreal(dFn)*dreal(dFn)+dimag(dFn)*dimag(dFn)
         enddo
         theta = -nume/deno
c     if (theta<0.0) theta=0.0
c     if (theta>1.0) theta=1.0
#if 1
         write(*,*) 'ite,theta=',ite,theta,' numerator/denominator=',nume,deno
#endif
         do i=1,n
c     p(n+1) = (1-alpha)*varrho(n) + alpha*varF(n)
c     varrho(n) = (1-theta)*rho(n-1)+theta*rho(n)
c     varF(n)   = (1-theta)*F(n-1)+theta*F(n)
            pout(i) = (1.0d0-alpha)*( (1.0d0-theta)*Fwork(i,2) + theta*Fwork(i,3) )
     .           +  alpha*( (1.0d0-theta)*Fwork(i,1)+theta*Fin(i) )  
         enddo
         do i=1,n
            Fwork(i,1) = Fin(i)
         enddo
         do i=1,n
            Fwork(i,2) = Fwork(i,3)
         enddo
         do i=1,n
            Fwork(i,3) = pout(i)
         enddo
      endif

      end subroutine mixing_diis3_c

#ifdef USEMODULE
      end module m_mixing_diis
#endif


c=================================

#if 0
c a sample program
#define  NSIZE   4
      program test
      use m_mixing_diis
      use m_mixing_simple

      integer,parameter:: size=NSIZE/2

      real(8):: fact(NSIZE)
      common/deriv/fact

      complex(8):: p(size),p2(size),f(size),fwork(size,3)
      complex(8):: df(size)

      integer::ite,maxite

      real(8)::fval
      real(8):: alpha,dlen

      integer::mixing_flag=2
c     =1  simple
c     =2  DIIS3

      write(*,*) 'alpha='
      read(*,*) alpha

      maxite=1000

      do i=1,NSIZE
         fact(i)=i
      enddo

      p(:)=-10.0
      p(1)=10.0
      p(2)=15.0

      ite=1

      if (mixing_flag==2) then
         p2=p
         call mixing_diis3(ite,size,alpha,p2,p,fwork,size)
         write(*,100)'err=',ite,dlen(size,p-p2),dlen(size,p)
         ite=ite+1
      endif

 100  format(a,i5,2e20.10,10f20.9)


      do while (ite<=maxite)

         fval = func(size,p)
         call dfunc(size,p,df)

         do i=1,size
            p2(i) = p(i)-fval/df(i)
         enddo

         if (mixing_flag==1) then
            call mixing_simple(ite,size,alpha,p,p2,p2)
            write(*,100) 'err=',ite,dlen(size,p2-p),dlen(size,p2),p2
            p=p2
         else if (mixing_flag==2) then
            call mixing_diis3(ite,size,alpha,p2,p,fwork,size)
            write(*,100) 'err=', ite,dlen(size,p-p2),dlen(size,p),p
         endif

         if (dlen(size,p)<1.0d-9) stop

         ite=ite+1
      enddo

      end

c------------
      real(8) function func(n,p)
      implicit none
      integer::n
      real(8):: p(NSIZE)

      real(8):: fact(NSIZE)
      common/deriv/fact

      integer::i
      func=0.0d0
      do i=1,NSIZE
         func = func+ fact(i)*p(i)**2
      enddo

      end  function func

c----------------
      subroutine dfunc(n,p,df)
      implicit none
      integer::n
      real(8):: p(NSIZE)
      real(8):: df(NSIZE)

      real(8):: fact(NSIZE)
      common /deriv/fact

      integer::i 

      do i=1,NSIZE
         df(i) = 2.0*fact(i)*p(i)
      enddo
      end subroutine dfunc

c-----------------

      real(8) function dlen(n,p)
      integer::i,n
      real(8)::p(NSIZE)
      dlen=0.0d0
      do i=1,NSIZE
         dlen=dlen+p(i)**2
      enddo
      dlen=dsqrt(dlen)
      end function dlen

#endif



