c Codifferential method for DC functions 
c by Adil Bagirov 2014

      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x0(maxvar),x(maxvar)
      COMMON /citer/maxiter,niter,/cnf/nf,/cngrad/ngrad
      character*30 outfi/'results.txt'/
      open(40,file=outfi)
c=======================================================
c  Input data:
c  n        - number of variables
c=======================================================
      n = 2
c======================================================
c nbundle - is the size of bundle, the maximum number of
c subgradients or discrete gradients which can be
c calculated at each iteration
c======================================================
      nbundle = n+3
c======================================================
c ngradient = 0 if you use approximation to subgradient
c ngradient = 1 if you use exact gradients of the objective
c               and constraint functions
c======================================================
      ngradient=1
c=======================================================
c  Starting point
c=======================================================
      x0(1)=-1.2d+00
      x0(2)=1.0d+00
c=======================================================
c  Calling CODIF
c=======================================================
       call codif(n,x0,nbundle,ngradient,x,fvalue,cputime,finit)
      WRITE(40,41) finit
  41  FORMAT('The value of the objective at initial point:  ',f16.8)
      WRITE(40,42) fvalue
  42  FORMAT('The value of the objective at final point:    ',f16.8)
      WRITE(40,43) nf
  43  FORMAT('The total number of the objective evaluations:',i12)
      WRITE(40,44) ngrad
  44  FORMAT('The total number of gradient evaluations:     ',i12)
      WRITE(40,46) cputime
  46  FORMAT('The CPU time:                                 ',f10.4)
c=======================================================
      CLOSE(40)
      STOP
      END

c======================================================================
      subroutine func1(x,f1)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar)
      COMMON /cka1/ka1
c==================================================================
      f1=abs(x(1)-1.0d+00)+2.0d+02*dmax1(0.0d+00,abs(x(1))-x(2))
c=================================================================
      return
      end

c======================================================================
      subroutine func2(x,f2)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar)
      COMMON /cka2/ka2
c===============================================================
      f2=1.0d+02*(abs(x(1))-x(2))
c======================================================
      return
      end

c======================================================================
      subroutine gradient1(x,grad)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar),grad(maxvar)
      COMMON /cka1/ka1
c===============================================================
      if(x(1).ge.1.0d+00) then
        grad(1)=1.0d+00
       else
        grad(1)=-1.0d+00
      end if
      grad(2)=0.0d+00
      if(abs(x(1)).ge.x(2)) then
       grad(2)=grad(2)-2.0d+02
       if(x(1).ge.0.0d+00) then
         grad(1)=grad(1)+2.0d+02
        else
         grad(1)=grad(1)-2.0d+02
       end if
      end if      
c===============================================================
      return
      end

c======================================================================
      subroutine gradient2(x,grad)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar),grad(maxvar)
      COMMON /cka2/ka2
c===============================================================
      if(x(1).ge.0.0d+00) then
        grad(1)=1.0d+02
       else
        grad(1)=-1.0d+02
      end if
      grad(2)=-1.0d+02
c===============================================================
      return
      end

c======================================================================
c Without scaling
c======================================================================
      subroutine codif(nvar,x0,nbundle,ngradient,x,f,cputime,finit)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar),x0(maxvar)
      COMMON /csize/m,/citer/maxiter,niter,/cnf/nf,/cgrad/ngradient1
     1 ,/cngrad/ngrad
      call cpu_time(time1)
      m=nvar
      nf=0
      niter=0
      ngrad=0
      nfgrad=0
      ngradient1=ngradient
      maxiter=10000
      slinit=1.0d+00
      do i=1,m
        x(i)=x0(i)
      end do
      call fv(x,f)
      finit=f
      call optimum(x,nbundle,slinit)
      call fv(x,f)
      call cpu_time(time2)
      cputime=time2-time1
c===========================================================
      return
      end

c============================================================
      subroutine optimum(x,nbundle,slinit)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000,maxit=100000)
      double precision x(maxvar),x1(maxvar),g(maxvar),v(maxvar)
     1 ,w1(maxdg,maxvar),w2(maxdg,maxvar),prod(maxdg,maxdg),z(maxdg)
     2 ,fvalues(maxit),prod1(maxdg,maxdg),prod3(maxdg,maxdg)
     3 ,prod2(maxdg),v1(maxvar),v2(maxvar),v3(maxvar)
      INTEGER ij(maxdg),ind(maxdg)
      common /csize/m,/citer/maxiter,niter,/cij/ij,jvertex,/cz/z
     1 ,/ckmin/kmin,/cmaxmin/fmax,fmin
