C
C                            Main driver program
C
C =============================================================================
C
      MODULE COMMON_DATA
	IMPLICIT NONE
	INTEGER,PARAMETER::DN=8
      REAL*8,DIMENSION(DN)::T1,T2
	INTEGER,DIMENSION(DN)::D1,D2,TD1,TD2
	INTEGER::N1,N2
	REAL*8,ALLOCATABLE,DIMENSION(:)::P1,P2,SPP1,SPP2
	INTEGER::NEW_XT
	END MODULE COMMON_DATA

      SUBROUTINE GENGRID
      USE COMMON_DATA
	IMPLICIT NONE
	INTEGER::I,J,K,L,IERROR
	REAL*8::TEMP(DN),TMIN
	REAL*8,ALLOCATABLE,DIMENSION(:)::DD,RR,SS
	INTEGER,ALLOCATABLE,DIMENSION(:)::KA,LA
	INTEGER,ALLOCATABLE,DIMENSION(:,:)::CN11,CN10,CN01,CN00
	REAL*8,PARAMETER::EPS=1.0E-08
	OPEN(UNIT=30,FILE='BIV.txt',STATUS='OLD',IOSTAT=IERROR)
	DO I=1,DN,1
       READ(30,*,IOSTAT=IERROR)T1(I),T2(I),D1(I),D2(I)
      END DO
	CLOSE(UNIT=30)
      TMIN=MINVAL(T1,D1==1)-1.0;TEMP=TMIN;TEMP(DN)=MAXVAL(T1)
	DO I=DN-1,1,-1
       IF (COUNT(D1==1 .AND. T1<TEMP(I+1))>0) THEN
	  TEMP(I)=MAXVAL(T1,D1==1 .AND. T1<TEMP(I+1))
	 ELSE
	  TEMP(I)=TMIN    
	 END IF
	END DO
	N1=COUNT(TEMP>TMIN)
	ALLOCATE(SPP1(N1),P1(N1),DD(N1),RR(N1),SS(N1))
	SPP1=TEMP(DN-N1+1:DN)
	DO I=1,N1,1
       DD(I)=REAL(COUNT(DABS(T1-SPP1(I))<EPS .AND. D1==1),8)
	 RR(I)=DD(I)+REAL(COUNT(T1>SPP1(I)),8)
	END DO
	P1(1)=DD(1)/RR(1);SS(1)=1.0-P1(1)
	DO I=2,N1-1,1
       P1(I)=SS(I-1)*DD(I)/RR(I)
	 SS(I)=SS(I-1)-P1(I)
	END DO
	P1(N1)=1.0-SUM(P1(1:(N1-1)))
	DEALLOCATE(DD,RR,SS)

      TMIN=MINVAL(T2,D2==1)-1.0;TEMP=TMIN;TEMP(DN)=MAXVAL(T2)
	DO I=DN-1,1,-1
       IF (COUNT(D2==1 .AND. T2<TEMP(I+1))>0) THEN
	  TEMP(I)=MAXVAL(T2,D2==1 .AND. T2<TEMP(I+1))
	 ELSE
	  TEMP(I)=TMIN    
	 END IF
	END DO
	N2=COUNT(TEMP>TMIN)
	ALLOCATE(SPP2(N2),P2(N2),DD(N2),RR(N2),SS(N2))
	SPP2=TEMP(DN-N2+1:DN)
	DO I=1,N2,1
       DD(I)=REAL(COUNT(DABS(T2-SPP2(I))<EPS .AND. D2==1),8)
	 RR(I)=DD(I)+REAL(COUNT(T2>SPP2(I)),8)
	END DO
	P2(1)=DD(1)/RR(1);SS(1)=1.0-P2(1)
	DO I=2,N2-1,1
       P2(I)=SS(I-1)*DD(I)/RR(I)
	 SS(I)=SS(I-1)-P2(I)
	END DO
	P2(N2)=1.0-SUM(P2(1:(N2-1)))
	DEALLOCATE(DD,RR,SS)

      ALLOCATE(KA(N1),LA(N2))
	DO I=1,N1,1
	 KA(I)=I
	END DO
	DO I=1,N2,1
       LA(I)=I
	END DO

	DO I=1,DN,1
       IF (D1(I)==1) THEN
	  TD1(I)=MINVAL(KA,DABS(SPP1-T1(I))<EPS)
	 ELSE
	  TD1(I)=MINVAL(KA,SPP1>=T1(I))
	 END IF

	 IF (D2(I)==1) THEN
	  TD2(I)=MINVAL(LA,DABS(SPP2-T2(I))<EPS)
	 ELSE
	  TD2(I)=MINVAL(LA,SPP2>=T2(I))
	 END IF
	END DO
	DEALLOCATE(KA,LA)

	ALLOCATE(CN11(N1,N2),CN10(N1,N2),CN01(N1,N2),CN00(N1,N2))
	DO K=1,N1,1
       DO L=1,N2,1
        CN11(K,L)=COUNT(TD1==K .AND. TD2==L .AND. D1==1 .AND. D2==1)
	  CN10(K,L)=COUNT(TD1==K .AND. TD2==L .AND. D1==1 .AND. D2==0)
	  CN01(K,L)=COUNT(TD1==K .AND. TD2==L .AND. D1==0 .AND. D2==1)
	  CN00(K,L)=COUNT(TD1==K .AND. TD2==L .AND. D1==0 .AND. D2==0)
	 END DO
	END DO
	WRITE(*,*)'%%%%%%%%%%%%%%%%%%%%%%%%'
	DO L=1,N2,1
	WRITE(*,*)N2-L+1,CN11(:,N2-L+1)
	END DO
      DO L=1,N2,1
      WRITE(*,*)N2-L+1,CN10(:,N2-L+1)
	END DO
      DO L=1,N2,1
	WRITE(*,*)N2-L+1,CN01(:,N2-L+1)
	END DO
	DO L=1,N2,1
      WRITE(*,*)N2-L+1,CN00(:,N2-L+1)
	END DO
	DEALLOCATE(CN11,CN10,CN01,CN00)
	END SUBROUTINE GENGRID



      program example
	USE COMMON_DATA
