
c     -*-fortran-*-
c==>projgrad.f
      program pgrad

c     A simple projected gradient algorithm written to demonstrate
c     the GAMS Callable Program Library (CPLIB)

c     this file is the main program of the solver.  It is called
c     from the command line or batch file, and just passes control
c     to cpmain.  Since we do fortran I/O, the initial call
c     to the solver system must be to a fortran program.

      call cpmain
      end


      subroutine corerq
      integer n, nnz, nwucor

      call cpgeti ('N',n)             ! size of problem
      call cpgeti ('NADIM',nnz)       ! nonzeros in Jacobian

      nwucor = 4 * n + 0 * nnz        ! we won't use Jacobian
      call cpputi ('ISTYPE',2)        ! solves general MCP's
      call cpputi ('NWUCOR', nwucor)  ! solvers memory requirement
      return
      end


      subroutine solver (work, nwucor)
      integer nwucor
      double precision work(nwucor)
      common /samplcom/ iolog, iosta, iterlim
      integer iolog, iosta, iterlim
      integer n

      call cpgeti ('IOLOG',iolog)     ! log file unit #
      call cpgeti ('IOSTA',iosta)     ! status file unit #
      call cpgeti ('ITERLIM',iterlim) ! iteration limit
      call cpgeti ('N',n)

      call cpputi ('STARTC', 0)
      write (iolog,1000)
      write (iosta,1000)
      call cpputi ('STOPC', 0)

      call projgrad (n, work(1), work(n+1), work(2*n+1), work(3*n+1))
      return
 1000 format ('Sample solver programmed by Steve Dirkse')
      end


      subroutine projgrad (n, lower, upper, z, F)
      integer n
      double precision lower(n), upper(n), z(n), F(n)
      common /samplcom/ iolog, iosta, iterlim
      integer iolog, iosta, iterlim
      integer iteration, i
      double precision metric, obj, stepsize, dt, nonsmoothnorm

      stepsize = 0.5                  ! should be set by user, but . .
      iteration = 0

      call cpbnds (z, lower, upper, n)
      call cpfunf (z, F, n)
      write (iolog, 1000)
      write (iolog, 1001)
      metric = nonsmoothnorm (n, z, F, lower, upper)
      call cpgetd ('OBJ', obj)
      write (iolog, 1010) iteration, metric, obj
      do while (iteration .lt. iterlim)
         if (metric .lt. 1e-6) then     ! convergence
            call cpsoln (z, n)
            call cpputi ('MODSTA', 1)   ! model solved
            call cpputi ('SOLSTA', 1)   ! solution computed
            return
         endif
         do 100 i = 1, n
            dt = z(i) - stepsize * F(i)
            dt = max (dt, lower(i))
            z(i) = min (dt, upper(i))
 100     continue
         iteration = iteration + 1
         call cpfunf (z, F, n)
         metric = nonsmoothnorm (n, z, F, lower, upper)
         call cpgetd ('OBJ', obj)
         write (iolog, 1010) iteration, metric, obj
      enddo
      call cpputi ('MODSTA', 7)      ! intermediate nonoptimal
      call cpputi ('SOLSTA', 2)      ! iteration interrupt
      return
 1000 format (/,' iterate', t19, 'residual norm', t38, 'CPLIB norm')
 1001 format (  ' -------', t19, '-------------', t38, '----------')
 1010 format (i6,t17,f15.7,t33,f15.7)
      end


      double precision function nonsmoothnorm (n, z, F, lower, upper)
      integer n
      double precision  z(n), F(n), lower(n), upper(n)
      integer i
      double precision dt, norm, zero

      zero = 0
      norm = 0
      do 100 i = 1, n
         if (z(i) .le. lower(i)) then
            dt = min (zero, F(i))
            norm = norm + dt**2
         elseif (z(i) .ge. upper(i)) then
            dt = max (zero, F(i))
            norm = norm + dt**2
         else
            norm = norm + f(i)**2
         endif
 100  continue
      nonsmoothnorm = sqrt(norm)
      return
      end