c====================================================================
      dist1=1.0d-07
      step0=-2.0d-01
      div=1.0d-01
      eps0=1.0d-07
      slmin=1.0d-10*slinit
      pwt=1.0d-07
      sdif=1.0d-05
      mturn=4
c====================================================================
      sl=slinit/div
      call fv(x,f2)
      fmax2=fmax
      fmin2=fmin
  1   sl=div*sl
      IF(sl.lt.slmin) return
      do i=1,m
       g(i)=1.0d+00/dsqrt(DBLE(m))
      end do
      nnew=0
c================================================================
   2  niter=niter+1
        IF(niter.gt.maxiter) RETURN
        nnew=nnew+1
        f1=f2
        fmax1=fmax2
        fmin1=fmin2
        fvalues(niter)=f1
c---------------------------------------------------------------
        if (nnew.gt.mturn) then
         mturn2=niter-mturn+1
         ratio1=(fvalues(mturn2)-f1)/(dabs(f1)+1.0d+00)
         IF(ratio1.LT.sdif) GO TO 1
        end if
        if (nnew.GE.(2*mturn)) then
         mturn2=niter-2*mturn+1
         ratio1=(fvalues(mturn2)-f1)/(dabs(f1)+1.0d+00)
         IF(ratio1.LT.(1.0d+01*sdif)) GO TO 1
        end if

        ndg1=0
        ndg2=0
        do i=1,maxdg
         ind(i)=1
        end do
c--------------------------------------------------------------
        do ndg=1,nbundle
            call dgrad(x,sl,g,v1,v2,f1,fmax1,fmin1,pwt)
            ndg10=ndg1
            ndg20=ndg2
            if (ndg.eq.1) then
             ndg1=ndg1+1
             ndg2=ndg2+1
             prod1(1,1)=0.0d+00
             do j=1,m+1
              w1(ndg1,j)=v1(j)
              prod1(1,1)=prod1(1,1)+v1(j)*v1(j)
             end do
             prod2(1)=0.0d+00
             prod3(1,1)=0.0d+00
             do j=1,m+1
              w2(ndg2,j)=v2(j)
              prod2(1)=prod2(1)+v2(j)*v2(j)
              prod3(1,1)=prod3(1,1)+v1(j)*v2(j)
             end do
            end if
c           to check whether hew subgradients are close to previous ones
            if (ndg.gt.1) then
             do i=1,ndg1
              d1=0.0d+00
              do j=1,m+1
               d1=d1+(w1(i,j)-v1(j))**2
              end do
              d1=dsqrt(d1)
              IF(d1.le.eps0) GO TO 11
             end do
             ndg1=ndg1+1
             dotprod=0.0d+00
             do j=1,m+1
              w1(ndg1,j)=v1(j)
              dotprod=dotprod+v1(j)*v1(j)
             end do
             do i=1,ndg1-1
              prod1(ndg,i)=0.0d+00
              do j=1,m+1
               prod1(ndg1,i)=prod1(ndg1,i)+w1(i,j)*v1(j)
              end do
              prod1(i,ndg1)=prod1(ndg1,i)
             end do
             prod1(ndg1,ndg1)=dotprod
             do i=1,ndg2
              prod3(ndg1,i)=0.0d+00
              do j=1,m+1
               prod3(ndg1,i)=prod3(ndg1,i)+w2(i,j)*v1(j)
              end do
             end do
11           CONTINUE

             do i=1,ndg2
              d2=0.0d+00
              do j=1,m+1
               d2=d2+(w2(i,j)-v2(j))**2
              end do
              d2=dsqrt(d2)
              IF(d2.le.eps0) GO TO 12
             end do
             ndg2=ndg2+1
             prod2(ndg2)=0.0d+00
             do j=1,m+1
              w2(ndg2,j)=v2(j)
              prod2(ndg2)=prod2(ndg2)+v2(j)*v2(j)
             end do
             do i=1,ndg1
              prod3(i,ndg2)=0.0d+00
              do j=1,m+1
               prod3(i,ndg2)=prod3(i,ndg2)+w1(i,j)*v2(j)
              end do
             end do
            end if

 12         CONTINUE
            IF((ndg10.eq.ndg1).AND.(ndg20.eq.ndg2)) GO TO 1

            ifin=0
            do i=1,ndg2
             ifin=ifin+ind(i)
            end do
            IF(ifin.eq.0) GO TO 1

            d7=0.0d+00
            do i=1,ndg2
             if (ind(i).eq.1) then
              do j=1,ndg1
               do k=1,ndg1
                prod(j,k)=prod1(j,k)-prod3(j,i)-prod3(k,i)+prod2(i)
               end do
              end do
              d5=prod(1,1)
              d6=prod(1,1)
              kmin=1
              do j=2,ndg1
               if (d5.gt.prod(j,j)) then
                d5=prod(j,j)
                kmin=j
               end if
               if (d6.lt.prod(j,j)) then
                d6=prod(j,j)
               end if
              end do
              toler=dmax1(eps0,dist1*d6)
              call wolfe(ndg1,prod)
              do k=1,m+1
               v3(k)=0.0d+00
               do k1=1,jvertex
                v3(k)=v3(k)+(w1(ij(k1),k)-w2(i,k))*z(k1)
               END do
              END do
              r=0.0d+00
              do k=1,m+1
               r=r+v3(k)*v3(k)
              end do
              r=dsqrt(r)
              if(r.lt.toler) ind(i)=0
              if (r.ge.toler) then
               if (d7.lt.r) then
                d7=r
                do k=1,m
                 g(k)=-v3(k)/d7
                end do
               end if
              end if
             end if
            end do
            IF(d7.le.eps0) GO TO 1
            r=d7