C
      implicit none
C
C     include the Ipopt return codes
C
      include 'IpReturnCodes.inc'
C
C     Size of the problem (number of variables and equality constraints)
C
      INTEGER NMAX,MMAX,IDX_STY
      parameter  (NMAX=DN*DN,MMAX=2*DN,IDX_STY=1)
	integer     N,     M,     NELE_JAC,     NELE_HESS
C
C     Space for multipliers and constraints
C
      double precision LAM(MMAX),G(MMAX)
C
C     Vector of variables
C
      double precision X(NMAX)
C
C     Vector of lower and upper bounds
C
      double precision X_L(NMAX), X_U(NMAX), Z_L(NMAX), Z_U(NMAX)
      double precision G_L(MMAX), G_U(MMAX)
C
C     Private data for evaluation routines
C     This could be used to pass double precision and integer arrays untouched
C     to the evaluation subroutines EVAL_*
C
      double precision DAT(MMAX)
      integer IDAT(1)
C
C     Place for storing the Ipopt Problem Handle
C
C     for 32 bit platforms
      integer IPROBLEM
      integer IPCREATE
CC     for 64 bit platforms:
C      integer*8 IPROBLEM
C      integer*8 IPCREATE
C
      integer IERR
      integer IPSOLVE, IPADDSTROPTION
      integer IPADDNUMOPTION, IPADDINTOPTION
      integer IPOPENOUTPUTFILE
C
      double precision F
      integer I
C
C     The following are the Fortran routines for computing the model
C     functions and their derivatives - their code can be found furhter
C     down in this file.
C
      external EV_F, EV_G, EV_GRAD_F, EV_JAC_G, EV_HESS

      CALL GENGRID
	N=N1*N2
      M=N1+N2-1
	NELE_JAC=N1*(2*N2-1)
	NELE_HESS=N*(N+1)/2
C
C     Set initial point and bounds:
C
      X=0.0
	DO I=1,N1,1
       DO J=1,N2,1
        X((I-1)*N2+J)=P1(I)*P2(J)
	 END DO
	END DO
      X_L=0.0;X_U=X_L
	X_U(1:N)=1.0
C
C     Set bounds for the constraints
C
	G_L=0.0;G_U=G_L
C
C     First create a handle for the Ipopt problem (and read the options
C     file)
C
      IPROBLEM = IPCREATE(N, X_L, X_U, M, G_L, G_U, NELE_JAC, NELE_HESS,
     1     IDX_STY, EV_F, EV_G, EV_GRAD_F, EV_JAC_G, EV_HESS)
      if (IPROBLEM.eq.0) then
         write(*,*) 'Error creating an Ipopt Problem handle.'
         stop
      endif
C
C     Open an output file
C
      IERR = IPOPENOUTPUTFILE(IPROBLEM, 'IPOPT.OUT', 5)
      if (IERR.ne.0 ) then
         write(*,*) 'Error opening the Ipopt output file.'
         goto 9000
      endif
