CFPP$ UNROLL (32) F c c TITLE: implementation of the finite-element method for model c with capacity utilization and moving costs, in ``The c Effects of Government purchases on Employment and c Hours,'' (Toni Braun, Ellen McGrattan). c c PROBLEM: find c1^l(k,e,z), n^l(k,e,z), c0^l(k,e,z), and c h^l(k,e,z) that solve the functional equations: c c U_c(c1^l,1-h^l)(1-a^l) - beta sum_{j=1:n} pi(l,j) c * int U_c(c1^l(x~),1-h^l(x~))(1-a^j) c *(1-delta+z~*F_k(k~,n^j(x~),h^j(x~))) c *f(eps) d(eps) = 0 (1) c U(c1^l,1-h^l) - U(c0^l,1)-P'(n^l)-U_c(c1^l,1-h^l) c *(c1^l-c0^l-z*F(k,n^l,h^l)+2*alpha c *(n^l-e))*(1-a^l)-2*alpha*beta c *sum_{j=1:n} pi(l,j) int (1-a^j) c *U_c(c1^j(x~),1-h^j(x~))*(n^j(x~)-e~) c *f(eps) d(eps) = 0, (2) c U_c(c0^l,1) - U_c(c1^l,1-h^l) = 0, (3) c n^l*U_l(c1^l,1-h^l) - U_c(c1^l,1-h^l)*F_h(k,n^l,h^l)=0 (4) c c for l=1,..,ns, where x~=(k~,e~,z~) and c k~ = z*F(k,n^l(k,e,z),h^l(k,e,z))-n^l(k,e,z)* c c1^l(k,e,z)-(1-n^l(k,e,z))*c0^l(k,e,z) c +(1-delta)*k-(cg^l+zeta*log(z)) c e~ = n^l(k,e,z) c z~ = z^psi*eps c c1(0,e,z) = 0 c c ALGORITHM: the Galerkin method with piecewise trilinear shape c functions N_a(x,y,z), e.g., c c1(x,y,z) = sum N_a(x,y,z)*c_a, c_a constant c n(x,y,z) = sum N_a(x,y,z)*n_a, n_a constant c c MODEL: stochastic growth model with government capital, c varying hours and employment, and moving costs; c the planner chooses c1[t],c0[t],i[t],n[t],h[t] c to maximize c c E[ sum beta^t {n[t]*U(c1[t],1-h[t])+(1-n[t]) c t *U(c0[t],1)-P(n[t])}*(1-a[t])|x[0]] c c subject to c c c[t]+i[t]+m[t] = z[t]*F(k[t]+kg[t],n[t],h[t])-ig[t]-cg[t] c c[t] = n[t]*c1[t]+(1-n[t])*c0[t] c m[t] = M(n[t],n[t-1]) c k[t+1] = (1-delta)*k[t]+i[t] c kg[t+1] = (1-delta)*kg[t]+ig[t] c z[t+1] = z[t]^psi*eps[t] c x[t] = (k[t],n[t-1],z[t]) c c FUNCTIONS: U(c,l) = [(c^gamma*l^(1-gamma))^(1-omega)-1]/(1-omega) c F(k,kg,n,h) = lam*z*(phi*k^rho+(1-phi)*kg^rho)^ ... c (theta/rho)*n^(1-theta)*h c M(n,e) = alpha*(n-e)*(n-e) c P(n) = [(n-pa)*pb]^3+abs((n-pa)*pb)^3 c c TEST CASE: set alpha=0 and compare results to Kydland-Prescott c c INPUTS: (a) cumc3da.inp: alpha,beta,delta,gamma,lam,omega, ... c pa,pb,phi,psi,rho,sig,theta c ig(j),cg(j),a(j), j=1:ns c mcp(i,j) i=1:ns, j=1:ns c c1(i,j) i=1:na, j=1:ns c n(i,j) i=1:na, j=1:ns c xa(i),nxp(i) i=1:nx c xa(nx+1) c ya(i),nyp(i) i=1:ny c ya(ny+1) c za(i),nzp(i) i=1:nz c za(nz+1) c where c1's,n's are initial guesses for policy rules c ig's,cg's,a's are values for Markov chain c mcp's are the transition probabilities c xa's,ya's,za's are grid pts. c nwp(i) are # of points per interval for c integrating over `w' interval, w=x,y,z c c OUTPUTS: (b) cumc3db.out: general results c (c) cumc3db.nxt: new input file c (d) cumc3d.dat: input file for cumc3d.m (Matlab file) c c SUBROUTINES: (LINPACK) blas.f dgefa.f, dgesl.f, dgeco.f c (KMN) dgefs.f, machcon.f, xerror.f c (McGrattan) qgausl.f c c REFERENCES: Eischen, Jeff, (notes on the Finite Element Method) c Kahaner, Moler, and Nash (1989), NUMERICAL METHODS c AND SOFTWARE (KMN) c Saad, Youcef (1993), ``Lecture notes: Scientific c Computation 8002,'' Univ. Minnesota c c SEE ALSO: cumc3da.f (which solves Ax=b assuming A is dense) c c Ellen McGrattan, 7-6-93 c Revised, 7-9-93, ERM c parameter (na=144,nx=15,ny=2,nz=2,ns=4,nas=1152,nxm=3,nym=3, c nzm=3,ne1=10,nw=10000,ldf=60000) implicit double precision (a-h,o-z) double precision f(nas), df(ldf), work(nw), c xa(nx+1), ya(ny+1), za(nz+1), xx(nxm,nx), wx(nxm,nx), c yy(nym,ny), wy(nym,ny), zz(nzm,nz), wz(nzm,nz), c x(nxm), u(nxm), y(nym), v(nym), z(nzm), w(nzm), c a0(na,ns), a1(na,ns), b0(na,ns), b1(na,ns), zzt(ne1), c e1(ne1), w1(ne1), wgt(ne1), sum(34), wrk(ne1*ns,32), c mcp(ns,ns), ginv(ns), cg(ns), army(ns), lam, n0, n0t, c mu, wrk2(nas) integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), kk(ne1), c ind(ne1*ns), ja(ldf), ia(ldf), nfound(8) c if (nas.gt.3000) then write(7,*) 'Increase nn in arnoldi.f to ', nas, c ' and update this message in cumc3db.f' write(*,*) 'Increase nn in arnoldi.f to ', nas, c ' and update this message in cumc3db.f' endif tiny = 0.0000001d0 pi = 3.14159265d0 crit = 0.000001d0 maxit = 20 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 open(unit=5, file='cumc3db.inp') open(unit=7, file='cumc3db.out') open(unit=9, file='cumc3db.nxt') open(unit=10,file='cumc3d.dat') read(5,*) alpha,beta,delta,gamma,lam,omega,pa,pb,phi,psi, c rho,sig,theta,zeta do 5 i=1,ns read(5,*) ginv(i),cg(i),army(i) 5 continue do 10 i=1,ns read(5,*) (mcp(i,j),j=1,ns) 10 continue do 20 i=1,na read(5,*) (a0(i,j),j=1,ns) 20 continue do 25 i=1,na read(5,*) (b0(i,j),j=1,ns) 25 continue do 30 i=1,nx read(5,*) xa(i),nxp(i) 30 continue read(5,*) xa(nx1) do 40 i=1,ny read(5,*) ya(i),nyp(i) 40 continue read(5,*) ya(ny1) do 50 i=1,nz read(5,*) za(i),nzp(i) 50 continue read(5,*) za(nz1) c delt1 = 1.d0-delta gamm1 = 1.d0-gamma ome1 = 1.d0-omega rho1 = 1.d0-rho thet1 = 1.d0-theta cwgt = gamma*ome1 cwgt1 = cwgt-1.d0 hwgt = gamm1*ome1 hwgt1 = hwgt-1.d0 sum1 = 0.d0 sum2 = 0.d0 call qgausl(ne1,-2.88d0*sig,2.88d0*sig,e1,w1) do 60 i=1,ne1 wgt(i) = dexp(-0.5d0*e1(i)*e1(i)/(sig*sig))/ c dsqrt(2.d0*pi*sig*sig)*w1(i) sum1 = sum1 + wgt(i) 60 continue do 70 i=1,ne1 wgt(i) = wgt(i)/sum1 sum2 = sum2 + wgt(i) 70 continue do 90 i=1,nx call qgausl(nxp(i),xa(i),xa(i+1),x,u) do 90 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 90 continue do 100 i=1,ny call qgausl(nyp(i),ya(i),ya(i+1),y,v) do 100 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 100 continue do 110 i=1,nz call qgausl(nzp(i),za(i),za(i+1),z,w) do 110 j=1,nzp(i) zz(j,i) = z(j) wz(j,i) = w(j) 110 continue ne = nx*ny*nz nxy = nx*ny dd = 1.d0 do 120 it=1,maxit do 130 i=1,nas f(i) = 0.d0 130 continue do 135 i=1,ldf df(i) = 0.d0 ia(i) = 0 ja(i) = 0 135 continue ndf = 0 do 140 l=1,ns army1 = 1.d0-army(l) do 150 n=1,ne iz = (n-1)/nxy+1 iy = (n-(iz-1)*nxy-1)/nx+1 ix = n-(iy-1)*nx-(iz-1)*nxy n1 = nxp(ix) n2 = nyp(iy) n3 = nzp(iz) x1 = xa(ix) x2 = xa(ix+1) dx = x2-x1 y1 = ya(iy) y2 = ya(iy+1) dy = y2-y1 z1 = za(iz) z2 = za(iz+1) dz = z2-z1 nd1= ix+nx1*(iy-1)+nx1*ny1*(iz-1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 220 i=1,n1 x(i) = xx(i,ix) u(i) = wx(i,ix) 220 continue do 230 j=1,n2 y(j) = yy(j,iy) v(j) = wy(j,iy) 230 continue do 240 k=1,n3 z(k) = zz(k,iz) w(k) = wz(k,iz) zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz do 240 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 240 i=1,n1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c1 = a0(nd1,l)*xp*yp*zp+a0(nd2,l)*xm*yp*zp+ c a0(nd3,l)*xm*ym*zp+a0(nd4,l)*xp*ym*zp+ c a0(nd5,l)*xp*yp*zm+a0(nd6,l)*xm*yp*zm+ c a0(nd7,l)*xm*ym*zm+a0(nd8,l)*xp*ym*zm n0 = b0(nd1,l)*xp*yp*zp+b0(nd2,l)*xm*yp*zp+ c b0(nd3,l)*xm*ym*zp+b0(nd4,l)*xp*ym*zp+ c b0(nd5,l)*xp*yp*zm+b0(nd6,l)*xm*yp*zm+ c b0(nd7,l)*xm*ym*zm+b0(nd8,l)*xp*ym*zm h = 1-gamm1*c1*(n0/x(i))**theta/(z(k)*lam*gamma) c0 = c1*(1.d0-h)**(hwgt/cwgt1) h1 = 1.d0-h mu = gamma*c1**cwgt1*h1**hwgt cg1 = cg(l)+zeta*dlog(z(k)) xt = delt1*x(i)+z(k)*lam*x(i)**theta*n0**thet1*h- c cg1-n0*c1-(1.d0-n0)*c0-alpha*(n0-y(j))* c (n0-y(j)) yt = n0 ii = 1 do 250 m=1,nx j1 = idint( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 250 continue ii = max(1,ii-1) jj = 1 do 260 m=1,ny j1 = idint( dmax1(yt-ya(m),0.d0)/(yt-ya(m))+.2) jj = jj + j1 260 continue jj = max(1,jj-1) CFPP$ UNROLL do 270 m=1,ne1 zzt(m) = dexp(psi*dlog(z(k))+e1(m)) kk(m) = 1 do 280 mm=1,nz j1 = idint( dmax1(zzt(m)-za(mm),0.d0)/ c (zzt(m)-za(mm))+.2) kk(m)= kk(m) + j1 280 continue kk(m) = max(1,kk(m)-1) 270 continue do 290 m=1,34 sum(m) = 0.d0 290 continue do 300 m=1,nas work(m) = 0.d0 wrk2(m) = 0.d0 300 continue d = (n0-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb dc0ca = -c0/c1*omega/cwgt1 dc0na = hwgt/cwgt1*c0*theta/n0 dhca = -h1/c1 dhna = -theta*h1/n0 dxca = z(k)*lam*x(i)**theta*n0**thet1*dhca-n0- c (1.d0-n0)*dc0ca dxna = z(k)*lam*x(i)**theta*n0**thet1*h*(thet1/ c n0+dhna/h)-c1+c0-(1.d0-n0)*dc0na-2.d0* c alpha*(n0-y(j)) do 310 mm=1,ne1 do 310 m=1,ns nct = (m-1)*ne1+mm zt = zzt(mm) nn = (kk(mm)-1)*nx*ny+(jj-1)*nx+ii x1t = xa(ii) x2t = xa(ii+1) dxt = x2t-x1t y1t = ya(jj) y2t = ya(jj+1) dyt = y2t-y1t z1t = za(kk(mm)) z2t = za(kk(mm)+1) dzt = z2t-z1t xpt = (x2t-xt)/dxt xmt = (xt-x1t)/dxt ypt = (y2t-yt)/dyt ymt = (yt-y1t)/dyt zpt = (z2t-zt)/dzt zmt = (zt-z1t)/dzt nd1t = (kk(mm)-1)*nx1*ny1+(jj-1)*nx1+ii nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 nd5t = nd1t+nx1*ny1 nd6t = nd5t+1 nd8t = nd5t+nx1 nd7t = nd8t+1 ind(nct) = nd1t c1t = a0(nd1t,m)*xpt*ypt*zpt+ c a0(nd2t,m)*xmt*ypt*zpt+ c a0(nd3t,m)*xmt*ymt*zpt+ c a0(nd4t,m)*xpt*ymt*zpt+ c a0(nd5t,m)*xpt*ypt*zmt+ c a0(nd6t,m)*xmt*ypt*zmt+ c a0(nd7t,m)*xmt*ymt*zmt+ c a0(nd8t,m)*xpt*ymt*zmt n0t = b0(nd1t,m)*xpt*ypt*zpt+ c b0(nd2t,m)*xmt*ypt*zpt+ c b0(nd3t,m)*xmt*ymt*zpt+ c b0(nd4t,m)*xpt*ymt*zpt+ c b0(nd5t,m)*xpt*ypt*zmt+ c b0(nd6t,m)*xmt*ypt*zmt+ c b0(nd7t,m)*xmt*ymt*zmt+ c b0(nd8t,m)*xpt*ymt*zmt ht = 1.d0-gamm1*c1t*(n0t/xt)**theta/ c (gamma*zt*lam) d = (n0t-pa)*pb dpnt = 3.d0*d*(d+dsign(d,0.d0))*pb d2pnt = 6.d0*(d+dsign(d,0.d0))*pb*pb dc1tx = (-a0(nd1t,m)*ypt*zpt+a0(nd2t,m)*ypt*zpt c +a0(nd3t,m)*ymt*zpt-a0(nd4t,m)*ymt*zpt c -a0(nd5t,m)*ypt*zmt+a0(nd6t,m)*ypt*zmt c +a0(nd7t,m)*ymt*zmt-a0(nd8t,m)*ymt*zmt) c /dxt dc1ty = (-a0(nd1t,m)*xpt*zpt-a0(nd2t,m)*xmt*zpt c +a0(nd3t,m)*xmt*zpt+a0(nd4t,m)*xpt*zpt c -a0(nd5t,m)*xpt*zmt-a0(nd6t,m)*xmt*zmt c +a0(nd7t,m)*xmt*zmt+a0(nd8t,m)*xpt*zmt) c /dyt dn0tx = (-b0(nd1t,m)*ypt*zpt+b0(nd2t,m)*ypt*zpt c +b0(nd3t,m)*ymt*zpt-b0(nd4t,m)*ymt*zpt c -b0(nd5t,m)*ypt*zmt+b0(nd6t,m)*ypt*zmt c +b0(nd7t,m)*ymt*zmt-b0(nd8t,m)*ymt*zmt) c /dxt dn0ty = (-b0(nd1t,m)*xpt*zpt-b0(nd2t,m)*xmt*zpt c +b0(nd3t,m)*xmt*zpt+b0(nd4t,m)*xpt*zpt c -b0(nd5t,m)*xpt*zmt-b0(nd6t,m)*xmt*zmt c +b0(nd7t,m)*xmt*zmt+b0(nd8t,m)*xpt*zmt) c /dyt out = beta*zt*lam*xt**theta*n0t**thet1*ht h1t = 1.d0-ht term1 = gamma*c1t**cwgt1*h1t**hwgt*(1.d0-army(m)) term2 = beta*delt1+theta*out/xt term3 = 2.d0*alpha*beta*(n0t-yt) dc1tca = dc1tx*dxca dn0tca = dn0tx*dxca dc1tna = dc1tx*dxna+dc1ty dn0tna = dn0tx*dxna+dn0ty dhtca =-h1t*(dc1tca/c1t+theta*dn0tca/n0t-theta* c dxca/xt) dhtna =-h1t*(dc1tna/c1t+theta*dn0tna/n0t-theta* c dxna/xt) dhtcb =-h1t/c1t dhtnb =-theta*h1t/n0t dt1ca = term1*(cwgt1*dc1tca/c1t-hwgt*dhtca/h1t) dt1cb = term1*(cwgt1/c1t-hwgt*dhtcb/h1t) dt1na = term1*(cwgt1*dc1tna/c1t-hwgt*dhtna/h1t) dt1nb =-term1*hwgt*dhtnb/h1t dt2ca = theta*out/xt*(-thet1*dxca/xt+thet1* c dn0tca/n0t+dhtca/ht) dt2cb = theta*out/xt*dhtcb/ht dt2na = theta*out/xt*(-thet1*dxna/xt+thet1* c dn0tna/n0t+dhtna/ht) dt2nb = theta*out/xt*(thet1/n0t+dhtnb/ht) dt3ca = 2.d0*alpha*beta*dn0tca dt3cb = 0.d0 dt3na = 2.d0*alpha*beta*(dn0tna-1.d0) dt3nb = 2.d0*alpha*beta sum(1) = sum(1) + mcp(l,m)*wgt(mm)*term1*term2 sum(2) = sum(2) + mcp(l,m)*wgt(mm)*term1*term3 drhsc = mcp(l,m)*wgt(mm)*(dt1ca*term2+ c term1*dt2ca) drhsn = mcp(l,m)*wgt(mm)*(dt1na*term2+ c term1*dt2na) c rhs of res1,res2 c d(res1) wrt a0(nd1),b0(nd1), .... sum(3) = sum(3) + drhsc*xp*yp*zp sum(4) = sum(4) + drhsn*xp*yp*zp sum(5) = sum(5) + drhsc*xm*yp*zp sum(6) = sum(6) + drhsn*xm*yp*zp sum(7) = sum(7) + drhsc*xm*ym*zp sum(8) = sum(8) + drhsn*xm*ym*zp sum(9) = sum(9) + drhsc*xp*ym*zp sum(10)= sum(10) + drhsn*xp*ym*zp sum(11)= sum(11) + drhsc*xp*yp*zm sum(12)= sum(12) + drhsn*xp*yp*zm sum(13)= sum(13) + drhsc*xm*yp*zm sum(14)= sum(14) + drhsn*xm*yp*zm sum(15)= sum(15) + drhsc*xm*ym*zm sum(16)= sum(16) + drhsn*xm*ym*zm sum(17)= sum(17) + drhsc*xp*ym*zm sum(18)= sum(18) + drhsn*xp*ym*zm drhsc = mcp(l,m)*wgt(mm)*(dt1ca*term3+ c term1*dt3ca) drhsn = mcp(l,m)*wgt(mm)*(dt1na*term3+ c term1*dt3na) c d(res2) wrt a0(nd1),b0(nd1), .... sum(19)= sum(19) + drhsc*xp*yp*zp sum(20)= sum(20) + drhsn*xp*yp*zp sum(21)= sum(21) + drhsc*xm*yp*zp sum(22)= sum(22) + drhsn*xm*yp*zp sum(23)= sum(23) + drhsc*xm*ym*zp sum(24)= sum(24) + drhsn*xm*ym*zp sum(25)= sum(25) + drhsc*xp*ym*zp sum(26)= sum(26) + drhsn*xp*ym*zp sum(27)= sum(27) + drhsc*xp*yp*zm sum(28)= sum(28) + drhsn*xp*yp*zm sum(29)= sum(29) + drhsc*xm*yp*zm sum(30)= sum(30) + drhsn*xm*yp*zm sum(31)= sum(31) + drhsc*xm*ym*zm sum(32)= sum(32) + drhsn*xm*ym*zm sum(33)= sum(33) + drhsc*xp*ym*zm sum(34)= sum(34) + drhsn*xp*ym*zm c drhsc = mcp(l,m)*wgt(mm)*(dt1cb*term2+ c term1*dt2cb) drhsn = mcp(l,m)*wgt(mm)*(dt1nb*term2+ c term1*dt2nb) c d(res1) wrt a0(nd1t),b0(nd1t), ... wrk(nct,1) = drhsc*xpt*ypt*zpt wrk(nct,2) = drhsn*xpt*ypt*zpt wrk(nct,3) = drhsc*xmt*ypt*zpt wrk(nct,4) = drhsn*xmt*ypt*zpt wrk(nct,5) = drhsc*xmt*ymt*zpt wrk(nct,6) = drhsn*xmt*ymt*zpt wrk(nct,7) = drhsc*xpt*ymt*zpt wrk(nct,8) = drhsn*xpt*ymt*zpt wrk(nct,9) = drhsc*xpt*ypt*zmt wrk(nct,10)= drhsn*xpt*ypt*zmt wrk(nct,11)= drhsc*xmt*ypt*zmt wrk(nct,12)= drhsn*xmt*ypt*zmt wrk(nct,13)= drhsc*xmt*ymt*zmt wrk(nct,14)= drhsn*xmt*ymt*zmt wrk(nct,15)= drhsc*xpt*ymt*zmt wrk(nct,16)= drhsn*xpt*ymt*zmt drhsc = mcp(l,m)*wgt(mm)*(dt1cb*term3+ c term1*dt3cb) drhsn = mcp(l,m)*wgt(mm)*(dt1nb*term3+ c term1*dt3nb) c d(res2) wrt a0(nd1t),b0(nd1t), ... wrk(nct,17)= drhsc*xpt*ypt*zpt wrk(nct,18)= drhsn*xpt*ypt*zpt wrk(nct,19)= drhsc*xmt*ypt*zpt wrk(nct,20)= drhsn*xmt*ypt*zpt wrk(nct,21)= drhsc*xmt*ymt*zpt wrk(nct,22)= drhsn*xmt*ymt*zpt wrk(nct,23)= drhsc*xpt*ymt*zpt wrk(nct,24)= drhsn*xpt*ymt*zpt wrk(nct,25)= drhsc*xpt*ypt*zmt wrk(nct,26)= drhsn*xpt*ypt*zmt wrk(nct,27)= drhsc*xmt*ypt*zmt wrk(nct,28)= drhsn*xmt*ypt*zmt wrk(nct,29)= drhsc*xmt*ymt*zmt wrk(nct,30)= drhsn*xmt*ymt*zmt wrk(nct,31)= drhsc*xpt*ymt*zmt wrk(nct,32)= drhsn*xpt*ymt*zmt 310 continue m = (l-1)*2*na i1 = m+2*nd1-1 i2 = i1+1 i3 = m+2*nd2-1 i4 = i3+1 i5 = m+2*nd3-1 i6 = i5+1 i7 = m+2*nd4-1 i8 = i7+1 i9 = m+2*nd5-1 i10 = i9+1 i11 = m+2*nd6-1 i12 = i11+1 i13 = m+2*nd7-1 i14 = i13+1 i15 = m+2*nd8-1 i16 = i15+1 dlhsc = -mu*army1*(cwgt1/c1-hwgt*dhca/h1) dlhsn = mu*army1*hwgt*dhna/h1 work(i1) = sum(3)+dlhsc*xp*yp*zp work(i2) = sum(4)+dlhsn*xp*yp*zp work(i3) = sum(5)+dlhsc*xm*yp*zp work(i4) = sum(6)+dlhsn*xm*yp*zp work(i5) = sum(7)+dlhsc*xm*ym*zp work(i6) = sum(8)+dlhsn*xm*ym*zp work(i7) = sum(9)+dlhsc*xp*ym*zp work(i8) = sum(10)+dlhsn*xp*ym*zp work(i9) = sum(11)+dlhsc*xp*yp*zm work(i10) = sum(12)+dlhsn*xp*yp*zm work(i11) = sum(13)+dlhsc*xm*yp*zm work(i12) = sum(14)+dlhsn*xm*yp*zm work(i13) = sum(15)+dlhsc*xm*ym*zm work(i14) = sum(16)+dlhsn*xm*ym*zm work(i15) = sum(17)+dlhsc*xp*ym*zm work(i16) = sum(18)+dlhsn*xp*ym*zm tem1 = c1-c0-z(k)*lam*(x(i)/n0)**theta*thet1*h+ c 2.d0*alpha*(n0-y(j)) dlhsc =-army1*(mu-gamm1/gamma*mu*c1/h1*dhca- c gamma*c0**cwgt1*dc0ca-tem1*(cwgt1*mu/c1- c hwgt*mu/h1*dhca)-mu*(1.d0-dc0ca-z(k)* c lam*(x(i)/n0)**theta*thet1*dhca)) dlhsn =-army1*(-gamm1/gamma*mu*c1/h1*dhna-gamma* c c0**cwgt1*dc0na-d2pn+tem1*hwgt*mu/h1* c dhna-mu*(-dc0na+z(k)*lam*(x(i)/n0)** c theta*thet1*h*(theta/n0-dhna/h)+2.d0* c alpha)) wrk2(i1) = sum(19)+dlhsc*xp*yp*zp wrk2(i2) = sum(20)+dlhsn*xp*yp*zp wrk2(i3) = sum(21)+dlhsc*xm*yp*zp wrk2(i4) = sum(22)+dlhsn*xm*yp*zp wrk2(i5) = sum(23)+dlhsc*xm*ym*zp wrk2(i6) = sum(24)+dlhsn*xm*ym*zp wrk2(i7) = sum(25)+dlhsc*xp*ym*zp wrk2(i8) = sum(26)+dlhsn*xp*ym*zp wrk2(i9) = sum(27)+dlhsc*xp*yp*zm wrk2(i10) = sum(28)+dlhsn*xp*yp*zm wrk2(i11) = sum(29)+dlhsc*xm*yp*zm wrk2(i12) = sum(30)+dlhsn*xm*yp*zm wrk2(i13) = sum(31)+dlhsc*xm*ym*zm wrk2(i14) = sum(32)+dlhsn*xm*ym*zm wrk2(i15) = sum(33)+dlhsc*xp*ym*zm wrk2(i16) = sum(34)+dlhsn*xp*ym*zm do 330 mm=1,ne1 do 330 m=1,ns nct = (m-1)*ne1+mm j1 = (m-1)*2*na+2*ind(nct)-1 j2 = j1+1 j3 = j1+2 j4 = j3+1 j5 = j3+2*nx1 j6 = j5+1 j7 = j1+2*nx1 j8 = j7+1 j9 = j1+2*nx1*ny1 j10 = j9+1 j11 = j3+2*nx1*ny1 j12 = j11+1 j13 = j11+2*nx1 j14 = j13+1 j15 = j9+2*nx1 j16 = j15+1 work(j1) = work(j1) + wrk(nct,1) work(j2) = work(j2) + wrk(nct,2) work(j3) = work(j3) + wrk(nct,3) work(j4) = work(j4) + wrk(nct,4) work(j5) = work(j5) + wrk(nct,5) work(j6) = work(j6) + wrk(nct,6) work(j7) = work(j7) + wrk(nct,7) work(j8) = work(j8) + wrk(nct,8) work(j9) = work(j9) + wrk(nct,9) work(j10) = work(j10) + wrk(nct,10) work(j11) = work(j11) + wrk(nct,11) work(j12) = work(j12) + wrk(nct,12) work(j13) = work(j13) + wrk(nct,13) work(j14) = work(j14) + wrk(nct,14) work(j15) = work(j15) + wrk(nct,15) work(j16) = work(j16) + wrk(nct,16) wrk2(j1) = wrk2(j1) + wrk(nct,17) wrk2(j2) = wrk2(j2) + wrk(nct,18) wrk2(j3) = wrk2(j3) + wrk(nct,19) wrk2(j4) = wrk2(j4) + wrk(nct,20) wrk2(j5) = wrk2(j5) + wrk(nct,21) wrk2(j6) = wrk2(j6) + wrk(nct,22) wrk2(j7) = wrk2(j7) + wrk(nct,23) wrk2(j8) = wrk2(j8) + wrk(nct,24) wrk2(j9) = wrk2(j9) + wrk(nct,25) wrk2(j10) = wrk2(j10) + wrk(nct,26) wrk2(j11) = wrk2(j11) + wrk(nct,27) wrk2(j12) = wrk2(j12) + wrk(nct,28) wrk2(j13) = wrk2(j13) + wrk(nct,29) wrk2(j14) = wrk2(j14) + wrk(nct,30) wrk2(j15) = wrk2(j15) + wrk(nct,31) wrk2(j16) = wrk2(j16) + wrk(nct,32) 330 continue res1 = (-mu*army1+sum(1))*u(i)*v(j)*w(k) res2 = ((-c1**cwgt*h1**hwgt/ome1+c0**cwgt/ c ome1+dpn+mu*(c1-c0-z(k)*lam*(x(i)/ c n0)**theta*thet1*h+2.d0*alpha*(n0- c y(j))))*army1+sum(2))*u(i)*v(j)*w(k) f(i1) = f(i1) + res1*xp*yp*zp f(i2) = f(i2) + res2*xp*yp*zp f(i3) = f(i3) + res1*xm*yp*zp f(i4) = f(i4) + res2*xm*yp*zp f(i5) = f(i5) + res1*xm*ym*zp f(i6) = f(i6) + res2*xm*ym*zp f(i7) = f(i7) + res1*xp*ym*zp f(i8) = f(i8) + res2*xp*ym*zp f(i9) = f(i9) + res1*xp*yp*zm f(i10) = f(i10) + res2*xp*yp*zm f(i11) = f(i11) + res1*xm*yp*zm f(i12) = f(i12) + res2*xm*yp*zm f(i13) = f(i13) + res1*xm*ym*zm f(i14) = f(i14) + res2*xm*ym*zm f(i15) = f(i15) + res1*xp*ym*zm f(i16) = f(i16) + res2*xp*ym*zm c c Store only nonzero elements of df: c tem1 = u(i)*v(j)*w(k) m1 = 0 do 340 mm=1,nas work(mm) = tem1*work(mm) if (work(mm).ge.tiny.or.work(mm).le.-tiny) then m1 = m1+1 iwork(m1)= mm endif 340 continue do 345 mm=1,m1 iw = iwork(mm) nfound(1) = 0 nfound(2) = 0 nfound(3) = 0 nfound(4) = 0 nfound(5) = 0 nfound(6) = 0 nfound(7) = 0 nfound(8) = 0 do 341 m2=1,ndf if (iw.eq.ja(m2)) then if (i1.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xp*yp*zp nfound(1) = 1 endif if (i3.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xm*yp*zp nfound(2) = 1 endif if (i5.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xm*ym*zp nfound(3) = 1 endif if (i7.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xp*ym*zp nfound(4) = 1 endif if (i9.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xp*yp*zm nfound(5) = 1 endif if (i11.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xm*yp*zm nfound(6) = 1 endif if (i13.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xm*ym*zm nfound(7) = 1 endif if (i15.eq.ia(m2)) then df(m2) = df(m2) + work(iw)*xp*ym*zm nfound(8) = 1 endif endif 341 continue if (nfound(1).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xp*yp*zp ja(ndf) = iw ia(ndf) = i1 endif if (nfound(2).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xm*yp*zp ja(ndf) = iw ia(ndf) = i3 endif if (nfound(3).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xm*ym*zp ja(ndf) = iw ia(ndf) = i5 endif if (nfound(4).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xp*ym*zp ja(ndf) = iw ia(ndf) = i7 endif if (nfound(5).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xp*yp*zm ja(ndf) = iw ia(ndf) = i9 endif if (nfound(6).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xm*yp*zm ja(ndf) = iw ia(ndf) = i11 endif if (nfound(7).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xm*ym*zm ja(ndf) = iw ia(ndf) = i13 endif if (nfound(8).eq.0) then ndf = ndf + 1 df(ndf) = work(iw)*xp*ym*zm ja(ndf) = iw ia(ndf) = i15 endif 345 continue m1 = 0 do 342 mm=1,nas wrk2(mm) = tem1*wrk2(mm) if (wrk2(mm).ge.tiny.or.wrk2(mm).le.-tiny) then m1 = m1+1 iwork(m1)= mm endif 342 continue do 347 mm=1,m1 iw = iwork(mm) nfound(1) = 0 nfound(2) = 0 nfound(3) = 0 nfound(4) = 0 nfound(5) = 0 nfound(6) = 0 nfound(7) = 0 nfound(8) = 0 do 343 m2=1,ndf if (iw.eq.ja(m2)) then if (i2.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xp*yp*zp nfound(1) = 1 endif if (i4.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xm*yp*zp nfound(2) = 1 endif if (i6.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xm*ym*zp nfound(3) = 1 endif if (i8.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xp*ym*zp nfound(4) = 1 endif if (i10.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xp*yp*zm nfound(5) = 1 endif if (i12.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xm*yp*zm nfound(6) = 1 endif if (i14.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xm*ym*zm nfound(7) = 1 endif if (i16.eq.ia(m2)) then df(m2) = df(m2) + wrk2(iw)*xp*ym*zm nfound(8) = 1 endif endif 343 continue if (nfound(1).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xp*yp*zp ja(ndf) = iw ia(ndf) = i2 endif if (nfound(2).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xm*yp*zp ja(ndf) = iw ia(ndf) = i4 endif if (nfound(3).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xm*ym*zp ja(ndf) = iw ia(ndf) = i6 endif if (nfound(4).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xp*ym*zp ja(ndf) = iw ia(ndf) = i8 endif if (nfound(5).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xp*yp*zm ja(ndf) = iw ia(ndf) = i10 endif if (nfound(6).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xm*yp*zm ja(ndf) = iw ia(ndf) = i12 endif if (nfound(7).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xm*ym*zm ja(ndf) = iw ia(ndf) = i14 endif if (nfound(8).eq.0) then ndf = ndf + 1 df(ndf) = wrk2(iw)*xp*ym*zm ja(ndf) = iw ia(ndf) = i16 endif 347 continue write(*,*) it,ndf 240 continue 150 continue 140 continue itask = 1 write(*,*) 'before arnoldi' call arnoldi(df, ia, ja, ldf, ndf, f, nas, work, iwork) write(*,*) 'after arnoldi' sum3 = 0.d0 do 350 j=1,ns do 350 i=1,na i1 = (j-1)*na*2+2*i-1 tem1 = f(i1) tem2 = f(i1+1) a1(i,j) = a0(i,j) - tem1 b1(i,j) = b0(i,j) - tem2 sum3 = sum3 + tem1*tem1+tem2*tem2 350 continue dd = dsqrt(sum3)/float(nas) write(*,*) ' at iteration ', it, ' the residual is ', dd if (dd.lt.crit) go to 999 do 360 j=1,ns do 360 i=1,na a0(i,j) = a1(i,j) b0(i,j) = b1(i,j) 360 continue 120 continue 999 continue c c print results c write(7,*) 'CAPACITY UTILIZATION-MOVING COSTS-3D-B RESULTS' write(7,*) write(7,*) ' Parameters ' write(7,'('' alpha = '',E9.3,'' '')') alpha write(7,'('' beta = '',E9.3,'' '')') beta write(7,'('' delta = '',E9.3,'' '')') delta write(7,'('' gamma = '',E9.3,'' '')') gamma write(7,'('' lam = '',E9.3,'' '')') lam write(7,'('' omega = '',E9.3,'' '')') omega write(7,'('' pa = '',E9.3,'' '')') pa write(7,'('' pb = '',E9.3,'' '')') pb write(7,'('' phi = '',E9.3,'' '')') phi write(7,'('' psi = '',E9.3,'' '')') psi write(7,'('' rho = '',E9.3,'' '')') rho write(7,'('' sig = '',E9.3,'' '')') sig write(7,'('' theta = '',E9.3,'' '')') theta write(7,'('' zeta = '',E9.3,'' '')') zeta write(7,*) write(7,*) ' Markov Chain ' write(7,*) ' transition matrix: ' if (ns.eq.1) then write(7,801) mcp(1,1) else do 370 i=1,ns if (ns.eq.2) write(7,802) (mcp(i,j),j=1,ns) if (ns.eq.3) write(7,803) (mcp(i,j),j=1,ns) if (ns.eq.4) write(7,804) (mcp(i,j),j=1,ns) if (ns.eq.5) write(7,805) (mcp(i,j),j=1,ns) if (ns.eq.6) write(7,806) (mcp(i,j),j=1,ns) if (ns.eq.7) write(7,807) (mcp(i,j),j=1,ns) if (ns.eq.8) write(7,808) (mcp(i,j),j=1,ns) 370 continue endif write(7,*) write(7,*) ' values of ginv, cg, army in each state: ' do 380 i=1,ns write(7,803) ginv(i),cg(i),army(i) 380 continue write(7,*) write(7,*) ' Integrals: ' write(7,'('' sum1 = '',E9.3,'' '')') sum1 write(7,'('' sum2 = '',E9.3,'' '')') sum2 write(7,*) c write(7,*) write(7,*) ' Coefficients for Shape functions (Consumption)' do 390 i=1,na if (ns.eq.1) write(7,801) (a0(i,j),j=1,ns) if (ns.eq.2) write(7,802) (a0(i,j),j=1,ns) if (ns.eq.3) write(7,803) (a0(i,j),j=1,ns) if (ns.eq.4) write(7,804) (a0(i,j),j=1,ns) if (ns.eq.5) write(7,805) (a0(i,j),j=1,ns) if (ns.eq.6) write(7,806) (a0(i,j),j=1,ns) if (ns.eq.7) write(7,807) (a0(i,j),j=1,ns) if (ns.eq.8) write(7,808) (a0(i,j),j=1,ns) 390 continue write(7,*) ' Coefficients for Shape functions (Employment)' do 395 i=1,na if (ns.eq.1) write(7,801) (b0(i,j),j=1,ns) if (ns.eq.2) write(7,802) (b0(i,j),j=1,ns) if (ns.eq.3) write(7,803) (b0(i,j),j=1,ns) if (ns.eq.4) write(7,804) (b0(i,j),j=1,ns) if (ns.eq.5) write(7,805) (b0(i,j),j=1,ns) if (ns.eq.6) write(7,806) (b0(i,j),j=1,ns) if (ns.eq.7) write(7,807) (b0(i,j),j=1,ns) if (ns.eq.8) write(7,808) (b0(i,j),j=1,ns) 395 continue write(7,*) write(7,*) ' Grid Points (xa,ya,za) ' write(7,*) ' xa ya za ' write(7,*) ' ------------------------------ ' do 400 k=1,nz1 do 401 j=1,ny1 do 402 i=1,nx1 write(7,813) xa(i),ya(j),za(k) 402 continue 401 continue 400 continue write(7,*) write(7,*) ' Number of iterations = ', it c c write file to be used as next input file: c write(9,807) alpha,beta,delta,gamma,lam,omega,pa write(9,807) pb,phi,psi,rho,sig,theta,zeta do 430 i=1,ns write(9,803) ginv(i),cg(i),army(i) 430 continue do 440 i=1,ns write(9,*) (mcp(i,j),j=1,ns) 440 continue do 450 i=1,na write(9,*) (a0(i,j),j=1,ns) 450 continue do 455 i=1,na write(9,*) (b0(i,j),j=1,ns) 455 continue do 460 i=1,nx write(9,811) xa(i),nxp(i) 460 continue write(9,801) xa(nx1) do 470 i=1,ny write(9,811) ya(i),nyp(i) 470 continue write(9,801) ya(ny1) do 480 i=1,nz write(9,811) za(i),nzp(i) 480 continue write(9,801) za(nz1) c c write file to be read into matlab: c write(10,801) alpha write(10,801) beta write(10,801) delta write(10,801) gamma write(10,801) lam write(10,801) omega write(10,801) pa write(10,801) pb write(10,801) phi write(10,801) psi write(10,801) rho write(10,801) sig write(10,801) theta write(10,801) zeta write(10,812) ns do 490 i=1,ns write(10,801) ginv(i) write(10,801) cg(i) write(10,801) army(i) 490 continue do 500 j=1,ns do 500 i=1,ns write(10,801) mcp(i,j) 500 continue write(10,812) nx1 write(10,812) ny1 write(10,812) nz1 do 510 i=1,nx1 write(10,801) xa(i) 510 continue do 520 i=1,ny1 write(10,801) ya(i) 520 continue do 530 i=1,nz1 write(10,801) za(i) 530 continue do 540 j=1,ns do 540 i=1,na write(10,801) a0(i,j) 540 continue do 550 j=1,ns do 550 i=1,na write(10,801) b0(i,j) 550 continue c 800 format (1X, I4, F12.6) 801 format (4X, E12.6) 802 format (2(4X, E12.5)) 803 format (3(4X, E12.5)) 804 format (4(4X, E12.5)) 805 format (5(4X, E12.5)) 806 format (6(2X, E12.5)) 807 format (7(2X, E12.5)) 808 format (8(2X, E12.5)) 811 format (4X, E12.6, 4X, I4) 812 format (4X, I4) 813 format (4X, 3(1X, E9.3)) stop end