GEIsystem.f

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

C File: GEIsystem.f

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 http://sites.google.com/site/susanschommersite/software

C and periodically updated versions ALGENCAN at

C the TANGO home page:

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

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

C GEI model

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

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

C states (prob) and promise in the asset, the GEI system is solved with

C logarithmic utility function

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

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

C incomplete markets, collateral and default penalties",

C Working Paper, 2008

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 Parameters of the subroutine:

C On Entry:

C This subroutine has no input parameters.

C On Return

C n integer,

C number of variables,

C x double precision x(n),

C initial point,

C l double precision l(n),

C lower bounds on x,

C u double precision u(n),

C upper bounds on x,

C m integer,

C number of constraints (excluding the bounds),

C lambda double precision lambda(m),

C initial estimation of the Lagrange multipliers,

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 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,truenposit,truennegat

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

double precision promise(10,10,10),prob(10)

common /base/agent,states,good,asset

common /more/yc,endow,promise,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)

truenposit=agent*states+agent*states*good+states*good+asset

truennegat=agent*asset

truen = truenposit+truennegat

C migual = number of equal constraints

migual=agent*states*good+agent*states+agent*asset

+ +states*good+asset+states

C mineq = number of inequality constraints

mineq=0

C Initial point

C write in "input" initial point,multiplier (1)

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)

C write in "randomizeinput" random initial point,multiplier (10)

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

do i = 1,truen

read(1,*) random,multiplier

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

end do

close(1)

C Input

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

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

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

open(unit = 24,file="inputpromise")

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

do contstates=2,states

do contgood=1,good

do contgood2=1,good

read (29,*) tmp,multiplier

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

end do

do contagent=1,agent

do contstates=1,states

do contgood=1,good

read (23,*) tmp,multiplier

read (32,*) random,mult

C endow(contagent,contstates,contgood)=tmp*multiplier

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

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

end do

do contstates=2,states

do contgood=1,good

do contasset=1,asset

read (24,*) tmp,multiplier

promise(contstates,contgood,contasset)=tmp*multiplier

end do

prob(1)=1

do contstates=2,states

read (27,*) tmp,multiplier

prob(contstates)=tmp*multiplier

end do

close(23)

close(24)

close(27)

close(29)

close(32)

C Lower and upper bounds

do i = 1,truenposit

l(i) = 0.0d0

end do

do i = truenposit+1,truen

l(i) = -100

end do

do i = 1,n

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 call writeconstraints(n,m,x,"initdata")

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.

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

C STOP HERE YOUR MODIFICATIONS OF SUBROUTINE INIP.

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

call randomjac(n,m)

end

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 Parameters of the subroutine:

C On Entry:

C This subroutine has no input parameters.

C On Return

C n integer,

C number of variables,

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

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

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 Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C On Return

C f double precision,

C objective function value at x,

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 ******************************************************************

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 Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C On Return

C g double precision g(n),

C gradient vector of the objective function evaluated at x,

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 ******************************************************************

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 Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C On Return

C nnzh integer,

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

C Hessian,

C hlin integer hlin(nnzh),

C see below,

C hcol integer hcol(nnzh),

C see below,

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 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 ******************************************************************

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 MOFIFY it

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

C modifications must be inserted.

C Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C j integer,

C index of the constraint to be computed,

C On Return

C c double precision,

C ind-th constraint evaluated at x,

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 promise(10,10,10),prob(10)

common /base/agent,states,good,asset

common /more/yc,endow,promise,prob

integer contagent,contstates,contgood,contasset

integer contstates2,contgood2

C ARRAY ARGUMENTS

double precision x(n)

integer lambda0,consumption0,price0,portfolio0,priceasset0

double precision modifierdur,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*states

price0=consumption0+agent*states*good

priceasset0=price0+states*good

portfolio0=priceasset0+asset

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-1)*states+1)

+ *x(price0+contgood)

do contstates2=2,states

do contgood2=1,good

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

+ *x(lambda0+(contagent-1)*states+contstates2)

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

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

end do

else

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

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

+ *x(lambda0+(contagent-1)*states+contstates)

+ *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*states) then

contagent=1

do while (i .gt. states)

i=i-states

contagent=contagent+1

end do

contstates=i

if (contstates .eq. 1) then

c=0

do contgood=1,good

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

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

+ *good+contgood)-endow(contagent,1,contgood))