C
C     Note: The following options are only examples, they might not be
C           suitable for your optimization problem.
C
C     Set a string option
C
      IERR = IPADDSTROPTION(IPROBLEM, 'mu_strategy', 'adaptive')
      if (IERR.ne.0 ) goto 9990
C
C     Set an integer option
C
      IERR = IPADDINTOPTION(IPROBLEM, 'max_iter', 1000)
      if (IERR.ne.0 ) goto 9990
C
C     Set a double precision option
C
      IERR = IPADDNUMOPTION(IPROBLEM, 'tol', 1.d-7)
      if (IERR.ne.0 ) goto 9990
C
C     As a simple example, we pass the constants in the constraints to
C     the EVAL_C routine via the "private" DAT array.
C
      DAT = 0.0
	DO I=1,N1,1
       DAT(I)=-P1(I)
	END DO
	DO I=1,N2-1,1
       DAT(N1+I)=-P2(I)
	END DO
	IDAT(1)=1

      IERR = IPSOLVE(IPROBLEM, X, G, F, LAM, Z_L, Z_U, IDAT, DAT)
C     
C     Output:
C
      if( IERR.eq.IP_SOLVE_SUCCEEDED ) then
         write(*,*)
         write(*,*) 'The solution was found.'
         write(*,*)
         write(*,*) 'The final value of the objective function is ',F
         write(*,*)
         write(*,*) 'The optimal values of X are:'
         write(*,*)
         do i = 1, N
            write(*,*) 'X  (',i,') = ',X(i)
         enddo
         write(*,*)
         write(*,*) 'The multipliers for the lower bounds are:'
         write(*,*)
         do i = 1, N
            write(*,*) 'Z_L(',i,') = ',Z_L(i)
         enddo
         write(*,*)
         write(*,*) 'The multipliers for the upper bounds are:'
         write(*,*)
         do i = 1, N
            write(*,*) 'Z_U(',i,') = ',Z_U(i)
         enddo
         write(*,*)
         write(*,*) 'The multipliers for the equality constraints are:'
         write(*,*)
         do i = 1, M
            write(*,*) 'LAM(',i,') = ',LAM(i)
         enddo
         write(*,*)
      else
         write(*,*)
         write(*,*) 'An error occoured.'
         write(*,*) 'The error code is ',IERR
         write(*,*)
      endif
C
 9000 continue
C
C     Clean up
C
	DEALLOCATE(P1,P2,SPP1,SPP2)
      call IPFREE(IPROBLEM)
      stop
C
 9990 continue
      write(*,*) 'Error setting an option'
      goto 9000

      end
C
C =============================================================================
C
C                    Computation of objective function
C
C =============================================================================
C
      subroutine EV_F(N, X, NEW_X, F, IDAT, DAT, IERR)
	USE COMMON_DATA
      implicit none
      integer N, NEW_X
      double precision F, X(N)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR
	INTEGER::I,J,K,L
	REAL*8::QQ(N1,N2)
	DO K=1,N1,1
       DO L=1,N2,1
        QQ(K,L)=X((K-1)*N2+L)
	 END DO
	END DO
	F=1.0
	DO I=1,DN,1
       IF (D1(I)==1 .AND. D2(I)==1) THEN
	  F=F-DLOG(QQ(TD1(I),TD2(I)))
	 ELSE IF (D1(I)==1 .AND. D2(I)==0) THEN
	  F=F-DLOG(SUM(QQ(TD1(I),TD2(I):N2)))
	 ELSE IF (D1(I)==0 .AND. D2(I)==1) THEN
	  F=F-DLOG(SUM(QQ(TD1(I):N1,TD2(I))))
	 ELSE
	  F=F-DLOG(SUM(QQ(TD1(I):N1,TD2(I):N2)))
	 END IF
	END DO
C	WRITE(*,*)'F=',F
      IERR = 0
      return
      end
