CFPP$ UNROLL (32) F c c PROGRAM: CUMC3DBi.F c c TITLE: implementation of the finite-element method for model c with capacity utilization and employment costs, in c ``The Effects of Government purchases on Employment c and 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)-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)-P'(n^l))*(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 -alpha(n^l(k,e,z)-e)^2-P(n^l(k,e,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 (NOTE: LQ solution used as initial guess for c_a's,n_a's) 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)}*(1-a[t])|x[0]] c c subject to c c c[t]+i[t]+m[t]+p[t] = z[t]*F(k[t]+kg[t],n[t],h[t]) c -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 p[t] = P(n[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],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 TEST CASE: set alpha=0 and compare results to Kydland-Prescott c c INPUTS: (a) cumc3dai.inp: alpha,beta,delta,gamma,lam,omega, ... c pa,pb,psi,sig,theta,zeta c ig(j),cg(j),a(j), j=1:ns c mcp(i,j) i=1:ns, j=1:ns c ssk,ssc,ssn,ssh c where alpha, etc. are parameters of utility/technology c ig's,cg's,a's are values for Markov chain c mcp's are the transition probabilities c ssk, etc. are initial guesses of the steady c for k, c, n, h (if 0, then program chooses) c c OUTPUTS: (a) cumc3d.dat: input file for cumc3d.m (Matlab file) 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 or cumc3db.f c c Ellen McGrattan, 10-24-93 c Revised, 10-24-93, ERM c parameter (na=126,nx=13,ny=2,nz=2,ns=4,nas=1008,nxm=3,nym=3,nzm=3, c ne1=10,nw=20000,nsv=4,ndv=4,nv=8,nt=1000) 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) integer iwork(nw), nxp(nx), nyp(ny), nzp(nz), kk(ne1), c ind(ne1*ns), istat(nt) c open(unit=5, file='cumc3dbi.inp') open(unit=10, file='cumc3db.dat') c read(5,*) alpha read(5,*) beta read(5,*) delta read(5,*) gamma read(5,*) lam read(5,*) omega read(5,*) pa read(5,*) pb read(5,*) psi read(5,*) sig read(5,*) theta read(5,*) zeta do 1 i=1,ns read(5,*) ginv(i),cg(i),army(i) 1 continue do 5 i=1,ns read(5,*) (mcp(i,j),j=1,ns) 5 continue read(5,*) ssk,ssc,sse,ssh read(5,*) lqonly c c 3-dimensional domain (xa,ya,za) defined by: c xa(1) = 0.3 xa(2) = 0.325 xa(3) = 0.35 xa(4) = 0.375 xa(5) = 0.4 xa(6) = 0.425 xa(7) = 0.45 xa(8) = 0.475 xa(9) = 0.5 xa(10) = 0.525 xa(11) = 0.55 xa(12) = 0.575 xa(13) = 0.6 xa(14) = 0.7 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 nx1 = nx+1 ny1 = ny+1 nz1 = nz+1 maxit = 20 maxit0 = 20 if (lqonly.eq.1) maxit0 = 0 c do 40 i=1,8 do 40 j=1,8 del(i,j) =0.d0 40 continue if (ssk.lt.tiny) ssk = 0.4 if (ssc.lt.tiny) ssc = 0.1 if (sse.lt.tiny) sse = 0.4 if (ssh.lt.tiny) ssh = 0.42 do 50 l=1,ns 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 pn = d*d*d+dabs(d*d*d) 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 gamma*ssc**cwgt1*(1.d0-ssh)**hwgt*(ssc-c0- c zkk*thet1*sse**(-theta)*ssh+dpn) foc(3) = sse*ssc+(1.d0-sse)*c0+delta*ssk+cg(l)- c zkk*ssh*sse**thet1+pn 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 pn = d*d*d+dabs(d*d*d) 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 gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh+dpn) foc3p = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1+pn 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 pn = d*d*d+dabs(d*d*d) 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 gamma*dssc**cwgt1*(1.d0-dssh)**hwgt*(dssc-c0- c zkk*thet1*dsse**(-theta)*dssh+dpn) foc3m = dsse*dssc+(1.d0-dsse)*c0+delta*dssk+cg(l)- c zkk*dssh*dsse**thet1+pn 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 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) 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 pn)/(1.d0-dsse) up = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) 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 pn)/(1.d0-dsse) um = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) 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)-pn)/(1.d0-dsse) upp = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) 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)-pn)/(1.d0-dsse) upm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) 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)-pn)/(1.d0-dsse) ump = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) 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)-pn)/(1.d0-dsse) umm = (1.d0-army(l))*((dsse*dssc**cwgt*(1.d0-dssh)**hwgt+ c (1.d0-dsse)*c0**cwgt)/ome1) 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,maxit0 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))-pn 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 pn = d*d*d+dabs(d*d*d) 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))-dpn 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) 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))+dpn 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+tem1*hwgt*mu/h1* c dhna-mu*(d2pn-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+mu*(dpn+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