subroutine NEWT(x,n,check,done) implicit double precision (a-h,o-z) c c Given an initial guess X(N) for a root in N dimensions, find the root c by a globally convergent Newton's method. The vector of functions to be c zeroed, called FVEC(N) below, is returned by a user-supplied subroutine c that must be called FUNCV and have the declaration c subroutine funcv(n,x,fvec) c The output quantity CHECK is false on a normal return and true if the c routine has converged to a local minimum of the function FMIN defined c below. In this case try restarting from a different initial guess. c c Parameters: NP is the maximum expected value of n; MAXITS is the maximum c number of iterations; TOLF sets the convergence criterion on function c values; TOLMIN sets the criterion for deciding whether spurious convergence c to a minimum of fmin has occurred; TOLX is the convergence criterion on dx; c STMX is the scaled maximum step length allowed in line searchec c INTEGER n,nn,NP,MAXITS LOGICAL check DOUBLE PRECISION x(n),fvec,TOLF,TOLMIN,TOLX,STPMX PARAMETER (NP=40,MAXITS=200,TOLF=1.d-4,TOLMIN=1.d-2*TOLF, 1 TOLX=1.d-3*TOLF,STPMX=100.d0) COMMON/newtv/ fvec(NP),nn SAVE /newtv/ INTEGER i,its,j,indx(NP) DOUBLE PRECISION d,den,f,fold,stpmax,sum,temp,test,fjac(NP,NP), 1 g(NP),p(NP),xold(NP),fmin EXTERNAL fmin nn=n f=fmin(x) test=0.d0 do 11 i=1,n if(abs(fvec(i)).gt.test) test=abs(fvec(i)) 11 continue if (test.lt.0.01d0*TOLF) return sum=0.d0 do 12 i=1,n sum=sum+x(i)**2 12 continue stpmax=STPMX*max(sqrt(sum),dble(n)) do 21 its=1,MAXITS f=fmin(x) call fdjac(n,x,fvec,NP,fjac) do 14 i=1,n sum=0.d0 do 13 j=1,n sum=sum+fjac(j,i)*fvec(j) 13 continue g(i)=sum 14 continue do 15 i=1,n xold(i)=x(i) 15 continue fold=f do 16 i=1,n p(i)=-fvec(i) 16 continue call ludcmp(fjac,n,NP,indx,d) call lubksb(fjac,n,NP,indx,p) c call lnsrch(n,xold,fold,g,p,x,f,stpmax,check,fmin) test=0.d0 do 17 i=1,n xold(i)=x(i) x(i)=x(i)+p(i) if (abs(fvec(i)).gt.test) test=abs(fvec(i)) 17 continue write (6,'(i5,1p6e12.4)') its,(fvec(iout),iout=1,n) write (6,'(10x,1p6e12.4)') (x(iout),iout=1,n) write (6,'(10x,1p6e12.4)') (p(iout),iout=1,n) c if (test.lt.TOLF) then if (test.lt.done) then check=.false. return endif c if (check) then c test=0.d0 c den=max(f,0.5d0*n) c do 18 i=1,n c temp=abs(g(i))*max(abs(x(i)),1.d0)/den c if (temp.gt.test) test=temp c18 continue c if (test.lt.TOLMIN)then c check=.true. c else c check=.false. c endif c return c endif test=0.d0 21 continue pause 'MAXITS exceeded in NEWT' end c subroutine fdjac(n,x,fvec,np,df) implicit double precision (a-h,o-z) c c Computes forward difference approximation to Jacobian. On input, X(N) is c the point at which the Jacobian is to be evaluated, FVEC(N) is the vector c of function values at the point, and NP is the physical dimension of the c Jacobian array DF(N,N) which is output. subroutine funcv(n,x,f) is a c fixed-name, user supplied routine that returns the vector of functions c at X. c c Parameters: NMAX is the maximum value of n; EPS is the approxmate square c root of the machine precision c c requires FUNCV c INTEGER n,np,NMAX DOUBLE PRECISION df(np,np),fvec(n),x(n),EPS PARAMETER (NMAX=40,EPS=1.d-6) c PARAMETER (NMAX=40,EPS=1.d-8) INTEGER i,j DOUBLE PRECISION h,temp,f(NMAX) do 12 j=1,n temp=x(j) h=EPS*abs(temp) if (h.eq.0.d0) h=EPS x(j)=temp+h h=x(j)-temp call funcv(n,x,f) x(j)=temp do 11 i=1,n df(i,j)=(f(i)-fvec(i))/h 11 continue 12 continue return end c SUBROUTINE lnsrch(n,xold,fold,g,p,x,f,stpmax,check,func) implicit double precision (a-h,o-z) c c Given an n-dimensional point xold(n), the value of the function and gradient c there, fold and g(n), and a direction p(n), finds a new point x(n) along the c direction p from xold where the function func has decreased sufficiently. The c new function value is returned in f. stpmax in an input quantity that limits c the length of the steps so that you do not try to evaluate the function in c regions where it is undefined or subject to overflow. p is usually the Newton c direction. The output quantity check is false on a normal exit. It is true c when x is too close to xold. In a minimization algorithm, this usually c signals convergence and can be ignored. However, in a zero-finding algorithm c the calling program should check whether the convergence is spurious. c c Parameters: ALF ensures sufficient decrease in function value; TOLX is the c convergence criterion on Dx c INTEGER n LOGICAL check DOUBLE PRECISION f,fold,stpmax,g(n),p(n),x(n),xold(n), 1 func,ALF,TOLX PARAMETER(ALF=1.d-4,TOLX=1.d-7) EXTERNAL func INTEGER i DOUBLE PRECISION a,alam,alam2,alamin,b,disc,f2,fold2,rhs1, 1 rhs2,slope,sum,temp,test,tmplam c check=.false. sum=0.0 do 11 i=1,n sum=sum+p(i)*p(i) 11 continue sum=sqrt(sum) if (sum.gt.stpmax) then do 12 i=1,n p(i)=p(i)*stpmax/sum 12 continue endif slope=0.d0 do 13 i=1,n slope=slope+g(i)*p(i) 13 continue test=0.d0 do 14 i=1,n temp=abs(p(i))/max(abs(xold(i)),1.d0) if (temp.gt.test)test=temp 14 continue alamin=TOLX/test alam=1.d0 1 continue do 15 i=1,n x(i)=xold(i)+alam*p(i) 15 continue f=func(x) if (alam.lt.alamin) then do 16 i=1,n x(i)=xold(i) 16 continue check=.true. return elseif (f.le.fold+ALF*alam*slope) then return else if(alam.eq.1.d0) then tmplam=-slope/(2.d0*(f-fold-slope)) else rhs1=f-fold-alam*slope rhs2=f2-fold2-alam2*slope a=(rhs1/alam**2-rhs2/alam2**2)/(alam-alam2) b=(-alam2*rhs1/alam**2+alam*rhs2/alam2**2)/ 1 (alam-alam2) if (a.eq.0.d0) then tmplam=-slope/(2.d0*b) else disc=b*b-3.d0*a*slope tmplam=(-b+sqrt(disc))/(3.d0*a) endif if (tmplam.gt.0.5d0*alam)tmplam=0.5d0*alam endif endif alam2=alam f2=f fold2=fold alam=max(tmplam,0.1d0*alam) goto 1 end c FUNCTION fmin(x) implicit double precision (a-h,o-z) c c Returns f=1/2 F.F at X. subroutine funcv(n,x,f) is a fixed-name, user c supplied routine that returns the vector of functions at X. The common c block newtv communicates the function valuse back to newt c c requires: funcv c INTEGER n,NP DOUBLE PRECISION fmin,x(*),fvec PARAMETER (NP=40) COMMON/newtv/ fvec(NP),n SAVE /newtv/ INTEGER i DOUBLE PRECISION sum call funcv(n,x,fvec) sum=0.d0 do 11 i=1,n sum=sum+fvec(i)**2 11 continue fmin=0.5d0*sum return end