C
C =============================================================================
C
C                Computation of gradient of objective function
C
C =============================================================================
C
      subroutine EV_GRAD_F(N, X, NEW_X, GRAD, IDAT, DAT, IERR)
	USE COMMON_DATA
      implicit none
      integer N, NEW_X
      double precision GRAD(N), X(N)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR
	INTEGER::I,J,K,L
	REAL*8::QQ(N1,N2)
	DO K=1,N1,1
       DO L=1,N2,1
        QQ(K,L)=X((K-1)*N2+L)
	 END DO
	END DO
	GRAD=0.0
      DO K=1,N1,1
       DO L=1,N2,1
	  DO I=1,DN,1
	   IF (TD1(I)==K.AND.TD2(I)==L.AND.D1(I)==1.AND.D2(I)==1) THEN
	    GRAD((K-1)*N2+L)=GRAD((K-1)*N2+L)-1.0/QQ(K,L)
	   ELSE IF(TD1(I)==K.AND.TD2(I)<=L.AND.D1(I)==1.AND.D2(I)==0) THEN
	    GRAD((K-1)*N2+L)=GRAD((K-1)*N2+L)-1.0/SUM(QQ(K,TD2(I):N2))
	   ELSE IF(TD1(I)<=K.AND.TD2(I)==L.AND.D1(I)==0.AND.D2(I)==1) THEN
	    GRAD((K-1)*N2+L)=GRAD((K-1)*N2+L)-1.0/SUM(QQ(TD1(I):N1,L))
	   ELSE IF(TD1(I)<=K.AND.TD2(I)<=L.AND.D1(I)==0.AND.D2(I)==0) THEN
	    GRAD((K-1)*N2+L)=GRAD((K-1)*N2+L)-
	1 	     1.0/SUM(QQ(TD1(I):N1,TD2(I):N2))  
	   END IF
	  END DO
	 END DO
	END DO
      IERR = 0
      return
      end
C
C =============================================================================
C
C                     Computation of equality constraints
C
C =============================================================================
C
      subroutine EV_G(N, X, NEW_X, M, G, IDAT, DAT, IERR)
	USE COMMON_DATA
      implicit none
      integer N, NEW_X, M
      double precision G(M), X(N)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR,I,J
	G=0.0
	DO I=1,N1,1
	 G(I)=DAT(I)
       DO J=1,N2,1
        G(I)=G(I)+X((I-1)*N2+J)
	 END DO
	END DO
      
	DO J=1,N2-1,1
	 G(N1+J)=DAT(N1+J)
       DO I=1,N1,1
        G(N1+J)=G(N1+J)+X((I-1)*N2+J)
	 END DO
	END DO
C	WRITE(*,*)'G=',G
C      G = SUM(X)
      IERR = 0
      return
      end
C
C =============================================================================
C
C                Computation of Jacobian of equality constraints
C
C =============================================================================
C
      subroutine EV_JAC_G(TASK, N, X, NEW_X, M, NZ, ACON, AVAR, A,
     1     IDAT, DAT, IERR)
	USE COMMON_DATA
      integer TASK, N, NEW_X, M, NZ
      double precision X(N), A(NZ)
      integer ACON(NZ), AVAR(NZ)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR, I,J
C
C     structure of Jacobian:
C
      integer AVAR1(N1*(2*N2-1)), ACON1(N1*(2*N2-1))
	DO I=1,N1,1
	 DO J=1,N2,1
        AVAR1((I-1)*N2+J)=(I-1)*N2+J
        ACON1((I-1)*N2+J)=I
	 END DO
	END DO

	DO J=1,N2-1,1
       DO I=1,N1,1
	  AVAR1(N1*N2+(I-1)*N2+J)=(I-1)*N2+J
	  ACON1(N1*N2+(I-1)*N2+J)=N1+J
	 END DO
	END DO
C
      if( TASK.eq.0 ) then
        do I = 1, N1*(2*N2-1)
          AVAR(I) = AVAR1(I)
          ACON(I) = ACON1(I)
        enddo
      else
        A=1.0
      endif
      IERR = 0
      return
      end
C
C =============================================================================
C
C                Computation of Hessian of Lagrangian
C
C =============================================================================
C
      subroutine EV_HESS(TASK, N, X, NEW_X, OBJFACT, M, LAM, NEW_LAM,
     1     NNZH, IRNH, ICNH, HESS, IDAT, DAT, IERR)
	USE COMMON_DATA
      implicit none
      integer TASK, N, NEW_X, M, NEW_LAM, NNZH
      double precision X(N), OBJFACT, LAM(M), HESS(NNZH)
      integer IRNH(NNZH), ICNH(NNZH)
      double precision DAT(*)
      integer IDAT(*)
      integer IERR
C
C     structure of Hessian:
C
      integer IRNH1(N*(N+1)/2), ICNH1(N*(N+1)/2)
      INTEGER::I,J,K,L,K1,L1,W,V,NN2
	REAL*8::QQ(N1,N2)
	
	NN2=N*(N+1)/2
	DO I=1,N,1
       DO J=1,I,1
        IRNH1(I*(I-1)/2+J)=I
	  ICNH1(I*(I-1)/2+J)=J
	 END DO
	END DO


      if( TASK.eq.0 ) then
         do I = 1, NN2,1
            IRNH(I) = IRNH1(I)
            ICNH(I) = ICNH1(I)
         enddo
      else
	 DO K=1,N1,1
        DO L=1,N2,1
         QQ(K,L)=X((K-1)*N2+L)
	  END DO
	 END DO
	 HESS=0.0
