CFPP$ UNROLL (32) F c c SUBROUTINE: CUMCMOM.F c c TITLE: moment calculation for simulated method of moments c estimation of parameters in cumc3da.f model c c PROBLEM: calculate equilibria as in cumc3da.f, simulate c time series, calculate certain moments that will c be compared with those in the data c c ALGORITHM: finite element method used to compute equilibria; c simulated method of moments used to estimate parameters 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*(k+kg)^theta*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 PARAMETERS: param = [alpha,beta,delta,gamma,lam,omega,pa,pb, ... c psi,sig,theta,zeta]' c c INPUTS: param: initial parameter vector, param c c OUTPUTS: mm: vector of moments c c OTHER (LINPACK) blas.f dgefa.f, dgesl.f, dgeco.f c SUBROUTINES: (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, cumc3db.f c c Ellen McGrattan, 9-29-93 c Revised, 9-29-93, ERM c subroutine cumcmom(param,mom) parameter (na=54,nx=5,ny=2,nz=2,ns=4,nas=432,nxm=3,nym=3, c nzm=3,ne1=10,nw=20000) 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 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), param(12), foc(4), dfoc(4,4), del(4,4), c mom integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), kk(ne1), c ind(ne1*ns) c alpha = param(1) beta = param(2) delta = param(3) gamma = param(4) lam = param(5) omega = param(6) pa = param(7) pb = param(8) psi = param(9) sig = param(10) theta = param(11) zeta = param(12) ginv(1) = 0.3e-04 ginv(2) = 0.433e-01 ginv(3) = 0.59e-03 ginv(4) = 0.367e-02 cg(1) = 0.08 cg(2) = 0.19 cg(3) = 0.10 cg(4) = 0.13 army(1) = 0.002 army(2) = 0.0674 army(3) = 0.011 army(4) = 0.020 mcp(1,1) = 0.962 mcp(1,2) = 0.000 mcp(1,3) = 0.019 mcp(1,4) = 0.019 mcp(2,1) = 0.000 mcp(2,2) = 0.660 mcp(2,3) = 0.000 mcp(2,4) = 0.340 mcp(3,1) = 0.043 mcp(3,2) = 0.000 mcp(3,3) = 0.895 mcp(3,4) = 0.062 mcp(4,1) = 0.000 mcp(4,2) = 0.115 mcp(4,3) = 0.400 mcp(4,4) = 0.485 xa(1) = 0.01 xa(2) = 0.1 xa(3) = 0.2 xa(4) = 0.6 xa(5) = 1.0 xa(6) = 2.0 ya(1) = 0.6 ya(2) = 0.7 ya(3) = 0.8 za(1) = 0.900000E+00 za(2) = 0.100000E+01 za(3) = 0.110000E+01 do 10 i=1,nx nxp(i) = 2 10 continue do 20 i=1,ny nyp(i) = 2 20 continue do 30 i=1,nz nzp(i) = 2 30 continue delt1 = 1.d0-delta gamm1 = 1.d0-gamma ome1 = 1.d0-omega thet1 = 1.d0-theta cwgt = gamma*ome1 cwgt1 = cwgt-1.d0 hwgt = gamm1*ome1 hwgt1 = hwgt-1.d0 tiny = 0.000000001d0 pi = 3.14159265d0 crit = 0.000001d0 maxit = 20 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 c do 40 l=1,ns ssg = ginv(l)/delta ssk = 0.8 ssc = 0.2 sse = 0.4 ssh = 0.6 do 45 i=1,4 do 46 j=1,4 del(i,j) =0.d0 46 continue del(i,i) = 1.0e-4 45 continue do 50 t=1,10 d = (sse-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb zkk = lam*ssk**theta c0 = ssc*(1.d0-ssh)**(hwgt/cwgt1) foc(1) = gamm1/gamma*sse**theta*ssc-zkk*(1.d0-ssh) foc(2) = ssc**cwgt*(1.d0-ssh)**hwgt/ome1-c0**cwgt/ome1- c dpn-gamma*ssc**cwgt1*(1.d0-ssh)**hwgt*(ssc-c0- c zkk*thet1*sse**(-theta)*ssh) foc(3) = sse*ssc+(1.d0-sse)*c0+delta*ssk+cg(l)- c zkk*ssh*sse**thet1 foc(4) = beta*(delt1+theta*zkk/ssk*sse**thet1*ssh)-1.d0 do 55 i=1,4 dssk = ssk + del(1,i) dssc = ssc + del(2,i) dsse = sse + del(3,i) dssh = ssh + del(4,i) zkk = lam*dssk**theta c0 = dssc*(1.d0-dssh)**(hwgt/cwgt1) d = (dsse-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb foc1p = gamm1/gamma*dsse**theta*dssc-zkk*(1.d0-dssh) foc2p = dssc**cwgt*(1.d0-dssh)**hwgt/ome1-c0**cwgt/ome1- c dpn-gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh) foc3p = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1 foc4p = beta*(delt1+theta*zkk/dssk*dsse**thet1*dssh)-1.d0 dssk = ssk - del(1,i) dssc = ssc - del(2,i) dsse = sse - del(3,i) dssh = ssh - del(4,i) zkk = lam*dssk**theta c0 = dssc*(1.d0-dssh)**(hwgt/cwgt1) d = (dsse-pa)*pb dpn = 3.d0*d*(d+dsign(d,0.d0))*pb foc1m = gamm1/gamma*dsse**theta*dssc-zkk*(1.d0-dssh) foc2m = dssc**cwgt*(1.d0-dssh)**hwgt/ome1-c0**cwgt/ome1- c dpn-gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh) foc3m = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1 foc4m = beta*(delt1+theta*zkk/dssk*dsse**thet1*dssh)-1.d0 dfoc(1,i) = 0.5*(foc1p-foc1m)/del(i,i) dfoc(2,i) = 0.5*(foc2p-foc2m)/del(i,i) dfoc(3,i) = 0.5*(foc3p-foc3m)/del(i,i) dfoc(4,i) = 0.5*(foc4p-foc4m)/del(i,i) 55 continue call dgeco(dfoc,4,4,iwork,rcond,work) if (rcond.eq.0.0) c call xerror( 'dgeco error -- singular matrix ',55,-4,0) call dgesl(dfoc,4,4,iwork,foc,0) ssk = ssk - foc(1) ssc = ssc - foc(2) sse = sse - foc(3) ssh = ssh - foc(4) 50 continue ssi = delta*ssk ssz = 0.d0 ssf = lam*ssk**theta*sse**thet1*ssh write(*,*) 'steady states of k+kg,c1,e,h' write(*,805) ssk,ssc,sse,ssh n = 1 do 65 k=1,nz1 do 65 j=1,ny1 do 65 i=1,nx1 a0(n,l) = 0.8*ssc/ssf*lam*za(k)*xa(i)**theta* c sse**thet1*ssh b0(n,l) = sse+.1-.1/ssk*xa(i) n = n+1 65 continue 40 continue c 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 do 130 j=1,nas df(i,j) = 0.d0 130 continue 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 tem3 = u(i)*v(j)*w(k) do 340 mm=1,nas tem1 = tem3*work(mm) tem2 = tem3*wrk2(mm) df(i1,mm) = df(i1,mm) + tem1*xp*yp*zp df(i2,mm) = df(i2,mm) + tem2*xp*yp*zp df(i3,mm) = df(i3,mm) + tem1*xm*yp*zp df(i4,mm) = df(i4,mm) + tem2*xm*yp*zp df(i5,mm) = df(i5,mm) + tem1*xm*ym*zp df(i6,mm) = df(i6,mm) + tem2*xm*ym*zp df(i7,mm) = df(i7,mm) + tem1*xp*ym*zp df(i8,mm) = df(i8,mm) + tem2*xp*ym*zp df(i9,mm) = df(i9,mm) + tem1*xp*yp*zm df(i10,mm) = df(i10,mm) + tem2*xp*yp*zm df(i11,mm) = df(i11,mm) + tem1*xm*yp*zm df(i12,mm) = df(i12,mm) + tem2*xm*yp*zm df(i13,mm) = df(i13,mm) + tem1*xm*ym*zm df(i14,mm) = df(i14,mm) + tem2*xm*ym*zm df(i15,mm) = df(i15,mm) + tem1*xp*ym*zm df(i16,mm) = df(i16,mm) + tem2*xp*ym*zm 340 continue 240 continue 150 continue 140 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 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 805 format (1X,5(1X, F12.5)) mom = ssk return end