GEsystem.f

C     =================================================================

C     File: GEsystem.f

C     =================================================================

C     =================================================================

C     Module: Subroutines that define the problem

C     =================================================================

C     Last update of any of the component of this module:

C     October 4, 2008.

C     Users are encouraged to download periodically updated versions of

C     this code at the Software home page:

C

C     http://sites.google.com/site/susanschommersite/software

C

C     and periodically updated versions ALGENCAN at

C     the TANGO home page:

C

C     www.ime.usp.br/~egbirgin/tango/

C     ******************************************************************

C     ******************************************************************

C     GE model

C     ----------------------------

C     Given a endowment (endow), durability technologies (yc) and probability of 

C     states (prob), the GE system is solved with logarithmic utility function

C

C      u^h = log(x^h_{0l})+ \sum_{s \in S}\ prob_s log(x^h_{sl})

C

C     [Sc] S. Schommer, "Computing general economic equilibrium with

C     incomplete markets, collateral and default penalties",

C     Working Paper, 2008

C

C     ******************************************************************

C     ******************************************************************

      subroutine inip(n,x,l,u,m,lambda,equatn,linear,coded,checkder)

      implicit none

C     This subroutine must set some problem data. For achieving this

C     objective YOU MUST MODIFY it according to your problem. See below

C     where your modifications must be inserted.

C    

C     Parameters of the subroutine:

C

C     On Entry:

C

C     This subroutine has no input parameters.

C

C     On Return

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              initial point,

C

C     l        double precision l(n),

C              lower bounds on x,

C

C     u        double precision u(n),

C              upper bounds on x,

C

C     m        integer,

C              number of constraints (excluding the bounds),

C

C     lambda   double precision lambda(m),

C              initial estimation of the Lagrange multipliers,

C

C     equatn   logical equatn(m)

C              for each constraint j, set equatn(j) = .true. if it is an

C              equality constraint of the form c_j(x) = 0, and set

C              equatn(j) = .false. if it is an inequality constraint of

C              the form c_j(x) <= 0,

C

C     linear   logical linear(m)

C              for each constraint j, set linear(j) = .true. if it is a

C              linear constraint, and set linear(j) = .false. if it is a

C              nonlinear constraint.

C     SCALAR ARGUMENTS

      integer m,n

      logical checkder

C     ARRAY ARGUMENTS

      logical coded(10),equatn(*),linear(*)

      double precision l(*),lambda(*),u(*),x(*),random,mult

C     LOCAL SCALARS

      integer i

      double precision multiplier,tmp

      double precision c

C     ******************************************************************

C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR PROBLEM

C     DATA:

C     ******************************************************************

C     Number of variables

      integer agent,states,good,asset,migual,mineq

      integer contagent,contstates

      integer contgood,contgood2,contasset,truen

      integer magain

      double precision yc(10,10,10),endow(10,10,10)

      double precision prob(10)

      common /base/agent,states,good,asset

      common /more/yc,endow,prob

C     write in "data" agent,states,good,asset

      open (unit=14,file="data")

      read (14,*) agent,states,good,asset

      close(14)

C     truen = number of variables (n)

      truen = agent+agent*states*good+states*good

C     migual = number of equal constraints

      migual=agent*states*good+agent

     +        +states*good+1

C     mineq = number of inequality constraints

C     mineq = 0

      mineq=0