c===========================================================
             do i=1,m
              x1(i)=x(i)+sl*g(i)
             end do
c===========================================================
             call fv(x1,f4)
             f3=(f4-f1)/sl
             decreas=step0*r
             if(f3.lt.decreas) then
                        call armijo(x,g,f1,f5,f5max,f5min,f4,sl,step,r)
                        f2=f5
                        fmax2=f5max
                        fmin2=f5min
                        do i=1,m
                         x(i)=x(i)+step*g(i)
                        end do
                        print 3, niter,ndg,f2,sl,r
                        sl=1.2d+00*sl
                        GO TO 2
 3                      format(I6,I4,3f16.8)
             end if
         END do
c=====================================================
      go to 1
      return
      end

c==============================================================
c  Subroutines Wolfe and Equations solves quadratic
c  programming problem, to find
c  descent direction, Step 3, Algorithm 2.
c===============================================================

      subroutine wolfe(ndg,prod)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000)
      common /csize/m,/w01/a,/cij/ij,jvertex,/cz/z,/ckmin/kmin
      INTEGER ij(maxdg)
      double precision z(maxdg),z1(maxdg),a(maxdg,maxdg)
     1 ,prod(maxdg,maxdg)
      j9=0
      jmax=500*ndg
      jvertex=1
      ij(1)=kmin
      z(1)=1.0d+00
c=======================================
c  To calculate X
c=======================================
 1    r=0.0d+00
      do i=1,jvertex
       do j=1,jvertex
        r=r+z(i)*z(j)*prod(ij(i),ij(j))
       end do
      end do
      IF(ndg.eq.1) RETURN
c========================================
c  To calculate <X,P_J> and J
c========================================
      t0=1.0d+12
      do i=1,ndg
        t1=0.0d+00
        do j=1,jvertex
          t1=t1+z(j)*prod(ij(j),i)
        end do
        if(t1.lt.t0) then
                     t0=t1
                     kmax=i
        end if
      end do
c========================================
c  First stopping criterion
c========================================
      rm=prod(kmax,kmax)
      do j=1,jvertex
       rm=dmax1(rm,prod(ij(j),ij(j)))
      end do
      r2=r-1.0d-12*rm
      if(t0.gt.r2) RETURN
c========================================
c  Second stopping criterion
c========================================
      do i=1,jvertex
       if(kmax.eq.ij(i)) RETURN
      end do
c========================================
c Step 1(e) from Wolfe's algorithm
c========================================
      jvertex=jvertex+1
      ij(jvertex)=kmax
      z(jvertex)=0.0d+00
c========================================
 2    do i=1,jvertex
       do j=1,jvertex
        a(i,j)=1.0d+00+prod(ij(i),ij(j))
       end do
      end do
      j9=j9+1
      if(j9.gt.jmax) RETURN
      call equations(jvertex,z1)
      do i=1,jvertex
       if(z1(i).le.1.0d-10) go to 3
      end do
      do i=1,jvertex
       z(i)=z1(i)
      end do
      go to 1
  3   teta=1.0d+00
      do i=1,jvertex
       z5=z(i)-z1(i)
       if(z5.gt.1.0d-10) teta=dmin1(teta,z(i)/z5)
      end do
      do i=1,jvertex
       z(i)=(1.0d+00-teta)*z(i)+teta*z1(i)
       if(z(i).le.1.0d-10) then
                          z(i)=0.0d+00
                          kzero=i
       end if
      end do
      j2=0
      do i=1,jvertex
       IF(i.ne.kzero) then
                     j2=j2+1
                     ij(j2)=ij(i)
                     z(j2)=z(i)
       END if
      end do
      jvertex=j2
      go to 2
      return
      end

      subroutine equations(n,z1)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000)
      common /w01/a
      double precision a(maxdg,maxdg),z1(maxdg),b(maxdg,maxdg)
      do i=1,n
       do j=1,n
        b(i,j)=a(i,j)
       end do
       b(i,n+1)=1.0d+00
      end do
      do i=1,n
       r=b(i,i)
       do j=i,n+1
        b(i,j)=b(i,j)/r
       end do
       do j=i+1,n
        do k=i+1,n+1
         b(j,k)=b(j,k)-b(i,k)*b(j,i)
        end do
       end do
      end do
      z1(n)=b(n,n+1)
      do i=1,n-1
        k=n-i
        z1(k)=b(k,n+1)
        do j=k+1,n
         z1(k)=z1(k)-b(k,j)*z1(j)
        END do
      end do
      z2=0.0d+00
      do i=1,n
       z2=z2+z1(i)
      end do
      do i=1,n
       z1(i)=z1(i)/z2
      end do
      return
      end