C    
C     CASE(1) K1=K AND L1=L
C
      DO K=1,N1,1
       DO L=1,N2,1
        K1=K;L1=L
        W=(K-1)*N2+L
	  V=(K1-1)*N2+L1
	  DO I=1,DN,1
	   IF (TD1(I)==K.AND.TD2(W)==L.AND.D1(I)==1.AND.D2(I)==1) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+1.0/QQ(K,L)**2
	   ELSE IF(TD1(I)==K.AND.TD2(I)<=L.AND.D1(I)==1.AND.D2(I)==0) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
	1	      1.0/SUM(QQ(K,TD2(I):N2))**2
	   ELSE IF(TD1(I)<=K.AND.TD2(I)==L.AND.D1(I)==0.AND.D2(I)==1) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
     1		  1.0/SUM(QQ(TD1(I):N1,L))**2
	   ELSE IF(TD1(I)<=K.AND.TD2(I)<=L.AND.D1(I)==0.AND.D2(I)==0) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
	1 	     1.0/SUM(QQ(TD1(I):N1,TD2(I):N2))**2  
	   END IF
	  END DO
	 END DO
	END DO
C    
C     CASE(2) K1=K AND L1<L
C
      DO K=1,N1,1
	 DO L=2,N2,1
        DO L1=1,L-1,1
         K1=K
         W=(K-1)*N2+L
	   V=(K1-1)*N2+L1
	   DO I=1,DN,1
	   
	  IF(TD1(I)==K.AND.TD2(I)<=L1.AND.D1(I)==1.AND.D2(I)==0) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
	1	      1.0/SUM(QQ(K,TD2(I):N2))**2
	  ELSE IF(TD1(I)<=K.AND.TD2(I)<=L1.AND.D1(I)==0.AND.D2(I)==0) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
	1 	     1.0/SUM(QQ(TD1(I):N1,TD2(I):N2))**2  
	  END IF

         END DO
	  END DO
	 END DO
	END DO 

C    
C     CASE(3) K1<K AND L1=L
C
      DO K=2,N1,1
	 DO L=1,N2,1
        DO K1=1,K-1,1
         L1=L
         W=(K-1)*N2+L
	   V=(K1-1)*N2+L1
	   DO I=1,DN,1
	   
	  IF(TD1(I)<=K1.AND.TD2(I)==L.AND.D1(I)==0.AND.D2(I)==1) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
     1		  1.0/SUM(QQ(TD1(I):N1,L))**2
	  ELSE IF(TD1(I)<=K1.AND.TD2(I)<=L.AND.D1(I)==0.AND.D2(I)==0) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
	1 	     1.0/SUM(QQ(TD1(I):N1,TD2(I):N2))**2  
	  END IF
	   
	   END DO
	  END DO
	 END DO
	END DO

C
C     CASE(4.1) K1<K L1<L
C
      DO K=2,N1,1
	 DO L=2,N2,1
	  DO K1=1,K-1,1
	   DO L1=1,L-1,1
          W=(K-1)*N2+L
	    V=(K1-1)*N2+L1
	    DO I=1,DN,1

	  IF(TD1(I)<=K1.AND.TD2(I)<=L1.AND.D1(I)==0.AND.D2(I)==0) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
	1 	     1.0/SUM(QQ(TD1(I):N1,TD2(I):N2))**2  
	  END IF

	    END DO
         END DO
        END DO
	 END DO
	END DO	 

C
C     CASE(4.2) K1<K L1>L
C
      DO K=2,N1,1
	 DO L=1,N2-1,1
	  DO K1=1,K-1,1
	   DO L1=L+1,N2,1
          W=(K-1)*N2+L
	    V=(K1-1)*N2+L1
	    DO I=1,DN,1

	  IF(TD1(I)<=K1.AND.TD2(I)<=L.AND.D1(I)==0.AND.D2(I)==0) THEN
	     HESS(W*(W-1)/2+V)=HESS(W*(W-1)/2+V)+
	1 	     1.0/SUM(QQ(TD1(I):N1,TD2(I):N2))**2  
	  END IF

	    END DO
         END DO
        END DO
	 END DO
	END DO	 
      endif
      IERR = 0
      return
      end