C     Initial point

      open(unit = 1,file="input")

      do i = 1,truen

          read(1,*) x(i),multiplier

      x(i) = x(i)*multiplier

      end do

      n=truen

      close(1)

      open(unit = 1,file="randomizeinput")

      do i = 1,truen

          read(1,*) random,multiplier

      x(i) = x(i)+random*multiplier

      end do

      close(1)

      open(unit = 29,file="inputyc")

      open(unit = 23,file="inputendow")

      open(unit = 27,file="inputprob")

      open(unit = 32,file="randomizeendow")

      do contstates=2,states

      do contgood=1,good

      do contgood2=1,good

      read (29,*) tmp,multiplier

      yc(contstates,contgood,contgood2)=tmp*multiplier

      end do

      end do

      end do

      do contagent=1,agent

      do contstates=1,states

      do contgood=1,good

      read (23,*) tmp,multiplier

      read (32,*) random,mult

      endow(contagent,contstates,contgood)=max(0.0d0,tmp

     +        *multiplier+(random-0.5)*mult)

      end do

      end do

      end do

      prob(1)=1

      do contstates=2,states

      read (27,*) tmp,multiplier

      prob(contstates)=tmp*multiplier

      end do

      close(23)

      close(27)

      close(29)

      close(32)

C     Lower and upper bounds

      do i = 1,n

          l(i) =    0.0d0

          u(i) =   1d10

      end do

C     Number of constraints (equalities plus inequalities)

      m = migual+mineq

C     Lagrange multipliers approximation. Most users prefer to use the

C     null initial Lagrange multipliers estimates. However, if the

C     problem that you are solving is "slightly different" from a

C     previously solved problem of which you know the correct Lagrange

C     multipliers, we encourage you to set these multipliers as initial

C     estimates. Of course, in this case you are also encouraged to use

C     the solution of the previous problem as initial estimate of the

C     solution. Similarly, most users prefer to use rho = 10 as initial

C     penalty parameters. But in the case mentioned above (good

C     estimates of solution and Lagrange multipliers) larger values of

C     the penalty parameters (say, rho = 1000) may be more useful. More

C     warm-start procedures are being elaborated.

      do i = 1,m

          lambda(i) =  0.0d0

      end do

C     For each constraint i, set equatn(i) = .true. if it is an equality

C     constraint of the form c_i(x) = 0, and set equatn(i) = .false. if

C     it is an inequality constraint of the form c_i(x) <= 0.

      do i = 1,mineq

          equatn(i) = .false.

      end do

      do i=mineq+1,m

      equatn(i) = .true.

      end do

C     For each constraint i, set linear(i) = .true. if it is a linear

C     constraint, otherwise set linear(i) = .false.

      do i = 1,m

          linear(i) = .false.

      end do

C     Indicate which subroutines did you code.

      coded( 1) = .true.  ! evalf   

      coded( 2) = .true.  ! evalg  

      coded( 3) = .true.  ! evalh  

      coded( 4) = .true.  ! evalc  

      coded( 5) = .true.  ! evaljac

      coded( 6) = .false.  ! evalhc 

      coded( 7) = .false. ! evalfc

      coded( 8) = .false. ! evalgjac

      coded( 9) = .false. ! evalhl

      coded(10) = .false. ! evalhlp

C     Set checkder = TRUE if you code some derivatives and you would

C     like them to be tested by finite differences. It is highly

C     recommended.

      checkder = .false.

      call writeconstraints(n,m,x,"initdata")

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE INIP.

C     ******************************************************************

      call randomjac(n,m)

      end

C     ******************************************************************

C     ******************************************************************

      subroutine randomjac(n,m)

C     This subroutine must compute the non zero coefficient's' jacobian.

C     For achieving these objective YOU MUST MODIFY it in the

C     way specified below. However, if you decide to use numerical

C     derivatives (we dont encourage this option at all!) you dont need

C     to modify evaljac

C

C

C     Parameters of the subroutine:

C

C     On Entry:

C

C     This subroutine has no input parameters.

C

C     On Return

C

C     n        integer,

C              number of variables,

C

C     m        integer,

