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, dgeco.f, dgedi.f dgefa.f, dgesl.f c SUBROUTINES: (McGrattan) qgausl.f c c REFERENCES: Eischen, Jeff, (notes on the Finite Element Method) 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, 10-2-93, ERM c parameter (na=63,nx=6,ny=2,nz=2,ns=2,nas=252,nxm=3,nym=3,nzm=3, c ne1=5,nw=20000,nsv=4,ndv=4,nv=8,nt=1000,nm=10) 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), foc(4), dfoc(4,4), del(8,8), c mom, amat(nsv,nsv), bmat(nsv,ndv), qmat(nsv,nsv), c rmat(ndv,ndv), wmat(nsv,ndv), aa(nsv,nsv), bb(nsv,nsv), c cc(nsv,nsv), dd(nsv,nsv), ee(nsv,nsv), qwr(nv,nv), c f0(ndv,nsv), f1(ndv,nsv), du(8), du2(8,8), ri(ndv,ndv), c rw(ndv,nsv), rb(ndv,ndv), det(2), lamhld integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), kk(ne1), c ind(ne1*ns), istat(nt) c open(unit=10, file='cumc3d.dat') beta = 0.99 delta = 0.0205 lam = 1.0 omega = 1.001 c alpha = 1e-02 gamma = 0.2204206508139047 gamma = 0.2 pa = 0.5400883836248350 pa = 0.3 pb = 1.031133977905275 psi = 0.983 psi = 0.95 sig = 0.0127 theta = 0.255 zeta = 0.1997631855911816 zeta = -0.01 ginv(1) = 0.0002 ginv(2) = 0.0004 cg(1) = 0.025 cg(2) = 0.04 c cg(1) = 0.03 cg(2) = 0.045 pa =0.48 gamma = 0.25 theta = 0.26 c write(*,*) ' alpha = ',alpha write(*,*) ' beta = ',beta write(*,*) ' delta = ',delta write(*,*) ' gamma = ',gamma write(*,*) ' lam = ',lam write(*,*) ' omega = ',omega write(*,*) ' pa = ',pa write(*,*) ' pb = ',pb write(*,*) ' psi = ',psi write(*,*) ' sig = ',sig write(*,*) ' theta = ',theta write(*,*) ' zeta = ',zeta write(*,*) ' gi1 = ',ginv(1) write(*,*) ' gi2 = ',ginv(2) write(*,*) ' gc1 = ',cg(1) write(*,*) ' gc2 = ',cg(2) c army(1) = 0.012975895 army(2) = 0.024033008 mcp(1,1) = 0.991397304 mcp(2,1) = 0.016234681 mcp(1,2) = 0.008602696 mcp(2,2) = 0.983765319 xa(1) = 2.0 xa(2) = 2.5 xa(3) = 3.0 xa(4) = 3.5 xa(5) = 4.0 xa(6) = 4.5 xa(7) = 5.0 ya(1) = 0.35 ya(2) = 0.5 ya(3) = 0.65 za(1) = 0.90000E+00 za(2) = 0.100000E+01 za(3) = 0.11000E+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.0000001d0 maxit = 20 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 c do 40 i=1,8 do 40 j=1,8 del(i,j) =0.d0 40 continue do 50 l=1,ns ssk = 3.0 ssc = 0.25 sse = 0.40 ssh = 0.42 del(1,1) = ssk*1.0e-5 del(2,2) = ssc*1.0e-5 del(3,3) = sse*1.0e-5 del(4,4) = ssh*1.0e-5 do 60 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 70 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) 70 continue call dgeco(dfoc,4,4,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when computing ', c 'steady state values' call dgesl(dfoc,4,4,iwork,foc,0) ssk = ssk - foc(1) ssc = ssc - foc(2) sse = sse - foc(3) ssh = ssh - foc(4) 60 continue write(*,1001) ' The steady state values: ',ssk,ssc,sse,ssh ssi = delta*ssk ssn = sse ssz = 0.d0 ssf = lam*ssk**theta*sse**thet1*ssh del(1,1) = ssk*1.0e-4 del(2,2) = ssn*1.0e-4 del(3,3) = 1.0e-5 del(4,4) = 1.0e-4 del(5,5) = ssc*1.0e-4 del(6,6) = ssi*1.0e-4 del(7,7) = sse*1.0e-4 del(8,8) = ssh*1.0e-4 c d = (sse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = ssc*(1.d0-ssh)**(hwgt/cwgt1) uu = (1.d0-army(l))*((sse*ssc**cwgt*(1.d0-ssh)**hwgt+ c (1.d0-sse)*c0**cwgt)/ome1-pn) du(4) = 0.d0 do 80 i=1,8 du(i) = 0.d0 do 80 j=1,8 du2(i,j) = 0.d0 80 continue do 90 i=1,8 dssk = ssk + del(1,i) dssn = ssn + del(2,i) dssz = ssz + del(3,i) dssc = ssc + del(5,i) dssi = ssi + del(6,i) dsse = sse + del(7,i) dssh = ssh + del(8,i) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)*(dsse-dssn))/ c (1.d0-dsse) up = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk - del(1,i) dssn = ssn - del(2,i) dssz = ssz - del(3,i) dssc = ssc - del(5,i) dssi = ssi - del(6,i) dsse = sse - del(7,i) dssh = ssh - del(8,i) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)*(dsse-dssn))/ c (1.d0-dsse) um = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) du(i)= 0.5*(up-um)/del(i,i) do 100 j=1,i dssk = ssk + del(1,i) + del(1,j) dssn = ssn + del(2,i) + del(2,j) dssz = ssz + del(3,i) + del(3,j) dssc = ssc + del(5,i) + del(5,j) dssi = ssi + del(6,i) + del(6,j) dsse = sse + del(7,i) + del(7,j) dssh = ssh + del(8,i) + del(8,j) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) upp = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk + del(1,i) - del(1,j) dssn = ssn + del(2,i) - del(2,j) dssz = ssz + del(3,i) - del(3,j) dssc = ssc + del(5,i) - del(5,j) dssi = ssi + del(6,i) - del(6,j) dsse = sse + del(7,i) - del(7,j) dssh = ssh + del(8,i) - del(8,j) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) upm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk - del(1,i) + del(1,j) dssn = ssn - del(2,i) + del(2,j) dssz = ssz - del(3,i) + del(3,j) dssc = ssc - del(5,i) + del(5,j) dssi = ssi - del(6,i) + del(6,j) dsse = sse - del(7,i) + del(7,j) dssh = ssh - del(8,i) + del(8,j) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) ump = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) dssk = ssk - del(1,i) - del(1,j) dssn = ssn - del(2,i) - del(2,j) dssz = ssz - del(3,i) - del(3,j) dssc = ssc - del(5,i) - del(5,j) dssi = ssi - del(6,i) - del(6,j) dsse = sse - del(7,i) - del(7,j) dssh = ssh - del(8,i) - del(8,j) d = (dsse-pa)*pb pn = d*d*d+dabs(d*d*d) c0 = (lam*dexp(dssz)*dssk**theta*dsse**thet1*dssh- c dssi-cg(l)-dsse*dssc-alpha*(dsse-dssn)* c (dsse-dssn))/(1.d0-dsse) umm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1-pn) du2(i,j) = 0.25*(upp-upm-ump+umm)/(del(i,i)*del(j,j)) du2(j,i) = du2(i,j) 100 continue 90 continue do 110 i=1,8 do 110 j=1,8 qwr(i,j) = 0.d0 110 continue do 120 j=1,8 qwr(1,j) = 0.5*du2(1,j) qwr(2,j) = 0.5*du2(2,j) qwr(3,j) = 0.5*du2(3,j) qwr(4,j) = 0.5*du2(4,j) qwr(5,j) = 0.5*du2(5,j) qwr(6,j) = 0.5*du2(6,j) qwr(7,j) = 0.5*du2(7,j) qwr(8,j) = 0.5*du2(8,j) 120 continue do 130 i=1,8 work(i) = 0.5*(du(i)-du2(i,1)*ssk-du2(i,2)*ssn-du2(i,5)* c ssc-du2(i,6)*ssi-du2(i,7)*sse-du2(i,8)*ssh) qwr(4,i) = qwr(4,i) + work(i) qwr(i,4) = qwr(i,4) + work(i) 130 continue qwr(4,4) = qwr(4,4) + uu-(work(1)+0.5*du(1))*ssk- c (work(2)+0.5*du(2))*ssn-(work(5)+0.5*du(5))*ssc- c (work(6)+0.5*du(6))*ssi-(work(7)+0.5*du(7))*sse- c (work(8)+0.5*du(8))*ssh do 140 i=1,4 do 140 j=1,4 qmat(i,j) = qwr(i,j) amat(i,j) = 0.d0 140 continue do 150 i=1,4 do 150 j=1,4 wmat(i,j) = qwr(i,j+4) bmat(i,j) = 0.d0 150 continue do 160 i=1,4 do 160 j=1,4 rmat(i,j) = qwr(i+4,j+4) 160 continue sb = dsqrt(beta) amat(1,1) = delt1*sb amat(3,3) = psi*sb amat(4,4) = sb bmat(1,2) = sb bmat(2,3) = sb c c Convert to the no-discount, no-cross product case: c do 170 j=1,4 ri(1,j) = rmat(1,j) ri(2,j) = rmat(2,j) ri(3,j) = rmat(3,j) ri(4,j) = rmat(4,j) 170 continue call dgeco(ri,ndv,ndv,iwork,rcond,work) call dgedi(ri,ndv,ndv,iwork,det,work,11) do 180 i=1,4 do 180 j=1,4 rw(i,j) = ri(i,1)*wmat(j,1)+ri(i,2)*wmat(j,2)+ c ri(i,3)*wmat(j,3)+ri(i,4)*wmat(j,4) 180 continue do 190 i=1,4 do 190 j=1,4 qmat(i,j) = qmat(i,j) -wmat(i,1)*rw(1,j)-wmat(i,2)* c rw(2,j)-wmat(i,3)*rw(3,j)-wmat(i,4)*rw(4,j) 190 continue do 200 i=1,4 do 200 j=1,4 amat(i,j) = amat(i,j) -bmat(i,1)*rw(1,j)-bmat(i,2)* c rw(2,j)-bmat(i,3)*rw(3,j)-bmat(i,4)*rw(4,j) 200 continue c c Doubling algorithm used to solve fixed point of Riccati equation c do 210 i=1,nsv do 210 j=1,ndv sum1 = 0.d0 do 220 k=1,ndv sum1 = sum1 + bmat(i,k)*ri(k,j) 220 continue wmat(i,j) = sum1 210 continue do 230 i=1,nsv do 230 j=1,nsv sum1 = 0.d0 do 240 k=1,ndv sum1 = sum1 + wmat(i,k)*bmat(j,k) 240 continue bb(i,j) = sum1 aa(i,j) = amat(i,j) 230 continue do 250 i=1,ndv do 250 j=1,nsv sum1 = 0.d0 do 260 k=1,nsv sum1 = sum1 + bmat(k,i)*qmat(k,j) 260 continue wmat(j,i) = sum1 250 continue do 270 i=1,ndv do 270 j=1,nsv sum1 = 0.d0 do 280 k=1,nsv sum1 = sum1 + wmat(k,i)*amat(k,j) 280 continue f0(i,j) = sum1 270 continue do 290 i=1,ndv do 290 j=1,ndv sum1 = 0.d0 do 300 k=1,nsv sum1 = sum1 + wmat(k,i)*bmat(k,j) 300 continue rb(i,j) = sum1+rmat(i,j) 290 continue call dgeco(rb,ndv,ndv,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when computing ', c 'linear-quadratic decision function' do 310 j=1,nsv call dgesl(rb,ndv,ndv,iwork,f0(1,j),0) 310 continue d1 = 1.d0 do 320 it=1,maxit do 330 i=1,nsv do 340 j=1,nsv sum1 = 0.d0 do 350 k=1,nsv sum1 = sum1 + bb(i,k)*qmat(k,j) 350 continue ee(i,j) = sum1 340 continue ee(i,i) = ee(i,i) + 1.d0 330 continue call dgeco(ee,nsv,nsv,iwork,rcond,work) call dgedi(ee,nsv,nsv,iwork,det,work,01) do 360 i=1,nsv do 360 j=1,nsv sum1 = 0.d0 do 370 k=1,nsv sum1 = sum1 + bb(i,k)*aa(j,k) 370 continue cc(i,j) = sum1 360 continue do 380 i=1,nsv do 380 j=1,nsv sum1 = 0.d0 do 390 k=1,nsv sum1 = sum1 + aa(i,k)*ee(k,j) 390 continue dd(i,j) = sum1 380 continue do 400 i=1,nsv do 400 j=1,nsv sum1 = 0.d0 do 410 k=1,nsv sum1 = sum1 + dd(i,k)*cc(k,j) 410 continue bb(i,j) = bb(i,j)+sum1 400 continue do 420 i=1,nsv do 420 j=1,nsv sum1 = 0.d0 do 430 k=1,nsv sum1 = sum1 + ee(i,k)*aa(k,j) 430 continue cc(i,j) = sum1 420 continue do 440 i=1,nsv do 440 j=1,nsv sum1 = 0.d0 do 450 k=1,nsv sum1 = sum1 + aa(k,i)*qmat(k,j) 450 continue dd(i,j) = sum1 440 continue do 460 i=1,nsv do 460 j=1,nsv sum1 = 0.d0 do 470 k=1,nsv sum1 = sum1 + dd(i,k)*cc(k,j) 470 continue qmat(i,j) = qmat(i,j)+sum1 460 continue do 480 i=1,nsv do 480 j=1,nsv sum1 = 0.d0 do 490 k=1,nsv sum1 = sum1 + ee(i,k)*aa(k,j) 490 continue cc(i,j) = sum1 480 continue do 500 i=1,nsv do 500 j=1,nsv sum1 = 0.d0 do 510 k=1,nsv sum1 = sum1 + aa(i,k)*cc(k,j) 510 continue dd(i,j) = sum1 500 continue do 520 i=1,nsv do 520 j=1,nsv aa(i,j) = dd(i,j) 520 continue do 530 i=1,ndv do 530 j=1,nsv sum1 = 0.d0 do 540 k=1,nsv sum1 = sum1 + bmat(k,i)*qmat(k,j) 540 continue wmat(j,i) = sum1 530 continue do 550 i=1,ndv do 550 j=1,nsv sum1 = 0.d0 do 560 k=1,nsv sum1 = sum1 + wmat(k,i)*amat(k,j) 560 continue f1(i,j) = sum1 550 continue do 570 i=1,ndv do 570 j=1,ndv sum1 = 0.d0 do 580 k=1,nsv sum1 = sum1 + wmat(k,i)*bmat(k,j) 580 continue rb(i,j) = sum1+rmat(i,j) 570 continue call dgeco(rb,ndv,ndv,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when computing ', c 'linear-quadratic decision function' do 590 j=1,nsv call dgesl(rb,ndv,ndv,iwork,f1(1,j),0) 590 continue sum1 = 0.d0 do 600 i=1,ndv do 600 j=1,nsv sum1 = sum1+(f1(i,j)-f0(i,j))*(f1(i,j)-f0(i,j)) 600 continue if (dsqrt(sum1)/float(ndv*nsv).lt.crit) go to 997 do 610 i=1,ndv do 610 j=1,nsv f0(i,j) = f1(i,j) 610 continue 320 continue 997 do 620 i=1,ndv do 620 j=1,nsv f1(i,j) = -f1(i,j)-rw(i,j) 620 continue n = 1 do 630 k=1,nz1 do 630 j=1,ny1 do 630 i=1,nx1 a0(n,l) = f1(1,1)*xa(i)+f1(1,2)*ya(j)+ c f1(1,3)*dlog(za(k))+f1(1,4) b0(n,l) = f1(3,1)*xa(i)+f1(3,2)*ya(j)+ c f1(3,3)*dlog(za(k))+f1(3,4) n = n+1 630 continue 50 continue c sum1 = 0.d0 sum2 = 0.d0 call qgausl(ne1,-2.88d0*sig,2.88d0*sig,e1,w1) do 640 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) 640 continue do 650 i=1,ne1 wgt(i) = wgt(i)/sum1 sum2 = sum2 + wgt(i) 650 continue do 660 i=1,nx call qgausl(nxp(i),xa(i),xa(i+1),x,u) do 660 j=1,nxp(i) xx(j,i) = x(j) wx(j,i) = u(j) 660 continue do 670 i=1,ny call qgausl(nyp(i),ya(i),ya(i+1),y,v) do 670 j=1,nyp(i) yy(j,i) = y(j) wy(j,i) = v(j) 670 continue do 680 i=1,nz call qgausl(nzp(i),za(i),za(i+1),z,w) do 680 j=1,nzp(i) zz(j,i) = z(j) wz(j,i) = w(j) 680 continue ne = nx*ny*nz nxy = nx*ny sumres = 1.d0 do 690 it=1,maxit do 700 i=1,nas f(i) = 0.d0 do 700 j=1,nas df(i,j) = 0.d0 700 continue do 710 l=1,ns army1 = 1.d0-army(l) do 720 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 730 i=1,n1 x(i) = xx(i,ix) u(i) = wx(i,ix) 730 continue do 740 j=1,n2 y(j) = yy(j,iy) v(j) = wy(j,iy) 740 continue do 750 k=1,n3 z(k) = zz(k,iz) w(k) = wz(k,iz) zp = (z2 - z(k))/dz zm = (z(k) - z1)/dz do 750 j=1,n2 yp = (y2 - y(j))/dy ym = (y(j) - y1)/dy do 750 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 CFPP$ UNROLL do 760 m=1,nx j1 = idint( dmax1(xt-xa(m),0.d0)/(xt-xa(m))+.2 ) ii = ii + j1 760 continue ii = max(1,ii-1) jj = 1 CFPP$ UNROLL do 770 m=1,ny j1 = idint( dmax1(yt-ya(m),0.d0)/(yt-ya(m))+.2) jj = jj + j1 770 continue jj = max(1,jj-1) CFPP$ UNROLL do 780 m=1,ne1 zzt(m) = dexp(psi*dlog(z(k))+e1(m)) kk(m) = 1 CFPP$ UNROLL do 790 mm=1,nz j1 = idint( dmax1(zzt(m)-za(mm),0.d0)/ c (zzt(m)-za(mm))+.2) kk(m)= kk(m) + j1 790 continue kk(m) = max(1,kk(m)-1) 780 continue do 800 m=1,34 sum(m) = 0.d0 800 continue do 810 m=1,nas work(m) = 0.d0 wrk2(m) = 0.d0 810 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)) CFPP$ UNROLL do 820 m=1,ns CFPP$ UNROLL do 820 mm=1,ne1 nct = (m-1)*ne1+mm zt = zzt(mm) 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 820 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 CFPP$ UNROLL do 830 m=1,ns CFPP$ UNROLL do 830 mm=1,ne1 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) 830 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 840 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 840 continue 750 continue 720 continue 710 continue call dgeco(df,nas,nas,iwork,rcond,work) if (rcond.eq.0.0) c write(*,*) 'WARNING: singular matrix when solving df*x=f' call dgesl(df,nas,nas,iwork,f,0) sum3 = 0.d0 do 850 j=1,ns do 850 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 850 continue sumres = dsqrt(sum3)/float(nas) write(*,*) ' at iteration ', it, ' the residual is ', sumres do 860 j=1,ns do 860 i=1,na a0(i,j) = a1(i,j) b0(i,j) = b1(i,j) 860 continue if (sumres.lt.crit) go to 999 690 continue 999 continue c c print results c 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) psi write(10,801) sig write(10,801) theta write(10,801) zeta write(10,812) ns do 870 i=1,ns write(10,801) ginv(i) write(10,801) cg(i) write(10,801) army(i) 870 continue do 880 j=1,ns do 880 i=1,ns write(10,801) mcp(i,j) 880 continue write(10,812) nx1 write(10,812) ny1 write(10,812) nz1 do 890 i=1,nx1 write(10,801) xa(i) 890 continue do 900 i=1,ny1 write(10,801) ya(i) 900 continue do 910 i=1,nz1 write(10,801) za(i) 910 continue do 920 j=1,ns do 920 i=1,na write(10,801) a0(i,j) 920 continue do 930 j=1,ns do 930 i=1,na write(10,801) b0(i,j) 930 continue c 801 format (4X, E12.6) 802 format (2(4X, F12.5)) 803 format (3(4X, F12.5)) 804 format (4(4X, F12.5)) 805 format (5(4X, F12.5)) 806 format (6(2X, F11.5)) 807 format (7(2X, F10.5)) 808 format (8(2X, F9.5)) 811 format (4X, E12.6, 4X, I4) 812 format (4X, I4) 813 format (4X, 3(1X, E9.3)) 1001 format (A27,4(F8.3,1X)) 1002 format (A27,E10.5,1X,2(F8.3,1X)) stop end