c***********************************************************************
c     Main program with nonlinear test
c***********************************************************************
      implicit none
      external objfun
      character*80 strdat
      integer*4 nolin,iestat
      
      strdat= 'ntest.dat'
      nolin=1 
      iestat=0
      call pprn1(nolin,strdat,iestat,objfun)
      write(*,*) 'exit=',iestat
      end 
      
C***********************************************************************
C     Objective function
C***********************************************************************
      subroutine objfun (nar,nat,x,f,g)
      implicit none
      integer*4 nar,nat
      real*8 f,x(0:nar,nat),g(0:nar,nat)
C     
c     constants per divisions
      real*8 coef1,coef2,coef3,coef23
      parameter(coef1= 1000.0d0)
      parameter(coef2= 1000.0d0)
      parameter(coef3= 1200.0d0)
      parameter(coef23= coef2*coef3)
c     variables locals
      integer*4 j,k
      real*8 xx,xxn,f1,f2,f3,xx2,xnx,u_u,r,r2,r3,r4,sqr

      f= 0.0d0
      do k= 1,nat
         do j= 0,nar
            g(j,k)= 0.0d0
         end do
      end do
c
      do k= 1,nat
c     inicialitzacions per cada article
         f1= 0.0d0
         f2= 0.0d0
         f3= 0.0d0
         u_u= 1.0d0
         xx= x(1,k)
c     bucle per article k (arc 1...n-1)
         do j= 1,nar-1
            xxn= x(j+1,k)
            xx2= xx*xx
            xnx= xx-xxn
            u_u= -u_u
            
            f1= f1+xx2
            sqr= dsqrt(1.0d0+xx2+xnx*xnx)
            f2= f2+sqr
            f3= f3+u_u*xx
            
            g(j,k)= g(j,k)+2.0d0*xx/coef1     !part f1
            r= 2.0d0*sqr*coef2
            g(j,k)= g(j,k)+(2.0d0*(xx+xnx))/r !part f2
            g(j+1,k)= g(j+1,k)-(2.0d0*xnx)/r  !part f2
            
            xx= xxn
         end do
c     tractament ultim arc
         u_u= -u_u
         f1= f1+xx*xx
         f3= f3+u_u*xx

         g(nar,k)= g(nar,k)+2.0d0*xx/coef1   !part f1         
c     valor funobj
         r= (f3+10.0d0)
         r2= r*r
         r4= r2*r2
         f= f + f1/coef1 + f2/coef2 + r4/coef23
c     gradient part f3
         r3= r2*r
         u_u= 1.0d0
         do j= 1,nar
            u_u= -u_u
            g(j,k)= g(j,k)+4.0d0*r3*u_u/coef23
         end do
      end do
      
      END