C              number of constraints (excluding the bounds),

      integer n,m,nnzjacbase(5000),curjacbase(5000),indjacbase(50000)

      integer i,j,l,u,nnzhbase(5000),curhbase(5000)

      integer hclinbase(50000),hccolbase(50000)

      double precision d,y(n)

      common /indices/indjacbase,nnzjacbase,curjacbase

      common /indicesh/nnzhbase,hclinbase,hccolbase,curhbase

      do i=1,n

      y(i)=rand(0)

      end do

      u=0

      do i=1,m

      curjacbase(i)=u

      nnzjacbase(i)=0

      do j=1,n

      call diff(n,y,i,j,d)

      if (d .ne. 0) then

      u=u+1

      nnzjacbase(i)=nnzjacbase(i)+1

      indjacbase(u)=j

      end if

      end do

      end do

      u=0

      do i=1,m

      curhbase(i)=u

      nnzhbase(i)=0

      do j=1,nnzjacbase(i)

      do l=1,j

      call diff2(n,y,i,indjacbase(curjacbase(i)+j),indjacbase

     +        (curjacbase(i)+l),d)

      if (d ** 2 .ge. 1d-20) then

      u=u+1

      nnzhbase(i)=nnzhbase(i)+1

      hclinbase(u)=indjacbase(curjacbase(i)+j)

      hccolbase(u)=indjacbase(curjacbase(i)+l)

      end if

      end do

      end do

      end do

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalf(n,x,f,flag)

      implicit none

C     This subroutine must compute the objective function. For achieving

C     this objective YOU MUST MODIFY it according to your problem. See

C     below where your modifications must be inserted.

C    

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     On Return

C

C     f        double precision,

C              objective function value at x,

C

C     flag     integer,

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of the objective

C              function. (For example, trying to compute the square root

C              of a negative number, dividing by zero or a very small

C              number, etc.) If everything was o.k. you must set it

C              equal to zero.

C     SCALAR ARGUMENTS

      integer flag,n

      double precision f

C     ARRAY ARGUMENTS

      double precision x(n)

C     Objective function

C     ******************************************************************

C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR OBJECTIVE

C     FUNCTION:

C     ******************************************************************

      flag = 0

      f = 0

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALF.

C     ******************************************************************

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalg(n,x,g,flag)

      implicit none

C     This subroutine must compute the gradient vector of the objective

C     function. For achieving these objective YOU MUST MODIFY it in the

C     way specified below. However, if you decide to use numerical

C     derivatives (we dont encourage this option at all!) you dont need

C     to modify evalg.

C

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     On Return

C

C     g        double precision g(n),

C              gradient vector of the objective function evaluated at x,

C

C     flag     integer,

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of any component

C              of the gradient vector. (For example, trying to compute

C              the square root of a negative number, dividing by zero or

C              a very small number, etc.) If everything was o.k. you

C              must set it equal to zero.

C     SCALAR ARGUMENTS

      integer flag,n

C     ARRAY ARGUMENTS

      double precision g(n),x(n)

C     LOCAL SCALARS

      integer i

C     Gradient vector

C     ******************************************************************

C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE GRADIENT

C     VECTOR OF YOUR OBJECTIVE FUNCTION:

C     ******************************************************************

      flag = 0

      do i = 1,n - 1

          g(i) = 0.0d0

      end do

      g(n) = 0.0d0

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALG.

C     ******************************************************************

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalh(n,x,hlin,hcol,hval,nnzh,flag)

      implicit none

C     This subroutine might compute the Hessian matrix of the objective

C     function. For achieving this objective YOU MAY MODIFY it according

C     to your problem. To modify this subroutine IS NOT MANDATORY. See

C     below where your modifications must be inserted.

C    

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     On Return

C

C     nnzh     integer,

C              number of perhaps-non-null elements of the computed

C              Hessian,

C

C     hlin     integer hlin(nnzh),

C              see below,

C

C     hcol     integer hcol(nnzh),

C              see below,

C

C     hval     double precision hval(nnzh),

C              the non-null value of the (hlin(k),hcol(k)) position

C              of the Hessian matrix of the objective function must

C              be saved at hval(k). Just the lower triangular part of

C              Hessian matrix must be computed,

C

C     flag     integer,

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of the Hessian

C              matrix of the objective function. (For example, trying

C              to compute the square root of a negative number,

C              dividing by zero or a very small number, etc.) If

