CFPP$ UNROLL (32) F c c TITLE: implementation of the finite-element method for model c two, with capacity utilization, in ``The Macroeconomics c of War and Peace,'' NBER Macroeconomics Annual. c c PROBLEM: find c1^l(k,kg,z), h^l(k,kg,z), c0^l(k,kg,z), c and n^l(k,kg,z) that solve the functional equation: c u_c(c^l,1-h^l)(1-a^l) - beta sum_{j=1:n} pi(l,j) c * int u_c(c^l(x~),1-h^l(x~))(1-a^j) c *(1-delta+z~*F_k(k~,kg~,n^j(x~),h^j(x~))) c *f(eps) d(eps), where x~=(k~,kg~,z~) c for l=1,..,nstates c k~ = z*F(k,kg,n^l(k,kg,z),h^l(k,kg,z)) c -n^l(k,kg,z)*c1^l(k,kg,z) c -(1-n^l(k,kg,z))*c0^l(k,kg,z) c +(1-delta)*k-ig^l-(cg^l+zeta*log(z)) c kg~ = (1-delta)*kg+ig^l c z~ = z^psi*eps c c1(0,0,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 c MODEL: stochastic growth model with government capital c and varying hours and employment; the planner c chooses c1[t],c0[t],i[t],n[t],h[t] 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)}*(1-a[t]) | x[0]] c c subject to c c c[t]+i[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 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],kg[t],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 c TEST CASE: need parameterization to get n[t]=1 c c INPUTS: (a) cu3da.inp: beta,delta,gamma,lam,omega,pa, ... c pb,phi,psi,rho,sig,theta c a0(i,j) i=1:na, j=1:n 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 a0's are initial guesses for policy rules 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) cu3da.out: general results c (c) cu3da.nxt: new input file c (d) cu3d.dat: input file for cu3d.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: cu3db.f (which solves Ax=b differently) c c Ellen McGrattan, 4-12-93 c Revised, 7-18-93, ERM c parameter (na=96,nx=15,ny=1,nz=2,ns=5,nas=480,nxm=3,nym=3, c nzm=3,ne1=10,nt=10000,nw=10000) implicit double precision (a-h,o-z) double precision f(nas), df(nas,nas), 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 c1(nt,ns), b(nt,ns), res(nt,ns), n0(nt,ns), h(nt,ns), c a0(na,ns), a1(na,ns), zzt(ne1), e1(ne1), w1(ne1), c wgt(ne1), sum(9), wrk(ne1*ns,8), mcp(ns,ns), ginv(ns), c cg(ns), army(ns), lam, n0t integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), e(nx*ny*nz,8), c kk(ne1), ind(ne1*ns) c tiny = 0.000000001d0 pi = 3.14159265d0 crit = 0.000001d0 maxit = 20 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 open(unit=5, file='cu3da.inp') open(unit=7, file='cu3da.out') open(unit=9, file='cu3da.nxt') open(unit=10,file='cu3d.dat') read(5,*) beta,delta,gamma,lam,omega,pa,pb,phi,psi,rho,sig, c 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 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 n = 0 do 80 k=1,nz do 80 j=1,ny do 80 i=1,nx ll = (k-1)*nx*ny+(j-1)*nx+i e(ll,1) = (k-1)*nx1*ny1+(j-1)*nx1+i e(ll,2) = i e(ll,3) = j e(ll,4) = k e(ll,5) = nxp(i) e(ll,6) = nyp(j) e(ll,7) = nzp(k) e(ll,8) = n n = n+nxp(i)*nyp(j)*nzp(k) 80 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 dd = 1.d0 do 120 it=1,maxit do 130 i=1,nas f(i) = 0.d0 do 130 j=1,nas df(i,j) = 0.d0 130 continue do 140 n=1,ne n1 = e(n,5) n2 = e(n,6) n3 = e(n,7) i1 = e(n,8) x1 = xa(e(n,2)) x2 = xa(e(n,2)+1) dx = x2-x1 y1 = ya(e(n,3)) y2 = ya(e(n,3)+1) dy = y2-y1 z1 = za(e(n,4)) z2 = za(e(n,4)+1) dz = z2-z1 nd1= e(n,1) nd2= nd1+1 nd4= nd1+nx1 nd3= nd4+1 nd5= nd1+nx1*ny1 nd6= nd5+1 nd8= nd5+nx1 nd7= nd8+1 do 150 j=1,n2 y(j) = yy(j,e(n,3)) 150 continue do 160 k=1,n3 z(k) = zz(k,e(n,4)) 160 continue c create vectors a,b,c for consumption, leisure, next capital do 170 i=1,n1 x(i) = xx(i,e(n,2)) xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx do 170 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 170 k=1,n3 zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz ll = i1+(k-1)*n1*n2+(j-1)*n1+i do 170 l=1,ns c1(ll,l) = 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 c1(ll,l) = dmax1(c1(ll,l),tiny) b(ll,l) = z(k)*lam*(phi*x(i)**rho+(1.d0-phi)* c y(j)**rho)**(theta/rho) 170 continue nn = n1*n2*n3 do 180 ll=i1+1,i1+nn do 180 l=1,ns hrs = 0.4d0 emp = 0.2d0 do 190 m=1,5 c0 = c1(ll,l)*(1.d0-hrs)**(hwgt/cwgt1) c pn = abs((emp*pb-pa*pb)**3.0)+(emp*pb-pa*pb)**3.0 d = (emp-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb fj11 = -d2pn/(c1(ll,l)**cwgt)-gamma*(1.d0-hrs)** c hwgt/c1(ll,l)*b(ll,l)*thet1*theta*emp** c (-theta-1.d0)*hrs fj22 = b(ll,l) fj12 = hwgt*(1.d0-hrs)**hwgt1/ome1-gamma*hwgt* c (1.d0-hrs)**hwgt1/c1(ll,l)*(c1(ll,l)-c0- c b(ll,l)*thet1*emp**(-theta)*hrs)-gamma* c (1.d0-hrs)**hwgt/c1(ll,l)*b(ll,l)*thet1* c emp**(-theta) fj21 = -gamm1*theta*emp**(-thet1)*c1(ll,l)/gamma det = fj11*fj22-fj12*fj21 f1 = gamm1/gamma*emp**theta*c1(ll,l)-b(ll,l)* c (1.d0-hrs) f2 = (1.d0-hrs)**hwgt/ome1-(c0/c1(ll,l))**cwgt/ c ome1-dpn/(c1(ll,l)**cwgt)-gamma*(1.d0-hrs)** c hwgt/c1(ll,l)*(c1(ll,l)-c0-b(ll,l)*thet1* c emp**(-theta)*hrs) hrs = hrs - (fj11*f1+fj21*f2)/det emp = emp - (fj12*f1+fj22*f2)/det hrs = dmax1(hrs,tiny) emp = dmax1(emp,tiny) 190 continue h(ll,l) = dmin1(hrs,1.d0) n0(ll,l) = dmin1(emp,1.d0) 180 continue 140 continue do 200 l=1,ns do 210 n=1,ne n1 = e(n,5) n2 = e(n,6) n3 = e(n,7) i1 = e(n,8) x1 = xa(e(n,2)) x2 = xa(e(n,2)+1) dx = x2-x1 y1 = ya(e(n,3)) y2 = ya(e(n,3)+1) dy = y2-y1 z1 = za(e(n,4)) z2 = za(e(n,4)+1) dz = z2-z1 nd1= e(n,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,e(n,2)) u(i) = wx(i,e(n,2)) 220 continue do 230 j=1,n2 y(j) = yy(j,e(n,3)) v(j) = wy(j,e(n,3)) 230 continue do 240 k=1,n3 z(k) = zz(k,e(n,4)) w(k) = wz(k,e(n,4)) 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 ll = (k-1)*n1*n2+(j-1)*n1+i+i1 xp = (x2 - x(i))/dx xm = (x(i) - x1)/dx c0 = c1(ll,l)*(1.d0-h(ll,l))**(hwgt/cwgt1) cg1 = cg(l)+zeta*dlog(z(k)) xt = b(ll,l)*n0(ll,l)**thet1*h(ll,l)+delt1* c x(i)-ginv(l)-cg1-n0(ll,l)*c1(ll,l)- c (1.d0-n0(ll,l))*c0 yt = (1-delta)*y(j)+ginv(l) ii = 1 do 250 m=1,nx j1 = int( 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 = int( 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 = int( 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,9 sum(m) = 0.d0 290 continue do 300 m=1,nas work(m) = 0.d0 300 continue xy = phi*xt**rho+(1.d0-phi)*yt**rho xr = xt**(rho-1.d0) h1 = 1.d0-h(ll,l) crat = c0/c1(ll,l) d = (n0(ll,l)-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb tem1 = gamm1/h1*(crat**cwgt-thet1*h1**hwgt1+ c ome1*dpn/(c1(ll,l)**cwgt)-d2pn*n0(ll,l)/ c (gamm1*theta*c1(ll,l)**cwgt)) tem2 = (d2pn*n0(ll,l)/theta+dpn*cwgt1)/ c (c1(ll,l)**(cwgt+1.d0)) dha = tem2/tem1 tem1 = c1(ll,l)/h1*(-gamm1*thet1*h(ll,l)*n0(ll,l)/ c (gamma*theta*h1)+n0(ll,l)/theta*(1.d0-crat)+ c gamm1*n0(ll,l)/gamma+crat*(1.d0-n0(ll,l))* c hwgt/cwgt1) tem2 = -gamm1*thet1*h(ll,l)*n0(ll,l)/(gamma*theta* c h1)+thet1*n0(ll,l)*(1.d0-crat)/theta-crat dxa = tem1*dha+tem2 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(e(nn,2)) x2t = xa(e(nn,2)+1) dxt = x2t-x1t y1t = ya(e(nn,3)) y2t = ya(e(nn,3)+1) dyt = y2t-y1t z1t = za(e(nn,4)) z2t = za(e(nn,4)+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 = e(nn,1) 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 c1t = dmax1(c1t,tiny) bt = zt*lam*(phi*xt**rho+(1.d0-phi)* c yt**rho)**(theta/rho) hrs = 0.4d0 emp = 0.2d0 do 320 mt=1,5 c0t = c1t*(1.d0-hrs)**(hwgt/cwgt1) d = (emp-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb d2pn = 6.d0*(d+dsign(d,0.d0))*pb*pb fj11 = -d2pn/(c1t**cwgt)-gamma*(1.d0-hrs)** c hwgt/c1t*bt*thet1*theta*emp**(-theta- c 1.d0)*hrs fj22 = bt dc0 = -c1t*(hwgt/cwgt1)*(1.d0-hrs)**(hwgt/ c cwgt1-1.d0) fj12 = hwgt*(1.d0-hrs)**hwgt1/ome1-gamma*hwgt* c (1.d0-hrs)**hwgt1/c1t*(c1t-c0t-bt* c thet1*emp**(-theta)*hrs)-gamma*(1.d0- c hrs)**hwgt/c1t*bt*thet1*emp**(-theta) fj21 = -gamm1*theta*emp**(-thet1)*c1t/gamma det = fj11*fj22-fj12*fj21 f1 = gamm1/gamma*emp**theta*c1t-bt* c (1.d0-hrs) f2 = (1.d0-hrs)**hwgt/ome1-(c0t/c1t)**cwgt/ c ome1-dpn/(c1t**cwgt)-gamma*(1.d0- c hrs)**hwgt/c1t*(c1t-c0t-bt*thet1*emp** c (-theta)*hrs) hrs = hrs - (fj11*f1+fj21*f2)/det emp = emp - (fj12*f1+fj22*f2)/det hrs = dmax1(hrs,tiny) emp = dmax1(emp,tiny) 320 continue ht = dmin1(hrs,1.d0) n0t = dmin1(emp,1.d0) d = (n0t-pa)*pb dpnt = 3.d0*d*(d+dsign(d,0.d0))*pb d2pnt = 6.d0*(d+dsign(d,0.d0))*pb*pb dc1t = (-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 out = beta*zt*lam*xy**(theta/rho)*n0t**thet1*ht trm1 = beta*delt1+out/xy*phi*theta*xr h1 = 1.d0-ht crat = c0t/c1t tem1 = gamm1/h1*(crat**cwgt-thet1*h1**hwgt1+ c ome1*dpnt/(c1t**cwgt)-d2pnt*n0t/ c (gamm1*theta*c1t**cwgt)) tem2 = (d2pnt*n0t/theta+dpnt*cwgt1)/(c1t** c (cwgt+1.d0)) tem3 =-d2pnt*phi*xr*n0t/(xy*c1t**cwgt) dhta = (tem2*dc1t+tem3)*dxa/tem1 dnta =-n0t*dhta/(theta*h1)-n0t*dc1t*dxa/(theta* c c1t)+phi*xr*n0t*dxa/xy dtm1a = phi*theta*xr*out/xy*((-rho1/xt+(theta/ c rho-1.d0)*phi*rho*xr/xy)*dxa+thet1/n0t* c dnta+dhta/ht) dhtb = tem2/tem1 dntb =-n0t*(dhtb/h1+1.d0/c1t)/theta dtm1b = phi*theta*xr*out/xy*(thet1/n0t*dntb+ c dhtb/ht) trm2 = gamma*c1t**cwgt1*h1**hwgt*(1.d0-army(m)) dtm2a = trm2*(cwgt1/c1t*dc1t*dxa-hwgt/h1*dhta) dtm2b = trm2*(cwgt1/c1t-hwgt/h1*dhtb) c sum(1) = sum(1) + mcp(l,m)*wgt(mm)*trm1*trm2 c drhs = mcp(l,m)*wgt(mm)*(dtm1a*trm2+trm1*dtm2a) sum(2) = sum(2) + drhs*xp*yp*zp sum(3) = sum(3) + drhs*xm*yp*zp sum(4) = sum(4) + drhs*xm*ym*zp sum(5) = sum(5) + drhs*xp*ym*zp sum(6) = sum(6) + drhs*xp*yp*zm sum(7) = sum(7) + drhs*xm*yp*zm sum(8) = sum(8) + drhs*xm*ym*zm sum(9) = sum(9) + drhs*xp*ym*zm c drhs = mcp(l,m)*wgt(mm)*(dtm1b*trm2+trm1*dtm2b) wrk(nct,1) = drhs*xpt*ypt*zpt wrk(nct,2) = drhs*xmt*ypt*zpt wrk(nct,3) = drhs*xmt*ymt*zpt wrk(nct,4) = drhs*xpt*ymt*zpt wrk(nct,5) = drhs*xpt*ypt*zmt wrk(nct,6) = drhs*xmt*ypt*zmt wrk(nct,7) = drhs*xmt*ymt*zmt wrk(nct,8) = drhs*xpt*ymt*zmt 310 continue dlhs = -gamma*c1(ll,l)**cwgt1*(1.d0-h(ll,l))**hwgt* c (1.d0-army(l))*(cwgt1/c1(ll,l)-hwgt/ c (1.d0-h(ll,l))*dha) m = (l-1)*na work(m+nd1) = sum(2)+dlhs*xp*yp*zp work(m+nd2) = sum(3)+dlhs*xm*yp*zp work(m+nd3) = sum(4)+dlhs*xm*ym*zp work(m+nd4) = sum(5)+dlhs*xp*ym*zp work(m+nd5) = sum(6)+dlhs*xp*yp*zm work(m+nd6) = sum(7)+dlhs*xm*yp*zm work(m+nd7) = sum(8)+dlhs*xm*ym*zm work(m+nd8) = sum(9)+dlhs*xp*ym*zm do 330 mm=1,ne1 do 330 m=1,ns nct = (m-1)*ne1+mm nd1t = ind(nct) nd2t = nd1t+1 nd4t = nd1t+nx1 nd3t = nd4t+1 nd5t = nd1t+nx1*ny1 nd6t = nd5t+1 nd8t = nd5t+nx1 nd7t = nd8t+1 mt = (m-1)*na work(mt+nd1t) = work(mt+nd1t) + wrk(nct,1) work(mt+nd2t) = work(mt+nd2t) + wrk(nct,2) work(mt+nd3t) = work(mt+nd3t) + wrk(nct,3) work(mt+nd4t) = work(mt+nd4t) + wrk(nct,4) work(mt+nd5t) = work(mt+nd5t) + wrk(nct,5) work(mt+nd6t) = work(mt+nd6t) + wrk(nct,6) work(mt+nd7t) = work(mt+nd7t) + wrk(nct,7) work(mt+nd8t) = work(mt+nd8t) + wrk(nct,8) 330 continue res(ll,l)= (-gamma*c1(ll,l)**cwgt1*(1.d0-h(ll,l))** c hwgt*(1.d0-army(l))+sum(1))*u(i)*v(j)*w(k) m = (l-1)*na f(m+nd1) = f(m+nd1) + res(ll,l)*xp*yp*zp f(m+nd2) = f(m+nd2) + res(ll,l)*xm*yp*zp f(m+nd3) = f(m+nd3) + res(ll,l)*xm*ym*zp f(m+nd4) = f(m+nd4) + res(ll,l)*xp*ym*zp f(m+nd5) = f(m+nd5) + res(ll,l)*xp*yp*zm f(m+nd6) = f(m+nd6) + res(ll,l)*xm*yp*zm f(m+nd7) = f(m+nd7) + res(ll,l)*xm*ym*zm f(m+nd8) = f(m+nd8) + res(ll,l)*xp*ym*zm c tem1 = u(i)*v(j)*w(k) do 340 mm=1,nas tem2 = tem1*work(mm) df(m+nd1,mm) = df(m+nd1,mm) + tem2*xp*yp*zp df(m+nd2,mm) = df(m+nd2,mm) + tem2*xm*yp*zp df(m+nd3,mm) = df(m+nd3,mm) + tem2*xm*ym*zp df(m+nd4,mm) = df(m+nd4,mm) + tem2*xp*ym*zp df(m+nd5,mm) = df(m+nd5,mm) + tem2*xp*yp*zm df(m+nd6,mm) = df(m+nd6,mm) + tem2*xm*yp*zm df(m+nd7,mm) = df(m+nd7,mm) + tem2*xm*ym*zm df(m+nd8,mm) = df(m+nd8,mm) + tem2*xp*ym*zm 340 continue 240 continue 210 continue 200 continue c do 344 j=1,nas c write(75,801) f(j) c do 344 i=1,nas c write(76,801) df(i,j) c344 continue c stop itask = 1 call dgefs (df, nas, nas, f, itask, ind1, work, iwork, rcond) sum3 = 0.d0 do 350 j=1,ns do 350 i=1,na tem1 = f((j-1)*na+i) a1(i,j) = a0(i,j) - tem1 sum3 = sum3 + tem1*tem1 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) 360 continue 120 continue 999 continue c c print results c write(7,*) 'CAPACITY UTILIZATION-3D-A RESULTS' write(7,*) write(7,*) ' Parameters ' write(7,'('' beta = '',F10.5,'' '')') beta write(7,'('' delta = '',F10.5,'' '')') delta write(7,'('' gamma = '',F10.5,'' '')') gamma write(7,'('' lam = '',F10.5,'' '')') lam write(7,'('' omega = '',F10.5,'' '')') omega write(7,'('' pa = '',F10.5,'' '')') pa write(7,'('' pb = '',F10.5,'' '')') pb write(7,'('' phi = '',F10.5,'' '')') phi write(7,'('' psi = '',F10.5,'' '')') psi write(7,'('' rho = '',F10.5,'' '')') rho write(7,'('' sig = '',F10.5,'' '')') sig write(7,'('' theta = '',F10.5,'' '')') theta write(7,'('' zeta = '',F10.5,'' '')') 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 = '',F6.2,'' '')') sum1 write(7,'('' sum2 = '',F6.2,'' '')') sum2 write(7,*) write(7,*) write(7,*) ' Results from dgefs.f ' if (ind1 .eq. -10) then write (7,*) ' error code =', ind1 else if (ind1 .lt. 0) then write (7,*) ' error code =', ind1 stop else write (7,*) ' number of accurate digits =', ind1 end if c write(7,*) write(7,*) ' Coefficients for Shape functions ' 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,*) 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,806) 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,807) (mcp(i,j),j=1,ns) 440 continue do 450 i=1,na write(9,807) (a0(i,j),j=1,ns) 450 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) 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 c 801 format (1X, F12.6) 802 format (2(1X, F12.6)) 803 format (3(1X, F12.6)) 804 format (4(1X, F12.6)) 805 format (5(1X, F12.6)) 806 format (6(1X, F11.5)) 807 format (7(1X, F10.5)) 808 format (8(1X, F9.5)) 811 format (1X, F12.6, 4X, I4) 812 format (9X, I4) 813 format (4X, 3(1X, F12.6)) stop end