c=====================================================================
c Subroutine dgrad calculates subgradients or discrete gradients
c=====================================================================
      subroutine dgrad(x,sl,g,dg1,dg2,f1,fmax1,fmin1,pwt)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000)
      double precision x1(maxvar),x2(maxvar),g(maxvar),x(maxvar)
     1 ,dg1(maxvar),dg2(maxvar)
      common /csize/m,/cgrad/ngradient1,/cngrad/ngrad

      sl2=1.0d-08
      f10=fmax1
      f20=fmin1
      do k=1,m
        x1(k)=x(k)+sl*g(k)
        x2(k)=x(k)+sl2*g(k)
      end do

      call fvmax(x1,r2)
      call fvmin(x2,r3)

      if (ngradient1.eq.0) then
       call dgrad2(x1,x2,dg1,dg2,r2,r3,pwt)
       ngrad=ngrad+1
      end if

      if (ngradient1.eq.1) then
       call gradient1(x1,dg1)
       call gradient2(x2,dg2)
       ngrad=ngrad+1
      end if

      d1=0.0d+00
      d2=0.0d+00
      do i=1,m
       d1=d1+dg1(i)*g(i)
       d2=d2+dg2(i)*g(i)
      end do
      dg1(m+1)=r2-f10-sl*d1
      dg2(m+1)=r3-f20-sl2*d2

      return
      end

c=====================================================================
c Subroutine dgrad calculates discrete gradients: Step 5, Algorithm 1
c=====================================================================
      subroutine dgrad2(x1,x2,v1,v2,r2,r3,pwt)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000)
      double precision x1(maxvar),x2(maxvar),v1(maxvar),v2(maxvar)
      common /csize/m
      alpha=1.0d+00
      t=pwt/alpha
      r4=r2
      r5=r3
      do k=1,m
        t=t*alpha
        r6=r4
        r7=r5
        x1(k)=x1(k)+t
        x2(k)=x2(k)+t
        call fvmax(x1,r4)
        call fvmin(x2,r5)
        nfgrad=nfgrad+1
        v1(k)=(r4-r6)/t
        v2(k)=(r5-r7)/t
      end do
      return
      end

c===========================================================
c Line search (Armijo-type), Step 5 Algorithm 2.
c===========================================================
      subroutine armijo(x,g,f1,f5,f5max,f5min,f4,sl,step,r)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000)
      common /csize/m,/cmaxmin/fmax,fmin
      double precision x(maxvar),g(maxvar),x1(maxvar)
      step=sl
      f5=f4
      f5max=fmax
      f5min=fmin
      step1=sl
      k=0
  1   k=k+1
      IF(k.gt.20) RETURN
      step1=2.0d+00*step1
      do i=1,m
       x1(i)=x(i)+step1*g(i)
      end do
      call fv(x1,f50)
      f30=f50-f1+5.0d-02*step1*r
      IF(f30.gt.0.0d+00) RETURN
      step=step1
      f5=f50
      f5max=fmax
      f5min=fmin
      GO TO 1
      return
      end

      subroutine fv(x,f)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar)
      COMMON /cnf/nf,/csize/m,/cmaxmin/fmax,fmin
c================================================================
      call func1(x,fmax)
      call func2(x,fmin)
      nf=nf+1
      f=fmax-fmin
c=================================================================
      return
      end

      subroutine fvmax(x,f)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar)
      COMMON /cnf/nf
c================================================================
      call func1(x,objf)
      nf=nf+1
      f=objf
c=================================================================
      return
      end

      subroutine fvmin(x,f)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar)
      COMMON /cnf/nf
c================================================================
      call func2(x,objf)
      f=objf
c=================================================================
      return
      end