C              everything was o.k. you must set it equal to zero.

C     SCALAR ARGUMENTS

      integer flag,n,nnzh

C     ARRAY ARGUMENTS

      integer hcol(*),hlin(*)

      double precision hval(*),x(n)

C     ******************************************************************

C     FROM HERE ON YOU MAY (OPTIONALLY) MODIFY THE SUBROUTINE TO SET THE

C     HESSIAN MATRIX OF YOUR OBJECTIVE FUNCTION:

C     ******************************************************************

      flag = 0

      nnzh = 0

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALH.

C     ******************************************************************

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evald(n,x,j,c)

C     This subroutine must compute the ind-th constraint of your

C     problem. For achieving this objective YOU MUST MODIFY it

C     according to your problem. See below the places where your

C     modifications must be inserted.

C

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     j        integer,

C              index of the constraint to be computed,

C

C     On Return

C

C     c        double precision,

C              ind-th constraint evaluated at x,

C

C     flag     integer

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of the

C              constraint. (For example, trying to compute the square

C              root of a negative number, dividing by zero or a very

C              small number, etc.) If everything was o.k. you must set

C              it equal to zero.

      integer j,flag,n,i,omega

      double precision c

      integer agent,states,good,asset,migual,mineq

      double precision yc(10,10,10),endow(10,10,10)

      double precision prob(10)

C      double precision alpha(10,10)

      common /base/agent,states,good,asset

      common /more/yc,endow,prob

      integer contagent,contstates,contgood,contasset

      integer contstates2,contgood2

C     ARRAY ARGUMENTS

      double precision x(n)

      integer lambda0,consumption0,price0

      double precision tmp,tmp2

C     Constraints

C     ******************************************************************

C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR

C     CONSTRAINTS:

C     ******************************************************************

      flag = 0

      lambda0=0

      consumption0=lambda0+agent

      price0=consumption0+agent*states*good

      i=j

C     First-order conditions consumptions

      if (i .le. agent*states*good) then

      contagent=1

      contstates=1

      do while (i .gt. states*good)

      i=i-(states*good)

      contagent=contagent+1

      end do

      do while (i .gt. good)

      i=i-good

      contstates=contstates+1

      end do

      contgood=i

      if (contstates .eq. 1) then

      c=1-x(consumption0+(contagent-1)

     +        *states*good+contgood)

     +        *x(lambda0+contagent)

     +        *x(price0+contgood)

      do contstates2=2,states

      do contgood2=1,good

      c=c+x(consumption0+(contagent-1)*states*good+contgood)

     +        *x(lambda0+contagent)

     +        *yc(contstates2,contgood,contgood2)*x(price0

     +        +(contstates2-1)*good+contgood2)

      end do

      end do

      else

      c=prob(contstates)-x(consumption0

     +        +(contagent-1)*states

     +        *good+(contstates-1)*good+contgood)

     +        *x(lambda0+contagent)

     +        *x(price0+(contstates-1)*good+contgood)

      end if

      return

      else

      i=i-agent*states*good

      end if

C      The budget constraint

      if (i .le. agent) then

      contagent=i

      c=0

      do contgood=1,good

      do contstates=1,states

      c=c+x(price0+(contstates-1)*good+contgood)

     +        *(x(consumption0+(contagent-1)*states

     +        *good+(contstates-1)*good+contgood)

     +        -endow(contagent,contstates,contgood))

      end do

      do contstates=2,states

      do contgood2=1,good

      c=c-x(price0+(contstates-1)*good+contgood)

     +        *yc(contstates,contgood,contgood2)

     +        *x(consumption0+(contagent-1)*states

     +        *good+contgood2)

      end do

      end do

      end do

      return

      else

      i=i-agent

      end if

