      program test_funobj

      integer size, dim
      parameter (size = 2, dim = 1)
      
      double precision lambda
      double precision Ker(size, size*dim)
      integer Y(size)
      common /c1/ lambda 
      common /c2/ Ker
      common /c3/ Y
      
      integer total
      double precision x(2*size+dim)
      double precision F, tmp, tmpReal
      double precision G(2*size+dim)
      double precision sigma
      double precision myX(size, dim)

      sigma = 0.5
      tmp = rand(15)
      do 20, i = 1, size
         do 10 , j = 1, dim
            myX(i,j) = rand(0)
            write (*,*) "myX", myX(i,j)
 10      continue
 20   continue


      do 30 i = 1 , size
         tmpReal = rand(0)
         tmp = 4*myX(i,1)-2.0
         tmp = exp(tmp)/(1+exp(tmp)) 
	 if (tmpReal .le. tmp) then
            Y(i) = 1
         else
            Y(i) = -1
         endif
         write(*,*) "inside Y is:", Y(i)
 30   continue
      
      
      do 60 k = 1, dim
         do 50, i = 1, size
            do 40, j = 1, size
               Ker(i, (k-1)*size+j) 
     $              = exp(-(myX(i,k)-myX(j,k))**2/(2*sigma**2))               
 40         continue
 50      continue
 60   continue
      
      
      lambda = 0.5
      total = 2*size+dim
      do 70 i=1,total
         x(i) = i
         write(*,*) " x = ", x(i)
 70   continue
      
      call lbsobj(total,x,F,G)
      end



      subroutine lbsobj (n, x, f, g)

      implicit           double precision (a-h,o-z)
      double precision   x(n), g(n)
      integer            size, dim
      parameter         (zero = 0.0d+0, size = 2, dim = 1)
      double precision lambda
      double precision B(size, size*dim)
      integer Y(size)  
      common /c1/ lambda
      common /c2/ B
      common /c3/ Y
      double precision  ft(size)
      double precision sum1, sum2

      
*     -------------
*     Normal entry.
*     -------------

      sum1 = zero
      sum2 = zero
      write (*,*) "lambda = ", lambda
      
C     input the kernel matrix Ker
C      open(1,file = 'kernel.dat', status = 'old',
C     $     access = 'sequential', form = 'formatted')
      do 110, i = 1, size
C         read(1,*) (B(i,j), j=1, size*dim)
         write(*,*) "kernel is:", (B(i,j), j=1,size*dim)
 110  continue
C      close(1)
 

C     input the data y
C      open(1,file = 'y.dat', status = 'old',
C     $     access = 'sequential', form = 'formatted')
      do 120, i = 1, size
C         read(1,*) Y(i)
         write(*,*) "outside Y(i)=", Y(i) 
 120   continue
C      close(1)
 
      
C     sum1 is the sum of eta
      do 130, i = 1,size
         sum1 = sum1 + x(size+i)
 130   continue
       write(*,*) "sum of eta is", sum1

C     sum2 is the sum of theta
      do 140, j = 1, dim
         sum2 = sum2 + x(2*size+j)
 140   continue
       write(*,*) "sum of theta is", sum2

C     m is between 1 to n, ft are the components of f 
      do 170, m = 1,size
         ft(m) = zero
         do 160, j = 1, dim
            do 150, i = 1, size
               ft(m) = ft(m) + x(2*size+j)*x(i)*B(m,i+(j-1)*size)
 150         continue
 160      continue
          write(*,*) "components of f = ", m, ft(m)
 170   continue 


C     f(x) is the nonlinear objective function of nonlinear variable x
      f = zero
      do 180 , m = 1, size
         f = f - (1+y(m))*ft(m)/2 + log(1+exp(ft(m)))
 180   continue
      f = f/size + lambda * sum1 * sum2
      write(*,*) "objective f = ", f


C     g has three different parts
C     part I:   index from 1 to size
C     part II:  index from (size+1) to 2*size
C     part III: index from 2*size+1 to 2*size+d
C     part I
      do 210 i = 1, size
         g(i) = zero
         do 200 m = 1, size
            do 190 j = 1, dim
               g(i) = g(i)+ ( -(1+y(m))/2 + exp(ft(m))/(1+exp(ft(m))))
     $              *x(2*size+j)*B(m,i+(j-1)*size)/size           
 190        continue
 200     continue
         write(*,*) "grad(c):",i,g(i)
 210  continue


C     part II
      do 220 i = 1,size
         g(size+i) = lambda * sum2
         write(*,*) "grad(eta) are: ",size+i,g(size+i)
 220  continue


C     part III
      do 250 j = 1, dim
         g(2*size+j) = zero
         do 240 m = 1, size
            do 230 i = 1, size
               g(2*size+j)=g(2*size+j)+
     $              (-(1+y(m))/2+exp(ft(m))/(1+exp(ft(m))))*x(i)
     $              *B(m,i+(j-1)*size)/size             
 230        continue
 240     continue
         g(2*size+j) = g(2*size+j) + lambda*sum1
         write(*,*) "grad(theta) are: ",2*size+j,g(2*size+j)
 250  continue


      
C      open(2,file = 'c.dat', status = 'old',
C     $     access = 'sequential', form = 'formatted')  
C      do 160, i=1, n
C         write(2,*) x(i)
C 1600   continue
C      close(2)

C      open(3,file = 'f.dat', status = 'old',
C     $     access = 'sequential', form = 'formatted')
C      write(3,*) f
C      close(3)       
      
C      open(4,file = 'g.dat', status = 'old',
C     $     access = 'sequential', form = 'formatted')
C      do 170, i=1, n
C        write(4,*) g(i)
C 170   continue
C      close(4)
      
*      return

C    end of lbsobj  (funobj for lbs)
      end








