      program test_funobj
      
      integer size, dim, basis
      include 'param.dat'
      
      double precision myX(size,dim)
      double precision linT(size,1+dim+dim*(dim-1)/2)
      double precision baseX(basis,dim)
      common /c1/ myX 
      common /c2/ linT
      common /c3/ baseX
      
      double precision trProb(size)
      integer Y(size)
      
      double precision x(1 + dim + dim*(dim-1)/2 + basis*dim +
     $     3*basis*dim*(dim-1)/2)
      double precision F
      double precision G(1 + dim + dim*(dim-1)/2 + basis*dim +
     $     3*basis*dim*(dim-1)/2)
      

C     generate original myX, Y array with true prob
      call simuXY(trProb,Y)

      n = 1+dim + dim*(dim-1)/2 + basis*dim + 3*basis*dim*(dim-1)/2

      l1 = 1+dim+dim*(dim-1)/2
      l2 = dim*basis
      l3 = basis*dim*(dim-1)
      l4 = basis*dim*(dim-1)/2

      do 80 j=1,l1
         x(j) = 4
         write(*,*) " x = ", x(j)
 80   continue
      do 90 j= 1+l1,l1+l2
         x(j) = 1
         write(*,*) " x = ", x(j)
 90   continue
      do 100 j= 1+l1+l2,l1+l2+l3
         x(j) = 2
         write(*,*) " x = ", x(j)
 100  continue
      do 110 j= 1+l1+l2+l3,l1+l2+l3+l4
         x(j) = 3
         write(*,*) " x = ", x(j)
 110  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 myX(size, dim)
      double precision linT(size, 1+dim+dim*(dim-1)/2)
      double precision baseX(basis, dim)
      common /c1/ myX
      common /c2/ linT
      common /c3/ baseX
C     must claim the function type in the main program
      double precision genKer

      integer leap, index, ltd
      double precision ft(size)
      double precision expf(size)

*     -------------
*     Normal entry.
*     -------------
C     print mode and n
      
      write(*,*)  '***********n=********',n
      
      
C     print x
      do 120, j = 1,n
         write(*,*) " x is :",j, x(j)
 120  continue
      

      ltd = 1 + dim + dim*(dim-1)/2

      write(*,*) "linT is"
      do 130, i = 1, size
         write(*,*) (linT(i,j), j = 1, ltd)
 130  continue
      
      
      do 230, i = 1,size
         ft(i) = zero
         tmp = zero
         do 140, j = 1,ltd
            ft(i) = ft(i) + linT(i,j)*x(j)
            tmp = tmp + linT(i,j)*x(j)
 140     continue
         write(*,*) "for i =", i,"     Td =", tmp
C     compute K1*c1 
         tmp = zero
         index = ltd
         do 160, j = 1,dim
            do 150, k = 1,basis
               ft(i) = ft(i) + genKer(myX(i,j),baseX(k,j))*
     $              x(index + (j-1)*basis+k)
               tmp = tmp + genKer(myX(i,j),baseX(k,j))*
     $              x(index + (j-1)*basis+k) 
 150        continue
 160     continue
         write(*,*) "for i=", i, "     K1*c1 = ", tmp
C     compute K2*c2
         tmp = zero
         leap = 1
         index = index + dim*basis
         do 190, l = 1,dim
            do 180, j = 1,dim
               if (j .ne. l) then
                  do 170, k = 1, basis
                     ft(i) = ft(i) + genKer(myX(i,l),baseX(k,l))*
     $                    (myX(i,j)-0.5)*(baseX(k,j)-0.5)*
     $                    x(index + leap)
                     tmp = tmp + genKer(myX(i,l),baseX(k,l))*
     $                    (myX(i,j)-0.5)*(baseX(k,j)-0.5)*
     $                    x(index + leap)
                     leap = leap + 1
 170              continue
               end if
 180        continue
 190     continue
         write(*,*) "for i=", i, "     K2*c2 = ", tmp
C     compute K3*c3
         tmp = zero
         leap = 1
         index = index + basis*dim*(dim-1)
         do 220, l = 1,dim-1
            do 210, j = l+1,dim
               do 200, k = 1,basis
                  ft(i) = ft(i) + genKer(myX(i,l),baseX(k,l))*
     $                 genKer(myX(i,j),baseX(k,j))*x(index+leap)
                  tmp = tmp + genKer(myX(i,l),baseX(k,l))*
     $                 genKer(myX(i,j),baseX(k,j))*x(index+leap)
                  leap = leap + 1
 200           continue
 210        continue
 220     continue
         write(*,*) "for i=", i, "    K3*c3 = ", tmp
         write(*,*) "components of f: f_", i, ft(i)
 230  continue
         
      
      
C     f(x) is the nonlinear objective function of nonlinear variable x
      f = zero
      do 240, i = 1, size
         f = f + log(1+exp(ft(i)))/size
         expf(i) = exp(ft(i))/(1+exp(ft(i)))
 240  continue
      write(*,*) "objective F(x) = ", f
      
         

C     g has index j from 1 to n = ltd + basis*dim + 3*dim*(dim-1)*/2       
      do 260, j = 1, ltd
         g(j) = zero
         do 250, i = 1, size 
            g(j) = g(j)+ expf(i)*linT(i,j)/size         
 250     continue
 260  continue

C     compute g(c1)= dF/dc1
      index = ltd
      do 290, j = 1, dim
         do 280, k = 1, basis 
            tmp = zero
            do 270 i = 1, size  
               tmp = tmp + expf(i)* genKer(myX(i,j),baseX(k,j))         
 270     continue
         g(index+(j-1)*basis+k) = tmp/size
 280  continue
 290  continue


C     compute g(c2) = dF/d(c2)
      leap = 1
      index = index + dim*basis
      do 330, l = 1,dim
         do 320, j = 1,dim
            if (j .ne. l) then
               do 310, k = 1, basis
                  tmp = zero
                  do 300, i = 1, size
                     tmp = tmp + expf(i)* genKer(myX(i,l),baseX(k,l))*
     $                       (myX(i,j)-0.5)*(baseX(k,j)-0.5)
 300              continue
                  g(index + leap) = tmp/size
                  leap = leap + 1
 310           continue
            end if
 320     continue
 330  continue


C     compute g(c3) = dF/d(c3)
      leap = 1
      index = index + basis*dim*(dim-1)
      do 370, l = 1,dim-1
         do 360, j = l+1,dim
            do 350, k = 1,basis
               tmp = zero
               do 340, i = 1,size
                  tmp = tmp + expf(i)*genKer(myX(i,l),baseX(k,l))*
     $                 genKer(myX(i,j),baseX(k,j))
 340           continue
               g(index + leap) = tmp/size
               leap = leap + 1               
 350        continue
 360     continue
 370  continue
            
      do 380 j = 1,ltd+basis*dim+3*basis*dim*(dim-1)/2
         write(*,*) "grad(d):",j,g(j)
 380  continue  
     

      end
C     end of lbsobj  (funobj for lbs)
      
  