C     Equilibrium commodity markets to clear

      if (i .le. states*good) then

      contstates=1

      do while (i .gt. good)

      i=i-good

      contstates=contstates+1

      end do

      contgood=i

      if (contstates .eq. 1) then

      c=0

      do contagent=1,agent

      c=c+x(consumption0+(contagent-1)*states*good+contgood)

     +        -endow(contagent,contstates,contgood)

      end do

      else

      c=0

      do contagent=1,agent

      c=c+x(consumption0+(contagent-1)*states*good+(contstates-1)*good

     +        +contgood)

     +        -endow(contagent,contstates,contgood)

      do contgood2=1,good

      c=c-yc(contstates,contgood,contgood2)

     +        *x(consumption0+(contagent-1)*states*good+contgood2)

      end do

      end do

      end if

      return

      else

      i=i-states*good

      end if

C     The simplex condition for prices

      c=-1

      do contgood=1,good

      do contstates=1,states

      c=c+x(price0+(contstates-1)*good+contgood)

      end do

      end do

      return

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalc(n,x,ind,c,flag)

      implicit none

C     This subroutine must compute the ind-th constraint of evald

C     For achieving this objective YOU MUST MODIFY it

C     according to your problem. See below the places where your

C     modifications must be inserted.

C

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     ind      integer,

C              index of the constraint to be computed,

C

C     On Return

C

C     c        double precision,

C              ind-th constraint evaluated at x,

C

C     flag     integer

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of the

C              constraint. (For example, trying to compute the square

C              root of a negative number, dividing by zero or a very

C              small number, etc.) If everything was o.k. you must set

C              it equal to zero.

C     SCALAR ARGUMENTS

      integer ind,flag,n,i,j

      double precision c

      integer agent,states,good,asset,m,migual,mineq

      integer contagent,contstates,contgood,contasset

      common /base/agent,states,good,asset

C     ARRAY ARGUMENTS

      double precision x(n),d(10000),teste

C     Constraints

C     ******************************************************************

C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET YOUR

C     CONSTRAINTS:

C     ******************************************************************

      flag = 0

      call evald(n,x,ind,c)

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALC.

C     ******************************************************************

      end

C     ******************************************************************

C     ******************************************************************

      subroutine writeconstraints(n,m,x,filename)

C     This subroutine can be used to write variables and constraints

      integer agent,states,good,asset,omega,n,m

      double precision yc(10,10,10),endow(10,10,10)

      double precision prob(10)

      common /base/agent,states,good,asset

      common /more/yc,endow,prob

      integer contagent,contstates,contgood,contasset

      integer contgood2

      character filename*8

C     ARRAY ARGUMENTS

      double precision x(n),c

      double precision utility

      integer lambda0,consumption0,price0

      integer yc0,endow0,prob0

      lambda0=0

      consumption0=lambda0+agent

      price0=consumption0+agent*states*good

      open(unit = 5,file=filename)

