c The code of Nonsmooth DCA method 
c by Adil Bagirov 2014

      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x0(maxvar),x(maxvar)
      COMMON /cnf/nf,/cngrad/ngrad,/cndca/ndca
      character*30 outfi/'results.txt'/
      open(40,file=outfi)
c=======================================================
c  n        - number of variables
c=======================================================
      n  = 2
c======================================================
c nbundle - is the size of bundle
c======================================================
      nbundle=n+3
c======================================================
c  Starting point
c=======================================================
      x0(1)=-1.2d+00
      x0(2)=1.0d+00
c=======================================================
c  Calling DCA
c=======================================================
      call DCA(n,x0,nbundle,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,45) ndca
  45  FORMAT('The total number of DCA cycles:               ',i12)
      WRITE(40,46) cputime
  46  FORMAT('The CPU time:                                 ',f10.4)
      WRITE(40,47) x 
  47  FORMAT('The solution:                                 ',f10.4)
c note: there is x(maxvar) variables in solution 
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 subfunc(x,objf)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000)
      double precision x(maxvar),v2(maxvar)
      COMMON /cv2/v2
c===============================================================
      call func1(x,f1)
      f2=v2(1)*x(1)+v2(2)*x(2)
      objf=f1-f2
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 DCA(nvar,x0,nbundle,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,/cngrad/ngrad
      call cpu_time(time1)
      m=nvar
      nf=0
      niter=0
      ngrad=0
      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 ,w(maxdg,maxvar),prod(maxdg,maxdg),z(maxdg),fvalues(maxit)
     2 ,v2(maxvar),w1(maxdg,maxvar),w2(maxdg,maxvar)
      INTEGER ij(maxdg)
      common /csize/m,/citer/maxiter,niter,/cij/ij,jvertex,/cz/z
     1 ,/ckmin/kmin,/cndca/ndca,/cv2/v2
c====================================================================
      dist1=1.0d-07
      step0=-2.0d-01
      div=2.0d-01
      eps0=1.0d-07
      slmin=1.0d-10*slinit
      sdif=1.0d-05
      mturn=4
      ndca=0
c====================================================================
      toler=eps0
100   ndca=ndca+1
      call gradient2(x,v2)
      if (ndca.gt.1) then
       if(ndg.gt.nbundle) ndg=nbundle
       do i=1,ndg
        do j=1,m
         w2(i,j)=w1(i,j)-v2(j)
        end do
       end do
       do i=1,ndg
        do j=1,i
         prod(i,j)=0.0d+00
         do k=1,m
          prod(i,j)=prod(i,j)+w2(i,k)*w2(j,k)
         end do
         prod(j,i)=prod(i,j)
        end do
       end do
       call wolfe(ndg,prod)
       do i=1,m
        v(i)=0.0d+00
        do j=1,jvertex
         v(i)=v(i)+w2(ij(j),i)*z(j)
        END do
       END do
       r=0.0d+00
       do i=1,m
        r=r+v(i)*v(i)
       end do
       r=dsqrt(r)
       if(r.le.toler) return
      end if      
c====================================================================
      sl=slinit/div
      call subfv(x,f2)
  1   sl=div*sl
      IF(sl.lt.slmin) go to 100
      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
        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
c--------------------------------------------------------------
        do ndg=1,nbundle
            call dgrad(x,sl,g,v)
            do i=1,m
             w1(ndg,i)=v(i)
             v(i)=v(i)-v2(i)
            end do
            dotprod=0.0d+00
            do i=1,m
             dotprod=dotprod+v(i)*v(i)
            end do
            r=dsqrt(dotprod)
            IF(r.lt.eps0) GO TO 1
            IF(ndg.eq.1) then
                         rmean=r
                         kmin=1
                         rmin=r
            END if
            IF(ndg.gt.1) then
                         rmin=dmin1(rmin,r)
                         IF(r.eq.rmin) kmin=ndg
                         rmean=((ndg-1)*rmean+r)/ndg
            END if
            toler=dmax1(eps0,dist1*rmean)
            do i=1,ndg-1
             prod(ndg,i)=0.0d+00
             do j=1,m
              prod(ndg,i)=prod(ndg,i)+w(i,j)*v(j)
             end do
             prod(i,ndg)=prod(ndg,i)
            end do
            prod(ndg,ndg)=dotprod
c====================================================================
            do i=1,m
             w(ndg,i)=v(i)
            end do
            call wolfe(ndg,prod)
c================================
            do i=1,m
             v(i)=0.0d+00
             do j=1,jvertex
              v(i)=v(i)+w(ij(j),i)*z(j)
             END do
            END do
c================================
            r=0.0d+00
            do i=1,m
             r=r+v(i)*v(i)
            end do
            r=dsqrt(r)
            if(r.lt.toler) GO TO 1
c===========================================================
             do i=1,m
              g(i)=-v(i)/r
              x1(i)=x(i)+sl*g(i)
             end do
c===========================================================
             call subfv(x1,f4)
             f3=(f4-f1)/sl
             decreas=step0*r
             if(f3.lt.decreas) then
                        call armijo(x,g,f1,f5,f4,sl,step,r)
                        f2=f5
                        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 norm of 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
      kzero=0
      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,dg)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000)
      double precision x1(maxvar),g(maxvar),x(maxvar),dg(maxvar)
      common /csize/m,/cngrad/ngrad
      do k=1,m
        x1(k)=x(k)+sl*g(k)
      end do
      ngrad=ngrad+1
      call gradient1(x1,dg)
      return
      end

c===========================================================
c Line search (Armijo-type), Step 5 Algorithm 2.
c===========================================================
      subroutine armijo(x,g,f1,f5,f4,sl,step,r)
      implicit double precision (a-h,o-z)
      PARAMETER(maxvar=2000, maxdg=1000)
      common /csize/m
      double precision x(maxvar),g(maxvar),x1(maxvar)
      step=sl
      f5=f4
      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 subfv(x1,f50)
      f30=f50-f1+5.0d-02*step1*r
      IF(f30.gt.0.0d+00) RETURN
      step=step1
      f5=f50
      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
c================================================================
      call func1(x,f1)
      call func2(x,f2)
      nf=nf+1
      f=f1-f2
c=================================================================
      return
      end

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