      program test_funobj
      
      integer size, dim, basis
      include 'param.dat'
      
      double precision Ker(size,dim*basis)
      double precision linT(size,dim+1)
      common /c1/ Ker 
      common /c2/ linT
      
      double precision myX(size, dim)
      double precision trProb(size)
      integer Y(size)
      
      
      double precision x(1+dim+basis*dim)
      double precision F
      double precision G(1+dim+basis*dim)
      
C     generate original myX, Y array with true prob
      call simuXY(myX,trProb,Y)

      n = 1+dim+basis*dim
      do 80 j=1,n
         x(j) = 1
         write(*,*) " x = ", x(j)
 80   continue
      
      call lbsobj(n,x,F,G)
      end
      
      
      subroutine lbsobj ( n, x, f, g )
C     subroutine lbsobj ( mode, n, x, f, g, nstate, nprob, z, nwcore )
      
      implicit           double precision (a-h,o-z)
      double precision   x(n), g(n)
      integer            size, dim, basis
      include 'param.dat'
      double precision B(size, dim*basis)
      double precision linT(size, dim+1)
      common /c1/ B
      common /c2/ linT
      double precision  ft(size)
      

*     -------------
*     Normal entry.
*     -------------
C     print mode and n
      mm = dim*basis
      
      write(*,*)  '***********n=********',n
      
      
C     print x
      do 120, j = 1,n
         write(*,*) " x is :", n, j, x(j)
 120  continue
      
      
      
C     ft are the components of f 
C     ft = T*x[1:dim+1]+B*x[dim+2:1+dim+mm]
      do 160, i = 1,size
         ft(i) = zero
         do 140, j = 1,dim+1
            ft(i) = ft(i) + linT(i,j)*x(j)
 140     continue
         do 150, j = 1,mm
            ft(i) = ft(i) + B(i,j)*x(1+dim+j)
 150     continue
         write(*,*) "components of f: f_", i, ft(i)
 160  continue
      

      
C     f(x) is the nonlinear objective function of nonlinear variable x
      f = zero
      do 180, i = 1, size
         f = f + log(1+exp(ft(i)))/size
 180  continue
      write(*,*) "objective F(x) = ", f
      
      
      
C     g has index j from 1 to n = 1+dim+mm     
      do 210 j = 1, dim+1
         g(j) = zero
         do 200 i = 1, size 
            g(j) = g(j)+ exp(ft(i))/(1+exp(ft(i)))*linT(i,j)/size         
 200     continue
 210  continue
      
      
      do 230 j = 1, mm
         g(1+dim+j) = zero
         do 220 i = 1, size 
            g(1+dim+j) = g(1+dim+j) + 
     $           exp(ft(i))/(1+exp(ft(i)))*B(i,j)/size         
 220     continue
 230  continue
         
      
      do 240 j = 1,dim+1
         write(*,*) "grad(d):",j,g(j)
 240  continue
      
      
      do 250 j = 1,mm
         write(*,*) "grad(c):",j,g(j+dim+1)
 250  continue
      end
      
C     end of lbsobj  (funobj for lbs)
      
      
      