C format , ,agent,states,good,asset,good2, ,

 1    format(f25.15,3x,3x,3x,3x,i3,3x,3x,a)

 2    format(f25.15,3x,3x,3x,i3,3x,3x,3x,a)

 3    format(f25.15,3x,3x,3x,i3,i3,3x,3x,a)

 4    format(f25.15,3x,3x,i3,3x,3x,3x,3x,a)

 5    format(f25.15,3x,3x,i3,3x,i3,3x,3x,a)

 6    format(f25.15,3x,3x,i3,i3,3x,3x,3x,a)

 7    format(f25.15,3x,3x,i3,i3,i3,3x,3x,a)

 8    format(f25.15,3x,i3,3x,3x,3x,3x,3x,a)

 9    format(f25.15,3x,i3,3x,3x,i3,3x,3x,a)

 10   format(f25.15,3x,i3,3x,i3,3x,3x,3x,a)

 11   format(f25.15,3x,i3,3x,i3,i3,3x,3x,a)

 12   format(f25.15,3x,i3,i3,3x,3x,3x,3x,a)

 13   format(f25.15,3x,i3,i3,3x,i3,3x,3x,a)

 14   format(f25.15,3x,i3,i3,i3,3x,3x,3x,a)

 15   format(f25.15,3x,i3,i3,i3,i3,3x,3x,a)

 16   format(f25.15,3x,3x,3x,3x,3x,3x,3x,a)

 17   format(f25.15,3x,3x,i3,i3,3x,i3,3x,a)

      omega=1

      do contagent=1,agent

      write (5,8) x(omega),contagent,"lambda"

      omega=omega+1

      end do

      do contagent=1,agent

      do contstates=1,states

      do contgood=1,good

      write (5,14) x(omega),contagent,contstates,contgood,"consumption"

      omega=omega+1

      end do

      end do

      end do

      do contstates=1,states

      do contgood=1,good

      write (5,6) x(omega),contstates,contgood,"price"

      omega=omega+1

      end do

      end do

      do contstates=2,states

      do contgood=1,good

      do contgood2=1,good

      write (5,17) yc(contstates,contgood,contgood2),

     +        contstates,contgood,contgood2,

     +        "yc"

      omega=omega+1

      end do

      end do

      end do

      do contagent=1,agent

      do contstates=1,states

      do contgood=1,good

      write (5,14) endow(contagent,contstates,contgood),contagent,

     +        contstates,contgood,

     +        "Endowment"

      omega=omega+1

      end do

      end do

      end do

      do contstates=2,states

      write (5,4) prob(contstates),contstates,"prob"

      omega=omega+1

      end do

      omega=1

      do contagent=1,agent

      do contgood=1,good

      do contstates=1,states

      call evald(n,x,omega,c)

      write(5,14) c,contagent,contstates,contgood,

     +        "Derutil"

      omega=omega+1

      end do

      end do

      end do

      do contagent=1,agent

      call evald(n,x,omega,c)

      write(5,8) c,contagent,

     +        "Restriction"

      omega=omega+1

      end do

      do contstates=1,states

      do contgood=1,good

      call evald(n,x,omega,c)

      write(5,6) c,contstates,contgood,

     +        "eqgood"

      omega=omega+1

      end do

      end do

      call evald(n,x,omega,c)

      write(5,16) c,

     +        "normalization"

      omega=omega+1

      do contagent=1,agent

      utility=0

      do contstates=1,states

      do contgood=1,good

      utility=utility+prob(contstates)

     +        *log(x(consumption0+(contagent-1)

     +        *states*good+(contstates-1)*good+contgood))

      end do

      end do

      write(5,8) utility,contagent,

     +        "Utility"

      end do

      close(5)

      end

C     ******************************************************************

C     ******************************************************************

      subroutine diff(n,x,ind,indjac,d)

C     This subroutine might compute the Jacobian

C      quadratic problem

      integer n,ind,indjac

      double precision c,d,x(n)

      x(indjac)=x(indjac)-1

      call evald(n,x,ind,c)

      d=c

      x(indjac)=x(indjac)+2

      call evald(n,x,ind,c)

      x(indjac)=x(indjac)-1

      d=(c-d)/2

      end

C     ******************************************************************

C     ******************************************************************

      subroutine diff2(n,x,ind,indlin,indrow,d)

C     This subroutine might compute the Hessian

C     quadratic problem

      integer n,ind,indlin,indrow

      double precision c,d,x(n),xl,xr

      xl=x(indlin)

      xr=x(indrow)

      x(indlin)=x(indlin)+0.5

      x(indrow)=x(indrow)+0.5

      call evald(n,x,ind,c)

      d=c

      x(indlin)=x(indlin)-1

      call evald(n,x,ind,c)

      d=d-c

      x(indrow)=x(indrow)-1

      call evald(n,x,ind,c)

      d=d+c

      x(indlin)=x(indlin)+1

      call evald(n,x,ind,c)

      x(indlin)=xl

      x(indrow)=xr

      d=d-c

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evaljac(n,x,ind,indjac,valjac,nnzjac,flag)

      implicit none

C     This subroutine must compute the gradient of the ind-th constraint.

C     For achieving these objective YOU MUST MODIFY it in the way

C     specified below.

