       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<n,
!    as well as linear and nonlinear equations (m=n).
! info   out: status sum of bits (LT 128) or single error code (GE 128)
!    Notation: pT q=transpose of p times q, ||q||=norm(q)=sqrt(qT q).
! CONVERGENCES
!     1 step less than tolerance (stptol).
!     2 rT r relative improvement less than tolerance (reltol).
!     8 rT r sum of squares less than absolute tolerance (abstol).
! FATAL ERRORS
!    128 input error.
!    129 Jacobian is zero--cannot proceed with unchanging function.
!    130 evaluation limit exceeded without termination.
!    133 hook search error: could not find a suitable mu in 30 tries.
! x0     in/out n-vector: initial guess and best position found.
! f0     out scalar: final sum of squares. Think chi-squared.
! ACCESSIBLE in GN_MOD:
! r0     out m-vector: final residuals at x0
! iter   out scalar: actual number of iterations taken
! nfcn   out scalar: actual number of function evaluations used
! mu0    out scalar: the "best" L-M parameter used
! del0   out scalar: the "best" trust region size
! D0     out n-vector: scale factor, input for iscale=2
! H0     out n-by-n matrix: the unaugmented Hessian
! Jac0   out m-by-n array: the "best" Jacobian
! LIMITATIONS:
!    Strict descent limits use to reproducible data.
!    x0, r0, and others must be 1-D vectors, not arrays.
!    Large m and n will require large storage and round-off errors.
!    Your residuals should not exceed sqrt(HUGE/m) to avoid overflows.
! RECOMMENDATIONS:
!    Use double-precision calculation, REAL*8, if possible.
!    Single precision can do 3 to 4 digits in x, at best.
!    Give a good guess or you may find a secondary minimum.
!    Scale parameters and residuals, so chi-square and relative
!    errors are meaningful.
! CREDITS:
!     * Dennis and Schnabel, "Numerical Methods for Unconstrained
!     Optimization and Nonlinear Equations" (1983),
!     for many ideas and the hook calculation, i.e., mu from phi.
!     * IMSL ZXSSQ for the Jacobian update attributed in
!     Fletcher, "Practical Methods of Optimization" (1980)
!     to Barnes (1965) and Broyden (1965).
!     * x(+) = x-(JacT Jac+mu I)**-1 JacT r augmented Gauss-Newton.
!     * See this for allocations and vector operators in F90.
!     <http://www.stanford.edu/class/me200c/tutorial_90/07_arrays.html>
!     * See "Fortran90 for Fortran77 Programmers" by Clive Page
!     <http://www.star.le.ac.uk/~cgp/f90course/f90.html> 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.
! Parabolic interpolation/extrapolation using f0, f and slope at x0.
! Update del to be distance to minimum, new del is delFac times old.
         lastdel=del
         if(abs(f-f0).le.max(reltol*max(f,f0),abstol)) info=2
! sfp is the step length s times the slope fp of f in the direction s.
! it would be the change in f for negligible curvature.
! sfp is negative because the step is in a decreasing direction.
         sfp=2*sum(g*s)
         rr=(f-f0)/sfp
! rr is bounded above by 1, when f < f0 and f behaves linearly,
! and unbounded below when 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
!-----------------------------------------------------------------------