call test end !----------------------------------------------------------------------- !##################################BEGIN GN############################# ! variables and parameters of GN that might need to be accessed module GN_MOD implicit none ! stptol in: step size for relative convergence test. ! reltol, abstol in: value relative/absolute convergence test. ! Setting accuracies too small wastes evaluations. ! Setting accuracies too large may not get the best solution. ! derivstp in: the step for derivatives. ! Must be large enough to get some change in the function. ! Must be small enough for some accuracy in the derivative. ! iprint in: degree of printout, if value is: ! 0=none, 1=final f0 and x0, 2=iterations, 3=each x tried, ! 4=all except Jacobian and Hessian, 5=all. ! ihist in: if > 0, write history file HistFile using unit ihist ! limit in: maximum number of all evaluations allowed, approximate. ! tuning constants: ! ZLOW,ZHIGH change bounds of trust region (del) (.5,2). ! Small changes in the path set by them greatly affect results. ! ZCP Cauchy step size (1). ! ZCPMIN minimum Cauchy step size (0.1). ! ZCPMAX maximum Cauchy step size (1000). ! MXBAD number bad steps with rank-1 updates (1). ! NewtStep1st=.true. start with Newton step. Use on linear problems. ! mu0,del0,r0,H0,Jac0,iHOOKmax out: final values ! dfj/dxi=Jac0(j,i) ! iHOOKmax is max number of iterations in hook search ! iscale=0, no scaling ! iscale=1, variable scaling ! iscale=2, fixed scaling based on D0, which must be allocated ! and filled by user before call to GN ! FIXED PARAMETER real*8 :: ONE=1 ! OUTPUT INFORMATION integer, save :: iHOOKmax,iter,nconsecMax,nfcn real*8, save :: del0,hnorm,mu0 real*8, save, dimension(:), allocatable :: D0,r0 real*8, save, dimension(:,:), allocatable :: H0,Jac0 ! USER PARAMETERS character*30, save :: HistFile='GNHIST.TXT' logical :: NewtStep1st=.false. integer :: ihist=0,iprint=0,iscale=0,limit=10000 real*8 :: abstol=1d-11,derivstp=1d-4,reltol=3d-7,stptol=1d-4 ! GN TUNING CONSTANTS integer :: MXBAD=1 real*8 :: ZCONSEC=1.005d0,ZHIGH=2d0,ZLOW=0.65d0 real*8 :: ZCP=.6d0,ZCPMAX=1000,ZCPMIN=0.18d0 end module GN_MOD !----------------------------------------------------------------------- !+augmented Gauss-Newton nonlinear least squares subroutine GN(fcn,m,n,info,x0,f0) ! AUTHORS: Kenneth Klare (kklare@gmail.com) ! and Guthrie Miller (guthriemiller@gmail.com) ! OPEN SOURCE: may be freely used, with author achnowledgment, ! for any peaceful, beneficial purpose ! PURPOSE: Minimize sum of squares of residuals r using ! augmented Gauss-Newton step and Levenberg-Marquardt trust region. ! Uses finite-difference derivatives and 1-D Jacobian updates. ! KEYWORDS: Nonlinear Least squares, nonlinear minimization, robust, ! optimization, finite-difference derivatives, ! augmented Gauss-Newton, Levenberg-Marquardt. ! ARGUMENTS: ! fcn in: A user subroutine to compute r(1:m) from x(1:n). ! It is invoked: call FCN(m,n,x,r), where ! x in. n parameters to evaluate residuals. ! r out. m residuals to be reduced. DO NOT SQUARE THEM. ! It must not change m, n, nor x. ! Scale the output residuals: ! r(i)=(fit(x,i)-y(i))/yerror(i), where y is the data. ! To reject the given x, set the r very large, like sqrt(largest). ! The starting point (x0) must valid and not be rejected. ! Declare FCN EXTERNAL so the compiler knows it. ! An iteration has function and derivative evaluations. ! Updates use function evaluations only. ! The program is designed to minimize evaluation count. ! m in: number of equations or fitting residuals, the size of r0. ! n in: number of variables, the size of x0. ! Program handles overdetermined, m>n, and underdetermined, m ! * See "Fortran90 for Fortran77 Programmers" by Clive Page ! 11/26/2001 ! DYNAMICALLY ALLOCATED: ! Jac scratch, m-by-n array: Jacobian of derivatives w.r.t. variables. ! H scratch, n-by-n array: JacT Jac Hessian without 2nd order derivs. ! L scratch, n-by-n array: lower-triangle factor, L LT = H+mu*I. ! r scratch, m vector of residuals. ChiSq=Sum(rT r). ! D,g,sN,s,x scratch, n vectors. ! They are scale factor, gradient, Newton step, step tried, position. ! You must scale by scaling factors D. ! Use dfj/dxi=Jac(j,i)/D(i) H(i,j)/D(i)/D(j) g(j)/D(j) sN(j)*D(j). ! INTERFACE: use GN_MOD implicit none external FCN,GN_CHOL,GN_LSOLVE intrinsic abs,epsilon,matmul,max,min intrinsic minval,sign,size,sqrt,sum,transpose integer info,m,n real*8 f0,GN_CauchyDist,x0(n) ! INTERNAL: logical hook,take,done integer i,iHOOK,j,nbad,nrank1,nconsec real*8 del,delFac,EPS,epsn12,f,sfp real*8 gnorm,hold,lastdel,mu,mulow,muhigh,muStart real*8 phi,phip,rr,scal,scal0,snewt,snorm,step,temp ! DYNAMIC: real*8, dimension(:,:), allocatable :: Jac,H,L real*8, dimension(:), allocatable :: D,g,r,s,sN,x if(min(limit,m,n).le.0) then write(*,*)'INPUT ERROR limit,m,n',limit,m,n info=128 return endif ! For Fortran90 or better, allocate module arrays and local arrays. ! To avoid fragmentation of the stack: allocate globals, then locals, ! deallocate locals in reverse order, then globals in reverse order. ! Globals remain allocated until program is run next or you do it. ! Deallocate/allocate global arrays in GN_MOD. ! with iscale=2, user needs to include the statements: ! use GN_MOD, only:D0 ! if(allocated(D0)) deallocate(D0) ! allocate(D0(n)) ! D0=... ! before calling GN. ! Scale vector D0 needs to have all n elements positive. if(iscale.eq.2) then if(.not.allocated(D0)) then info=128 write(*,*)'INPUT ERROR scale vector D0 not allocated' elseif(size(D0).lt.n.or.minval(D0).le.0) then info=128 write(*,*)'INPUT ERROR bad scale vector D0',D0 endif endif if(info.eq.128) return if(iscale.ne.2) then if(allocated(D0)) deallocate(D0) allocate(D0(n)) endif if(allocated(Jac0)) deallocate(Jac0) if(allocated(H0)) deallocate(H0) if(allocated(r0)) deallocate(r0) allocate(r0(m),H0(n,n),Jac0(m,n)) ! Allocate local arrays. allocate(Jac(m,n),H(n,n),L(n,n)) allocate(r(m),D(n),g(n),s(n),sN(n),x(n)) EPS=epsilon(ONE) epsn12=sqrt(n*EPS) x=x0 D=1 if(iscale.eq.2) D=D0(:n) nfcn=0 nrank1=n hook=.false. take=.false. nconsec=0 nconsecMax=0 H=0 if(ihist.gt.0) then open(ihist,file=HistFile,access="append") write(ihist,*)'iter,f,x' endif if(iprint.ge.4) then write(*,*) write(*,*)'GN: new Problem: m,n,NewtStep1st=',m,n,NewtStep1st write(*,*)'derivstp,stptol,reltol,abstol', 1 derivstp,stptol,reltol,abstol if(iscale.eq.2) write(*,*) 'fixed scale D based on user D0',D endif ! REPEAT UNTIL info.ne.0 and have a new Jacobian. ! Do not indent main loop to keep width smaller--------------------- MainLoop: do iter=1,limit ! termination on nfcn >= limit info=0 nfcn=nfcn+1 call FCN(m,n,x,r) f=sum(r**2) if(iter.gt.1) then take=f.lt.f0 if(take) then nconsec=nconsec+1 nconsecMax=max(nconsecMax,nconsec) else nconsec=0 endif endif if(iprint.ge.4) write(*,*) if(iprint.ge.2) write(*,*) 1 'after call FCN: iter,nfcn,f,info,take',iter,nfcn,f,info,take if(iprint.ge.3) write(*,*) 'x',x if(iprint.ge.4) write(*,*) 'r',r if(iter.eq.1) then take=.true. elseif(iter.eq.2.and.NewtStep1st.and..not.take) then ! CAUCHY STEP is distance to Cauchy point in Newton direction. ! Here if Newton step didn't work. hook=.false. del=GN_CauchyDist(n,L,g,gnorm) else ! UPDATE TRUST REGION. ! Strict descent requires f f0 delfac=0.5/max(1-rr,EPS) del=snorm*min(max(delFac,ZLOW),ZHIGH) ! we also increase trust radius based on number of consecutive takes. ! Test problem #36 illustrates the need for this. if(iprint.ge.4)then write(*,*)'Update trust region del,lastdel,rr,delFac' 1 ,del,lastdel,rr,delFac if(nconsec.gt.0) write(*,*)' nconsec, ZCONSEC**nconsec' 1 ,' Newdel=',nconsec,ZCONSEC**nconsec,del*ZCONSEC**nconsec endif del=del*ZCONSEC**nconsec endif if(take) then ! SAVE THE BEST. if(iprint.ge.4) write(*,*) 'saving best' nbad=0 x0=x f0=f mu0=mu del0=del D0=D H0=H Jac0=Jac if(ihist.gt.0) write(ihist,*) iter,f0,x0 elseif(iter.ne.2.or..not.NewtStep1st) then ! Overshot the mark. Try, then get a new Jacobian. hook=.true. if(nrank1.gt.0) nbad=nbad+1 if(iprint.ge.4) write(*,*)'Overshot nrank1,nbad',nrank1,nbad endif ! Some convergence tests. if(iter.gt.1) then if(snorm.le.stptol) info=info+1 endif if(f.lt.abstol) info=info+8 done=(info.ne.0.and.nrank1.eq.0).or.info.ge.8 if(.not.done.and.nfcn.ge.limit) info=130 if(iter.gt.1.and.(done.or.info.ge.128)) exit MainLoop ! ---Jacobian dfi/dxj, gradient g=JacT r, Hessian H=JacT Jac+2nd order ! Jacobian update when not stale nor final, else a full Jacobian. if(nbad.le.MXBAD.and.info.eq.0.and.nrank1.lt.n) then if(.not.take) goto 40 ! Rank-1 update to Jacobian. Jac(new) = Jac+((r-r0-Jac s) sT)/(sT s). if(iprint.ge.4) write(*,*) 1 'Rank=1 Jacobian update: nrank1,nbad,snorm',nrank1,nbad,snorm nrank1=nrank1+1 Jac=Jac+matmul(reshape(r-r0-matmul(Jac,s),(/m,1/)), 1 reshape(s/snorm**2,(/1,n/))) !outer product if(take) r0=r else ! Full Jacobian. ! Step away from zero to avoid crossing it. if(iprint.ge.4) write(*,*)'Full Jacobian: nrank1,nbad,take' 1 ,nrank1,nbad,take nrank1=0 nbad=0 if(take) r0=r nfcn=nfcn+n if(iscale.eq.1) then ! variable scale scal0=sum((s*D)**2) D=max(0.5*(abs(x0)+D),stptol) if(iprint.ge.4) write(*,*) 'variable scale D',D scal=sum((s*D)**2) if(scal.gt.0) del=del*sqrt(scal0/scal) endif do j=1,n hold=x0(j) step=sign(derivstp,hold) x0(j)=x0(j)+step*D(j) call FCN(m,n,x0,r) Jac(:,j)=(r-r0)/step x0(j)=hold enddo endif ! Gradient and Hessian. g=matmul(transpose(Jac),r0) H=matmul(transpose(Jac),Jac) if(iprint.ge.5) then do i=1,m write(*,*) 'J',Jac(i,:) enddo do j=1,n write(*,*) 'H',H(j,:) enddo endif gnorm=sqrt(sum(g**2)) if(iprint.ge.4) write(*,*) 'gnorm,g',gnorm,g ! L1 norm (max of row sums of abs) of H(i,j). ! H=JacT Jac symmetric. hnorm=0 do j=1,n hnorm=max(hnorm,sum(abs(H(:,j)))) enddo ! Get a small number for further augmentation, check underflow. hnorm=hnorm*epsn12 if(hnorm.eq.0) then info=129 exit MainLoop endif ! Find bad rows and ignore them. do j=1,n if(H(j,j).le.0) then H(j,:)=0 H(:,j)=0 H(j,j)=1 endif enddo ! --- GET NEWTON STEP, H sN = -g --- mu=0 ! solve (H + mu I)s = g for step s, possibly augmenting diagonal of H call GN_CHOL(n,mu,H,g,L,sN,hnorm) snewt=sqrt(sum(sN**2)) if(snewt.le.del) hook=.false. if(iprint.ge.4) write(*,*) 'snewt,mu,sN',snewt,mu,sN if(iter.eq.1) then if(NewtStep1st) then ! Try Newton step first. del=snewt s=sN snorm=snewt if(iprint.ge.4) write(*,*)'Taking Newton step first' goto 100 else ! more conservative approach, try Cauchy step first hook=.false. del=GN_CauchyDist(n,L,g,gnorm) endif endif 40 continue ! --- CHOOSE NEWTON OR HOOK STEPS ------------------- if(.not.hook) then ! Allow Newton step to be up to the current trust radius del temp=1 if(snewt.gt.del) temp=del/snewt s=sN*temp snorm=snewt*temp if(iprint.ge.4) write(*,*) 'Step in Newton direction:' 1 ,' snewt,del,snorm',snewt,del,snorm else ! --- Hook search --- ! Find step of length del by finding mu ! that gives snorm = ||s|| = ||(H + mu I)**-1 g|| = del. ! Because del = ||(H + mu I)**-1 g|| and H is positive, ! del is less than ||g||/mu. muhigh=gnorm/del mulow=0 ! mulow <= mu <= muhigh muStart=mu ! REPEAT UNTIL abs(del-||s||)<.05 or mulow>muhigh. HookLoop: do iHOOK=1,30 call GN_CHOL(n,mu,H,g,L,s,hnorm) mulow=max(mulow,mu) snorm=sqrt(sum(s**2)) if(abs(snorm-del).le.0.05d0*del.or.mulow.ge.muhigh) 1 exit HookLoop phi=snorm-del ! phi<0 indicates mu is too large, use this mu to update muhigh if(phi.lt.0) muhigh=min(mu,muhigh) if(phi.gt.0) mulow=max(mu,mulow) ! Want sT (H + mu I)**-1 s = ||L**-1 s||**2. ! Start with x = L**-1 s. call GN_LSOLVE(n,L,x,s) phip=0 if(snorm.gt.0) phip=-sum(x**2)/snorm ! As mu increases the step size decreases, so phip must be negative. if(phip.lt.0) mu=mu-(snorm/del)*(phi/phip) mu=max(min(mu,muhigh),mulow) enddo HookLoop ! End of REPEAT. if(iprint.ge.4.or.iHOOK.ge.30) then write(*,*) 'iHOOK,iHOOKmax,muStart,mu,lastdel,del,snorm' 1 ,iHOOK,iHOOKmax,muStart,mu,lastdel,del,snorm info = 133 if(iHOOK.ge.30) exit MainLoop endif iHOOKmax=max(iHOOKmax,iHOOK) endif 100 continue ! TAKE THE STEP. if(iprint.ge.4) write(*,*) 'Taking step s',s x=x0+s*D enddo MainLoop ! termination on nfcn >= limit if(iprint.ge.1) then write(*,'(" iter m n nf f info take",i4,2i3,i5 1 ,1pg21.14,i4,L2)')iter,m,n,nfcn,f,info,take write(*,*) 'x0',x0 endif ! memory freed in reverse order from borrowed as in next line. deallocate(x,sN,s,g,D,r,L,H,Jac) if(iprint.ge.4) then write(*,*)'GN: finished Problem' write(*,*) endif end subroutine GN !----------------------------------------------------------------------- !+Cauchy distance is distance along gradient to minimum. ! Cauchy step is Cauchy distance in the Newton direction. ! Cauchy distance = -||g||**2/(gT H g) ||g||, gT H g = ||LT g||**2. ! Use ZCP * Cauchy distance. ! ||v|| = sqrt(vT v) is L2 or Euclidean norm of vector v. real*8 function GN_CauchyDist(n,L,g,gnorm) use GN_MOD, only: ZCP,ZCPMIN,ZCPMAX,iprint implicit none intrinsic max,min,sum real*8 L(n,n),g(n),gnorm,del,temp integer j,n if(gnorm.eq.0) then del=ZCPMIN return endif temp=0 ! calculate temp = gT H g/||g||**2 = ||LT g||**2/||g||**2 do j=1,n temp=temp+(sum(L(j:,j)*g(j:))/gnorm)**2 enddo temp=ZCP*gnorm/temp del=max(ZCPMIN,min(temp,n*ZCPMAX)) if(iprint.ge.4) write(*,*) 1 'Cauchy step,del,stpmax',temp,del,n*ZCPMAX GN_CauchyDist=del end function GN_CauchyDist !----------------------------------------------------------------------- !+Cholesky decomposition and solution, (H + (mu+add)*I) s = -g. ! Decomposition to lower triangle L without addition if it works. ! Break as soon as OK. subroutine GN_CHOL(n,mu,H,g,L,s,add) implicit none intrinsic sqrt,min,sum integer n,i,iadd,j real*8 add,mu,tmp,H(n,n),g(n),L(n,n),s(n) loop1: do iadd=1,n do j=1,n do i=j,n L(i,j)=H(j,i)-sum(L(i,:j-1)*L(j,:j-1)) enddo tmp=L(j,j)+mu if(tmp.le.0) then mu=mu+add cycle loop1 endif tmp=sqrt(tmp) L(j,j)=tmp L(j+1:,j)=L(j+1:,j)/tmp enddo exit enddo loop1 ! forward row reduction and backward substitution do i=1,n s(i)=(-g(i)-sum(L(i,:i-1)*s(:i-1)))/L(i,i) enddo do i=n,1,-1 s(i)=(s(i)-sum(L(i+1:,i)*s(i+1:)))/L(i,i) enddo end subroutine GN_CHOL !----------------------------------------------------------------------- !+solve L x = y. x and y may be the same vector, L lower triangle. subroutine GN_LSOLVE(n,L,x,y) implicit none intrinsic sum integer n,i real*8 L(n,n),x(n),y(n) do i=1,n x(i)=(y(i)-sum(L(i,:i-1)*x(:i-1)))/L(i,i) enddo end subroutine GN_LSOLVE !###############################END GN################################## !----------------------------------------------------------------------- module test_MOD implicit none integer,save:: icase,icasefmax,icasexmax,icase1=1,icase2=37,jxmax 1 ,m,n,nfcn,nfcnsum,nfcntot real*8,save:: ferr,ferrmax 1 ,xans(30),xerr,xerrmax,xerrmax1,xstart(30) real*8,save :: ans end module test_MOD !----------------------------------------------------------------------- subroutine test_init use test_MOD use GN_MOD, only:iscale,D0 implicit none integer j real*8 del,dx,xi,xmax,xmin select case(icase) case(1) n=2 m=2 xstart(1)=-1.2d0 xstart(2)=1 ans=0 xans(1)=1 xans(2)=1 case(2) n=2 m=2 xstart(1)=0.5d0 xstart(2)=-2 ans=48.9842d0 xans(1)=11.413d0 xans(2)=-0.89680d0 case(3) n=2 m=2 xstart(1)=0 xstart(2)=1 ans=0 xans(1)=0.10982d-04 xans(2)=9.1061d0 case(4) n=2 m=3 xstart(1)=1 xstart(2)=1 ans=0 xans(1)=1d6 xans(2)=2d-6 case(5) n=2 m=3 xstart(1)=1 xstart(2)=1 ans=0 xans(1)=3 xans(2)=0.5d0 case(6) n=2 m=10 xstart(1)=.3d0 xstart(2)=.4d0 ans=124.362d0 xans=.257825d0 case(7) n=3 m=3 xstart(1)=-1 xstart(2)=0 xstart(3)=0 ans=0 xans(1)=1 xans(2)=0 xans(3)=0 case(8) n=3 m=15 xstart(1)=1 xstart(2)=1 xstart(3)=1 ans=8.21487d-3 xans(1)=0.82411d-01 xans(2)=1.1330d0 xans(3)=2.3437d0 case(9) n=3 m=15 xstart(1)=0.4d0 xstart(2)=1 xstart(3)=0 do j=1,n xstart(j)=0.1d0 enddo ans=1.12793d-8 xans(1)=0.39896d0 xans(2)=1 xans(3)=0 case(10) n=3 m=16 xstart(1)=.02d0 xstart(2)=4000 xstart(3)=250 ans=87.9458d0 xans(1)=0.56096d-02 xans(2)=6181.3d0 xans(3)=345.22d0 case(11) n=3 m=99 xstart(1)=5 xstart(2)=2.5d0 xstart(3)=.15d0 ans=0 xans(1)=50 xans(2)=25 xans(3)=1.5d0 case(12) n=3 m=9 xstart(1)=0 xstart(2)=10 xstart(3)=20 ans=0 xans(1)=1 xans(2)=10 xans(3)=1 case(13) n=4 m=4 xstart(1)=3 xstart(2)=-1 xstart(3)=0 xstart(4)=1 ans=0 xans(1)=0 xans(2)=0 xans(3)=0 xans(4)=0 case(14) n=4 m=6 xstart(1)=-3 xstart(2)=-1 xstart(3)=-3 xstart(4)=-1 ans=0 xans(1)=1 xans(2)=1 xans(3)=1 xans(4)=1 case(15) n=4 m=11 xstart(1)=.25d0 xstart(2)=.39d0 xstart(3)=.415d0 xstart(4)=.39d0 ans=3.075056d-4 xans(1)=0.193d0 xans(2)=0.191d0 xans(3)=0.123d0 xans(4)=0.136d0 case(16) n=4 m=20 xstart(1)=25 xstart(2)=5 xstart(3)=-5 xstart(4)=-1 ans=85822.2d0 xans(1)=-11.594d0 xans(2)=13.204d0 xans(3)=-0.40339d0 xans(4)=0.23608 case(17) n=5 m=33 xstart(1)=.5d0 xstart(2)=1.5d0 xstart(3)=-1 xstart(4)=.01d0 xstart(5)=.02d0 ans=5.46d-5 xans(1)=0.37541d0 xans(2)=1.9358d0 xans(3)=-1.4647d0 xans(4)=0.12868d-01 xans(5)=0.22123d-01 case(18) n=6 m=13 xstart(1)=1 xstart(2)=2 xstart(3)=3 xstart(4)=4 ! More had xstart(3)=xstart(4)=1 ! These new starting values, when changed slightly ! give same final values. xstart(5)=1 xstart(6)=1 ans=0 xans(1)=4 xans(2)=10 xans(3)=3 xans(4)=5 xans(5)=1 xans(6)=1 case(19) n=11 m=65 xstart(1)=1.3d0 xstart(2)=.65d0 xstart(3)=.65d0 xstart(4)=.7d0 xstart(5)=.6d0 xstart(6)=3 xstart(7)=5 xstart(8)=7 xstart(9)=2 xstart(10)=4.5d0 xstart(11)=5.5d0 ans=4.01d-2 xans(1)=1.3100d0 xans(2)=0.432d0 xans(3)=0.63366d0 xans(4)=0.59943d0 xans(5)=0.754d0 xans(6)=0.90429d0 xans(7)=1.3658d0 xans(8)=4.8237d0 xans(9)=2.3987d0 xans(10)=4.57d0 xans(11)=5.6753d0 case(20) n=9 m=31 xstart=0 ans=1.39976d-6 xans(1)=-0.15328d-04 xans(2)=0.99979d0 xans(3)=0.14768E-01 xans(4)=0.14632d0 xans(5)=1.0009d0 xans(6)=-2.6177d0 xans(7)=4.1043d0 xans(8)=-3.1435d0 xans(9)=1.0526d0 case(21) n=12 m=12 do j=1,6 xstart(2*j-1)=-1.2d0 xstart(2*j)=1 enddo ans=0 do j=1,n xans(j)=1 enddo case(22) n=12 m=12 do j=1,3 xstart(4*j-3)=3 xstart(4*j-2)=-1 xstart(4*j-1)=0 xstart(4*j)=1 enddo ans=0 do j=1,n xans(j)=0 enddo case(23) n=4 m=5 do j=1,n xstart(j)=j xans(j)=0.25d0 enddo ans=2.2499775d-5 case(24) n=4 m=8 xstart=0.5d0 ans=9.376293d-6 xans(1)=0.2d0 xans(2)=0.191d0 xans(3)=0.480d0 xans(4)=0.519d0 case(25) n=9 m=11 do j=1,n xstart(j)=1-j/9.d0 xans(j)=1 enddo ans=0 case(26) n=9 m=n ! More had xstart(j)=1/n ! These new starting values, when changed slightly ! give same final values. xstart=0 xstart(1)=1 ans=0 xans(1)=0.56477q-01 xans(2)=0.60351q-01 xans(3)=0.58282q-01 xans(4)=0.62772q-01 xans(5)=0.65681q-01 xans(6)=0.69317q-01 xans(7)=0.20963q0 xans(8)=0.16722q0 xans(9)=0.12150q0 case(27) n=9 m=9 do j=1,n xstart(j)=0.5d0 xans(j)=1 enddo ans=0 case(28) n=9 m=9 do j=1,n xi=j xstart(j)=xi/10*(xi/10-1) enddo ans=0 xans(1)=-0.472d-1 xans(2)=-0.886d-1 xans(3)=-0.123d0 xans(4)=-0.149d0 xans(5)=-0.166d0 xans(6)=-0.171d0 xans(7)=-0.161d0 xans(8)=-0.133d0 xans(9)=-0.814d-1 case(29) n=9 m=9 do j=1,n xi=j xstart(j)=xi/10*(xi/10-1) enddo ans=0 xans(1)=-0.472d-1 xans(2)=-0.886d-1 xans(3)=-0.123d0 xans(4)=-0.149d0 xans(5)=-0.166d0 xans(6)=-0.171d0 xans(7)=-0.161d0 xans(8)=-0.133d0 xans(9)=-0.814d-1 case(30) n=9 m=9 do j=1,n xstart(j)=-1 enddo ans=0 xans(1)=-0.571d0 xans(2)=-0.682d0 xans(3)=-0.702d0 xans(4)=-0.704d0 xans(5)=-0.701d0 xans(6)=-0.692d0 xans(7)=-0.666d0 xans(8)=-0.596d0 xans(9)=-0.416d0 case(31) n=9 m=9 do j=1,n xstart(j)=-(-1)**j enddo ans=0 xans(1)=-0.427d0 xans(2)=-0.428d0 xans(3)=-0.477d0 xans(4)=-0.52d0 xans(5)=-0.558d0 xans(6)=-0.59248d0 xans(7)=-0.625d0 xans(8)=-0.592d0 xans(9)=-0.591d0 case(32) n=9 m=12 do j=1,n xstart(j)=1 xans(j)=-1 enddo ans=m-n case(33) n=9 m=12 do j=1,n xstart(j)=1 enddo ans=m*(m-1)/(4*m+2.d0) case(34) n=9 m=12 do j=1,n xstart(j)=1 enddo ans=(m**2+3*m-6)/(4*m-6.d0) case(35) n=12 m=9 ! More had n<=m del=1.d0/n xmin=-1+del xmax=1-del dx=(xmax-xmin)/(n-1.d0) do j=1,n xstart(j)=xmin+(j-1)*dx enddo ! More had xstart(j)=j/(n+1) ! This problem illustrates n>m ! also, starting parameter values can be changed slightly ! without affecting final values. ans=0 xans(1)=-0.932d0 xans(2)=-0.728d0 xans(3)=-0.6d0 xans(4)=-0.421d0 xans(5)=-0.24088d0 xans(6)=-0.78642d-01 xans(7)=-xans(6) xans(8)=-xans(5) xans(9)=-xans(4) xans(10)=-xans(3) xans(11)=-xans(2) xans(12)=-xans(1) case(36) ! Modified Beale (problem 5) that required many evaluations n=2 m=3 xans(1)=9833.5 xans(2)=0.99984d0 xstart(1)=3 xstart(2)=0.5d0 ans=1.232674 case(37) ! Rosenblock (problem 1) changed to single function (m=1) n=2 m=1 xstart(1)=-1.2d0 xstart(2)=1 ans=0 xans(1)=1 xans(2)=1 end select end !----------------------------------------------------------------------- subroutine test_fcn(m,n,x,r) use test_MOD, only: icase,nfcn implicit none integer i,j,j2,j4,m,n real*8 a,f1,f2,h,r(m),test_ex,test_cheb,theta,ti,total,ui,vi 1 ,wi,x(n),yi real*8 bard(15) data bard 1 /.14d0,.18d0,.22d0,.25d0,.29d0,.32d0,.35d0, 1 .39d0,.37d0,.58d0,.73d0,.96d0,1.34d0,2.10d0,4.39d0/ real*8 gaussi(15) data gaussi 1 /.0009d0,.0044d0,.0175d0,.0540d0,.1295d0 1 ,.2420d0,.3521d0,.3989d0,.3521d0,.2420d0,.1295d0 1 ,.0540d0,.0175d0,.0044d0,.0009d0/ integer meyer(16) data meyer 1 /34780,28610,23650,19630,16370,13720,11540,9744, 1 8261,7030,6005,5147,4427,3820,3307,2872/ real*8 kowyi(11) data kowyi 1 /.1957d0,.1947d0,.1735d0,.1600d0,.0844d0, 1 .0627d0,.0456d0,.0342d0,.0323d0,.0235d0,.0246d0/ real*8 kowui(11) data kowui 1 /4.d0,2.d0,1.d0,.5d0,.25d0, 1 .1670d0,.1250d0,.1000d0,.0833d0,.0714d0,.0625d0/ real*8 osbor1(33) data osbor1 1 /.844d0,.908d0,.932d0,.936d0,.925d0,.908d0,.881d0,.850d0, 1 .818d0,.784d0,.751d0,.718d0,.685d0,.658d0,.628d0,.603d0, 1 .580d0,.558d0,.538d0,.522d0,.506d0,.490d0,.478d0,.467d0, 1 .457d0,.448d0,.438d0,.431d0,.424d0,.420d0,.414d0,.411d0,.406d0/ real*8 osbor2(65) data osbor2 1 /1.366d0,1.191d0,1.112d0,1.013d0,.991d0, 1 .885d0,.831d0,.847d0,.786d0,.725d0, 1 .746d0,.679d0,.608d0,.655d0,.616d0, 1 .606d0,.602d0,.626d0,.651d0,.724d0, 1 .649d0,.649d0,.694d0,.644d0,.624d0, 1 .661d0,.612d0,.558d0,.533d0,.495d0, 1 .500d0,.423d0,.395d0,.375d0,.372d0, 1 .391d0,.396d0,.405d0,.428d0,.429d0, 1 .523d0,.562d0,.607d0,.653d0,.672d0, 1 .708d0,.633d0,.668d0,.645d0,.632d0, 1 .591d0,.559d0,.597d0,.625d0,.739d0, 1 .710d0,.729d0,.720d0,.636d0,.581d0, 1 .428d0,.292d0,.162d0,.098d0,.054d0/ real*8, save :: ONE=1 real*8 pi select case(icase) case(1) ! Rosenblock n=m=2 r(1)= 10*(x(2)-x(1)**2) r(2)= 1-x(1) case(2) ! Freudenstein and Roth n=m=2 r(1)=-13+x(1)+((5-x(2))*x(2)-2)*x(2) r(2)=-29+x(1)+((x(2)+1)*x(2)-14)*x(2) case(3) ! Powell badly scaled n=m=2 r(1) = 10000*x(1)*x(2)-1 r(2) = test_ex(-x(1))+test_ex(-x(2))-1.0001q0 case(4) ! Brown badly scaled n=2 m=3 r(1) = x(1)-1000000 r(2) = x(2)-2*ONE/10**6 r(3) = x(1)*x(2)-2 case(5) ! Beale n=2 m=3 r(1) = 1.5q0-x(1)*(1-x(2)) r(2) = 2.25q0-x(1)*(1-x(2)**2) r(3) = 2.625q0-x(1)*(1-x(2)**3) case(6) ! Jennrick and Sampson n=2 do i=1,m r(i)=2+2*i-test_ex(i*x(1))-test_ex(i*x(2)) enddo case(7) ! Helical valley n=m=3 pi=4*atan(ONE) theta=atan2(x(2),x(1))/(2*pi) r(1) = 10*(x(3)-10*theta) r(2) = 10*(sqrt(x(1)**2+x(2)**2)-1) r(3) = x(3) case(8) ! Bard n=3 m=15 do i=1,m ui = i vi = 16-i if(ui.lt.vi) then wi=ui else wi=vi endif r(i) = bard(i)-(x(1)+ui/(vi*x(2)+wi*x(3))) enddo case(9) ! Gaussian n=3 m=15 do i=1,m ti=ONE*(8-i)/2 r(i) = x(1)*test_ex(-x(2)*(ti-x(3))**2/2)-gaussi(i) enddo case(10) ! Meyer n=3 m=16 do i=1,m ti=45+5*i r(i)= x(1)*test_ex(x(2)/(ti+x(3)))-meyer(i) enddo case(11) ! Gulf research and development n=3 m<101 do i=1,m ti=i*ONE/100 yi=abs((-50*log(ti))**(2*ONE/3)+25-x(2)) r(i)=test_ex(-yi**x(3)/x(1))-ti enddo case(12) ! Box three dimensional n=3 do i=1,m ti=ONE*i/10 r(i)= test_ex(-ti*x(1))-test_ex(-ti*x(2)) 1 -x(3)*(test_ex(-ti)-test_ex(-10*ti)) enddo case(13) ! Powell singular n=m=4 r(1)=x(1)+10*x(2) r(2)=sqrt(5*ONE)*(x(3)-x(4)) r(3)=(x(2)-2*x(3))**2 r(4)=sqrt(10*ONE)*(x(1)-x(4))**2 case(14) ! Wood n=4 m=6 r(1)=10.d0*(x(2)-x(1)**2) r(2)=1-x(1) r(3)=sqrt(90*ONE)*(x(4)-x(3)**2) r(4)=1-x(3) r(5)=sqrt(10*ONE)*(x(2)+x(4)-2) r(6)=(x(2)-x(4))/sqrt(10*ONE) case(15) ! Kowalik and Osborne n=4 m=11 do i=1,m r(i) = kowyi(i) 1 -x(1)*(kowui(i)**2+kowui(i)*x(2)) 1 /(kowui(i)**2+kowui(i)*x(3)+x(4)) enddo case(16) ! Brown and Dennis n=4 do i=1,m ti=ONE*i/5 r(i) = (x(1)+ti*x(2)-test_ex(ti))**2 1 +(x(3)+x(4)*sin(ti)-cos(ti))**2 enddo case(17) ! Osborne 1 n=5 m=33 do i=1,m ti=10*ONE*(i-1) r(i) = osbor1(i)-(x(1)+x(2)*test_ex(-ti*x(4)) 1 +x(3)*test_ex(-ti*x(5))) enddo case(18) ! Biggs EXP6 n=6 do i=1,m ti=ONE*i/10 yi = test_ex(-ti)-5*test_ex(-10*ti)+3*test_ex(-4*ti) r(i) = x(3)*test_ex(-ti*x(1))-x(4)*test_ex(-ti*x(2)) 1 +x(6)*test_ex(-ti*x(5))-yi enddo case(19) ! Osborne 2 n=11 m=65 do i=1,m ti=(i-1)*ONE/10 r(i) = osbor2(i)-(x(1)*test_ex(-ti*x(5)) 1 +x(2)*test_ex(-(ti-x(9))**2*x(6)) 1 +x(3)*test_ex(-(ti-x(10))**2*x(7)) 1 +x(4)*test_ex(-(ti-x(11))**2*x(8))) enddo case(20) ! Watson n>1 m=31, n=9 r(30)=x(1) r(31)=x(2)-x(1)**2-1 do i=1,29 yi=0 ui=x(1) ti=i*ONE/29 do j=2,n ui=ui+x(j)*ti**(j-1) yi=yi+(j-1)*x(j)*ti**(j-2) enddo r(i)=yi-ui**2-1 enddo case(21) ! Extended Rosenblock n=m=2k do j2=1,6 r(2*j2-1)=10*(x(2*j2)-x(2*j2-1)**2) r(2*j2)=1-x(2*j2-1) enddo case(22) ! Extended Powell singular n=m=4k do j4=1,m/4 r(4*j4-3)=x(4*j4-3)+10*x(4*j4-2) r(4*j4-2)=sqrt(5*ONE)*(x(4*j4-1)-x(4*j4)) r(4*j4-1)=(x(4*j4-2)-2*x(4*j4-1))**2 r(4*j4)=sqrt(10*ONE)*(x(4*j4-3)-x(4*j4)) enddo case(23) ! Penalty I m=n+1 total=0 do j=1,n r(j)=sqrt(ONE/10**5)*(x(j)-1) total=total+x(j)**2 enddo r(m)=total-ONE/4 case(24) ! Penalty II m=2n r(1)=x(1)-2*ONE/10 total=0 do j=1,n total=total+(n-j+1)*x(j)**2 enddo r(2)=total-1 a=sqrt(ONE/10**5) do j=2,n yi = test_ex(j*ONE/10)+test_ex((j-1)*ONE/10) f1 = a*(test_ex(x(j)/10)+test_ex(x(j-1)/10)-yi) f2 = a*(test_ex(x(j)/10)-test_ex(-ONE/10)) r(2*j-1)=f1 r(2*j)=f2 enddo case(25) ! Variably dimensioned m=n+2 do j=1,n r(j)=x(j)-1 enddo total=0 do j=1,n total=total+j*(x(j)-1) enddo r(n+1)=total r(n+2)=total**2 case(26) ! Trigonmetric m=n total=0 do j=1,n total=total+cos(x(j)) enddo do i=1,m r(i)=n-total+i*(1-cos(x(i)))-sin(x(i)) enddo case(27) ! Brown almost-linear m=n total=0 do j=1,n total=total+x(j) enddo do j=1,n-1 r(j)=total-(n+1)+x(j) enddo r(n)=x(1)-1 case(28) ! Discrete boundary value m=n h=ONE/(n+1) do j=1,n r(j)=2*x(j)+h**2*(x(j)+j*h+1)**3/2 if(j.gt.1) r(j)=r(j)-x(j-1) if(j.lt.n) r(j)=r(j)-x(j+1) enddo case(29) ! Discrete integral equation m=n h=ONE/(n+1) yi=0 do j=1,n ui=0 do i=j+1,n ti=i*h ui=ui+(1-ti)*(x(i)+ti+1)**3 enddo ti=j*h yi=yi+ti*(x(j)+ti+1)**3 r(j)=x(j)+h*((1-ti)*yi+ti*ui)/2 enddo case(30) ! Broyden tridiagonal m=n do j=1,n r(j)=(3-2*x(j))*x(j)+1 if(j.gt.1) r(j)=r(j)-x(j-1) if(j.lt.n) r(j)=r(j)-2*x(j+1) enddo case(31) ! Broyden banded m=n do j=1,n yi=0 do i=j-5,j+1 if(i.gt.1.and.i.lt.n) then if(i.ne.j) yi=yi+x(i)*(1+x(i)) endif enddo r(j)=x(j)*(2+5*x(j)**2)+1-yi enddo case(32) ! Linear function--full rank total=0 do j=1,n total=total+x(j) enddo yi=2*total/m do j=1,n r(j)=x(j)-yi-1 enddo do i=n+1,m r(i)=-yi-1 enddo case(33) ! Linear--rank 1 total=0 do j=1,n total=total+x(j)*j enddo do i=1,m r(i)=i*total-1 enddo case(34) ! Linear function--rank 1 with zero columns and rows total=0 do j=1,n-2 total=total+x(j)*j enddo do i=1,m-2 r(i)=i*total-1 enddo r(m-1)=-1 r(m)=-1 case(35) ! Chebyquad 1/n Sum j=1,n Ti(x(j)) - Int 0 to 1 Ti(x) dx = 0 do i=1,m yi=0 do j=1,n yi=yi+test_cheb(i,x(j)) enddo r(i)=yi/n if(mod(i,2).eq.0) r(i)=r(i)+ONE/(i**2-1) enddo case(36) ! Modified Beale (problem 5) that required many evaluations r(1) = ONE/2-x(1)*(1-x(2)) r(2) = 3-x(1)*(1-x(2)**2) r(3) = 5-x(1)*(1-x(2)**3) case(37) ! Rosenblock (problem 1) changed to single residual (m=1) r(1)=sqrt((10*(x(2)-x(1)**2))**2+(1-x(1))**2) end select nfcn=nfcn+1 end !----------------------------------------------------------------------- real*8 function test_ex(x) ! substitute exponential to avoid overflow implicit none intrinsic exp real*8 x ! in order to prevent single precision overflow if(x.lt.87d0) then test_ex=exp(x) return endif test_ex=exp(87d0) end !----------------------------------------------------------------------- real*8 function test_cheb(i,x) ! substitute chebychev function with penalty for abs(x) > 1 implicit none intrinsic acos,cos integer i real*8 x if(abs(x).gt.1) then test_cheb=1000000 return endif test_cheb=cos(i*acos(x)) end !----------------------------------------------------------------------- character*(*) function GN_infoTxt(info) ! decoded info returned as text function implicit none integer info,j intrinsic btest character*100 txt1 logical btest character*25 mess(8) data mess/ 1 "Step < stptol" 1 ,"RelativeImprov < reltol" 1 ," " 1 ,"F < abstol" 1 ," " 1 ,"H sing (no linear uncert)" 1 ," " 1 ,"Input error"/ ! 1 step less than tolerance (stptol). ! 2 relative improvement less than tolerance (reltol). ! 4 ! 8 sum of squares less than absolute tolerance (abstol). ! 32 H singular, linear uncertainty calculation not possible. ! 64 ! 128 input error. ! 129 Jacobian zero. ! 130 evaluation limit exceeded. ! 133 hook search error. write(GN_infoTxt,'(i3,":")') info txt1='' if(info.eq.0) then txt1='no fit' elseif(info.eq.129) then txt1='Jacobian zero' elseif(info.eq.130) then txt1='evaluation limit exceeded' elseif(info.eq.133) then txt1='hook search error' else do j=1,8 if(btest(info,j-1)) then if(txt1.eq.'') then txt1=trim(mess(j)) else txt1=trim(txt1)//';'//mess(j) endif endif enddo endif GN_infoTxt=trim(GN_infoTxt)//trim(txt1) END !----------------------------------------------------------------------- subroutine test ! Test GN using More' test problems. use GN_MOD, only: abstol,D0,del0,derivstp 1 ,HistFile,ihist,iHOOKmax,iprint,iscale,iter 1 ,limit,mu0,MXBAD,nconsecMax,NewtStep1st,reltol,stptol 1 ,ZCONSEC,ZCP,ZCPMIN,ZHIGH,ZLOW use test_MOD implicit none character*100 GN_infoTxt real*8 csq,ftol,sum,x(30),xtol integer info,irun,j external test_fcn data xtol/0.01d0/,ftol/5d-5/ ! changes in default values from GN_MOD limit=1000000 ! changes in default values from test_MOD icase1=1 icase2=36 nfcnsum=0 do irun=1,3 write(*,*) if(irun.eq.1) then write(*,*)'START NEAR ANSWER' elseif(irun.eq.2) then write(*,*)'NORMAL START' elseif(irun.eq.3) then write(*,*)'PERTURBED NORMAL START' endif write(*,*)'case m n nconsec iter nfcn f' 1 ,' mu del f X MaxXerr info' nfcntot=0 xerrmax=0 ferrmax=0 do icase=icase1,icase2 call test_init if(irun.eq.1) xstart=xans-0.001q0 ! start near answer x=xstart ! normal starting point for problem icase if(irun.eq.3) x=xstart+.04q0 ! perturbed normal start nfcn=0 call GN(test_fcn,m,n,info,x,csq) select case(icase) case(33) ! n-1 dimension continuum of answers sum=0 do j=1,n-1 xans(j)=x(j) sum=sum+j*x(j) enddo xans(n)=(3/(2*m+1d0)-sum)/n case(34) ! n-1 dimension continuum of answers xans(n)=x(n) xans(n-1)=x(n-1) sum=0 do j=1,n-3 xans(j)=x(j) sum=sum+j*x(j) enddo xans(n-2)=(3/(2*m-3d0)-sum)/(n-2) case(36) ! this problem is approximately dependent only on x(1)*(1-x(2)) xans(2)=x(2) xans(1)=1.54d0/(1-x(2)) end select ! update maximum parameter errors xerrmax1=0 do j=1,n if(xans(j).gt.0) then xerr=min(abs(x(j)-xans(j)) 1 ,abs(x(j)-xans(j))/abs(xans(j))) else xerr=abs(x(j)-xans(j)) endif if(xerr.gt.xerrmax) then icasexmax=icase jxmax=j xerrmax=xerr endif if(xerr.gt.xerrmax1) then xerrmax1=xerr endif enddo ! update maximum function error if(ans.gt.0) then ferr=min(csq-ans,(csq-ans)/ans) else ferr=csq-ans endif if(ferr.gt.ferrmax) then icasefmax=icase ferrmax=ferr endif write(*,'(i5,i3,i3,3i8,g22.14,2(1pe10.3), 1 1x,l1,1x,l1,1x,1pe9.2,a100)') 1 icase,m,n,nconsecMax,iter,nfcn,csq,mu0,del0 1 ,ferr.lt.ftol,xerrmax1.lt.xtol,xerrmax1,GN_infoTxt(info) nfcntot=nfcntot+nfcn enddo if(xerrmax.gt.xtol.or.ferr.gt.ftol) nfcntot=nfcntot+1000000 ! number of millions indicates number of problems with bad results write(*,*)'nfcntot=',nfcntot nfcnsum=nfcnsum+nfcntot write(*,*)'maximum parameter error=',xerrmax 1 ,' for icase=',icasexmax,' j=',jxmax write(*,*)'maximum function error=',ferrmax 1 ,' for icase=',icasefmax write(*,*)'derivstp,stptol,limit,abstol,reltol=' 1 ,derivstp,stptol,limit,abstol,reltol write(*,*)'ZLOW,ZHIGH,ZCP,ZCPMIN,ZCONSEC,MXBAD=' 1 ,ZLOW,ZHIGH,ZCP,ZCPMIN,ZCONSEC,MXBAD write(*,*)'iscale,NewtStep1st,iHOOKmax=' 1 ,iscale,NewtStep1st,iHOOKmax enddo write(*,*)'nfcnsum=',nfcnsum write(*,*)'ftol,xtol=',ftol,xtol end !-----------------------------------------------------------------------