C

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     ind      integer,

C              index of the constraint whose gradient will be computed,

C

C     On Return

C

C     nnzjac   integer,

C              number of perhaps-non-null elements of the computed

C              gradient,

C

C     indjac   integer indjac(nnzjac),

C              see below,

C

C     valjac   double precision valjac(nnzjac),

C              the non-null value of the partial derivative of the

C              ind-th constraint with respect to the indjac(k)-th

C              variable must be saved at valjac(k).

C

C     flag     integer

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of the

C              constraint. (For example, trying to compute the square

C              root of a negative number, dividing by zero or a very

C              small number, etc.) If everything was o.k. you must set

C              it equal to zero.

C     SCALAR ARGUMENTS

      integer flag,ind,n,m,nnzjac

C     ARRAY ARGUMENTS

      integer indjac(n),j

      double precision x(n),valjac(n)

      integer nnzjacbase(5000),indjacbase(50000)

      integer curjacbase(5000)

      common /indices/indjacbase,nnzjacbase,curjacbase

C     Sparse gradient vector of the ind-th constraint

C     ******************************************************************

C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO SET THE GRADIENTS 

C     OF YOUR CONSTRAINTS:

C     ******************************************************************

      flag = 0

      nnzjac=nnzjacbase(ind)

      do j=1,nnzjac

      indjac(j)=indjacbase(curjacbase(ind)+j)

      call diff(n,x,ind,indjac(j),valjac(j))

      end do

      end

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALJAC.

C     ******************************************************************

C     ******************************************************************

C     ******************************************************************

      subroutine evalhc(n,x,ind,hclin,hccol,hcval,nnzhc,flag)

      implicit none

C     This subroutine might compute the Hessian matrix of the ind-th

C     constraint. For achieving this objective YOU MAY MODIFY it

C     according to your problem. To modify this subroutine IS NOT

C     MANDATORY. See below where your modifications must be inserted.

C    

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     ind      integer,

C              index of the constraint whose Hessian will be computed,

C

C     On Return

C

C     nnzhc    integer,

C              number of perhaps-non-null elements of the computed

C              Hessian,

C

C     hclin    integer hclin(nnzhc),

C              see below,

C

C     hccol    integer hccol(nnzhc),

C              see below,

C

C     hcval    double precision hcval(nnzhc),

C              the non-null value of the (hclin(k),hccol(k)) position

C              of the Hessian matrix of the ind-th constraint must

C              be saved at hcval(k). Just the lower triangular part of

C              Hessian matrix must be computed,

C

C     flag     integer,

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of the Hessian

C              matrix of the ind-th constraint. (For example, trying

C              to compute the square root of a negative number,

C              dividing by zero or a very small number, etc.) If

C              everything was o.k. you must set it equal to zero.

C     SCALAR ARGUMENTS

      integer flag,ind,n,nnzhc,i

      integer nnzhbase(5000),hclinbase(50000),hccolbase(50000)

      integer curhbase(5000)

C     ARRAY ARGUMENTS

      integer hccol(*),hclin(*)

      double precision hcval(*),x(n)

      common /indicesh/nnzhbase,hclinbase,hccolbase,curhbase

C     ******************************************************************

C     FROM HERE ON YOU MAY (OPTIONALLY) MODIFY THE SUBROUTINE TO SET THE

C     HESSIAN'S OF YOUR CONSTRAINTS:

C     ******************************************************************

      flag = 0

      nnzhc = nnzhbase(ind)

      do i=1,nnzhc

      hccol(i)=hccolbase(curhbase(ind)+i)

      hclin(i)=hclinbase(curhbase(ind)+i)

      call diff2(n,x,ind,hclin(i),hccol(i),hcval(i))

      end do

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE EVALHC.

C     ******************************************************************

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalfc(n,x,f,m,c,flag)

      implicit none

C     SCALAR ARGUMENTS

      integer flag,m,n

      double precision f