end do

do contasset=1,asset

c=c+x(priceasset0+contasset)

+ *x(portfolio0+(contagent-1)*asset+contasset)

end do

else

c=0

do contgood=1,good

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

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

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

+ -endow(contagent,contstates,contgood))

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

do contasset=1,asset

do contgood=1,good

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

+ *promise(contstates,contgood,contasset)

+ *x(portfolio0+(contagent-1)*asset+contasset)

end do

end if

return

else

i=i-agent*states

end if

C First-order conditions portfolio

if (i .le. agent*asset) then

contagent=1

do while (i .gt. asset)

i=i-asset

contagent=contagent+1

end do

contasset=i

c=-x(lambda0+(contagent-1)*states+1)*x(priceasset0+contasset)

do contgood=1,good

do contstates=2,states

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

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

+ *promise(contstates,contgood,contasset)

end do

return

else

i=i-agent*asset

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 if

return

else

i=i-states*good

end if

C Equilibrium asset markets to clear

if (i .le. asset) then

c=0

contasset=i

do contagent=1,agent

c=c+x(portfolio0+(contagent-1)*asset+contasset)

end do

return

else

i=i-asset

end if

C The simplex condition for prices

if (i .le. states) then

contstates=i

c=-1

do contgood=1,good

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

end do

return

else

i=i-states

end if

end

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 MOFIFY it

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

C modifications must be inserted.

C Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C ind integer,

C index of the constraint to be computed,

C On Return

C c double precision,

C ind-th constraint evaluated at x,

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 ******************************************************************

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 promise(10,10,10),prob(10)

common /base/agent,states,good,asset

common /more/yc,endow,promise,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,portfolio0,priceasset0

integer yc0,endow0,promise0,prob0

lambda0=0

consumption0=lambda0+agent*states

price0=consumption0+agent*states*good

priceasset0=price0+states*good

portfolio0=priceasset0+asset

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

do contstates=1,states

write (5,12) x(omega),contagent,contstates,"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

do contstates=1,states

do contgood=1,good

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

omega=omega+1

end do

do contasset=1,asset

write (5,1) x(omega),contasset,"priceasset"

omega=omega+1

end do

do contagent=1,agent

do contasset=1,asset

write (5,9) x(omega),contagent,contasset,"portfolio"

omega=omega+1

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

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

do contstates=2,states

do contgood=1,good

do contasset=1,asset

write (5,7) promise(contstates,contgood,contasset),contstates,

+ contgood,contasset,

+ "Promise"

omega=omega+1

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

do contagent=1,agent

do contstates=1,states

call evald(n,x,omega,c)

write(5,12) c,contagent,contstates,

+ "Restriction"

omega=omega+1

end do

do contagent=1,agent

do contasset=1,asset

call evald(n,x,omega,c)

write(5,9) c,contagent,contasset,

+ "Derportfolio"

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

do contasset=1,asset

call evald(n,x,omega,c)

write(5,1) c,contasset,

+ "eqasset"

omega=omega+1

end do

do contstates=1,states

call evald(n,x,omega,c)

write(5,4) c,contstates,

+ "normalization"

omega=omega+1

end do

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

write(5,8) utility,contagent,

+ "Utility"

end do

close(5)

end

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 ******************************************************************

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 ******************************************************************

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 Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C ind integer,

C index of the constraint whose gradient will be computed,

C On Return

C nnzjac integer,

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

C gradient,

C indjac integer indjac(nnzjac),

C see below,

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 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 ******************************************************************

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 Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C ind integer,

C index of the constraint whose Hessian will be computed,

C On Return

C nnzhc integer,

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

C Hessian,

C hclin integer hclin(nnzhc),

C see below,

C hccol integer hccol(nnzhc),

C see below,

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 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 ******************************************************************

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 ******************************************************************

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 ******************************************************************

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 ******************************************************************

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 Parameters of the subroutine:

C On Entry:

C n integer,

C number of variables,

C x double precision x(n),

C current point,

C m integer,

C number of constraints,

C lambda double precision lambda(m),

C vector of Lagrange multipliers,

C p double precision p(n),

C vector of the matrix-vector product,

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 On Return

C hp double precision hp(n),

C Hessian-vector product,

C goth logical,

C see above,

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 ******************************************************************

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 Parameters of the subroutine:

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 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