C     ARRAY ARGUMENTS

      double precision c(m),x(n)

      flag = - 1

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalgjac(n,x,g,m,jcfun,jcvar,jcval,jcnnz,flag)

      implicit none

C     SCALAR ARGUMENTS

      integer flag,jcnnz,m,n

C     ARRAY ARGUMENTS

      integer jcfun(*),jcvar(*)

      double precision g(n),jcval(*),x(n)

      flag = - 1

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalhl(n,x,m,lambda,scalef,scalec,hllin,hlcol,hlval,

     +hlnnz,flag)

      implicit none

C     SCALAR ARGUMENTS

      integer flag,hlnnz,m,n

      double precision scalef

C     ARRAY ARGUMENTS

      integer hlcol(*),hllin(*)

      double precision hlval(*),lambda(m),scalec(m),x(n)

      flag = - 1

      end

C     ******************************************************************

C     ******************************************************************

      subroutine evalhlp(n,x,m,lambda,sf,sc,p,hp,goth,flag)

      implicit none

C     SCALAR ARGUMENTS

      logical goth

      integer flag,m,n

      double precision sf

C     ARRAY ARGUMENTS

      double precision hp(n),lambda(m),p(n),sc(m),x(n)

C     This subroutine might compute the product of the Hessian of the

C     Lagrangian times vector p (just the Hessian of the objective

C     function in the unconstrained or bound-constrained case).

C    

C     Parameters of the subroutine:

C

C     On Entry:

C

C     n        integer,

C              number of variables,

C

C     x        double precision x(n),

C              current point,

C

C     m        integer,

C              number of constraints,

C

C     lambda   double precision lambda(m),

C              vector of Lagrange multipliers,

C

C     p        double precision p(n),

C              vector of the matrix-vector product,

C

C     goth     logical,

C              can be used to indicate if the Hessian matrices were

C              computed at the current point. It is set to .false.

C              by the optimization method every time the current

C              point is modified. Suggestion: if its value is .false.

C              then compute the Hessian's, save them in a common

C              structure and set goth to .true.. Otherwise, just use

C              the Hessian's saved in the common block structure,

C

C     On Return

C

C     hp       double precision hp(n),

C              Hessian-vector product,

C

C     goth     logical,

C              see above,

C             

C     flag     integer,

C              You must set it to any number different of 0 (zero) if

C              some error occurred during the evaluation of the

C              Hessian-vector product. (For example, trying to compute

C              the square root of a negative number, dividing by zero

C              or a very small number, etc.) If everything was o.k. you

C              must set it equal to zero.

      flag = - 1

      end

C     ******************************************************************

C     ******************************************************************

      subroutine endp(n,x,l,u,m,lambda,equatn,linear)

      implicit none

C     This subroutine can be used to do some extra job after the solver

C     has found the solution,like some extra statistics, or to save the

C     solution in some special format or to draw some graphical

C     representation of the solution. If the information given by the

C     solver is enough for you then leave the body of this subroutine

C     empty.

C    

C     Parameters of the subroutine:

C

C     The parameters of this subroutine are the same parameters of

C     subroutine inip. But in this subroutine there are not output

C     parameter. All the parameters are input parameters.

C     SCALAR ARGUMENTS

      integer m,n,omega

C     ARRAY ARGUMENTS

      logical equatn(m),linear(m)

      double precision l(n),lambda(m),u(n),x(n)

C     ******************************************************************

C     FROM HERE ON YOU MUST MODIFY THE SUBROUTINE TO COMPLEMENT THE

C     INFORMATION RELATED TO THE SOLUTION GIVEN BY THE SOLVER

C     ******************************************************************

C     ******************************************************************

C     STOP HERE YOUR MODIFICATIONS OF SUBROUTINE ENDP

C     ******************************************************************

      call writeconstraints(n,m,x,"solution")

C     write in "newinput" solution for initial point

      open (unit=11,file="newinput")

      do omega=1,n

      write (11,21) x(omega),",1"

      end do

 21   format(d20.15,a)

      end