C  TEST ALG. 16 GRID SEARCH  J.C.NASH  JULY 1978
      REAL VAL(20)
      EXTERNAL FUNS 
      INTEGER IPR,COUNT,NBIS,NIN
      LOGICAL NOCOM 
      REAL U,V,TOL
      NIN=15
      IPR=16
   1  READ(NIN,901)COUNT,NBIS,U,V,TOL
 901  FORMAT(2I5,3F10.5)
      WRITE(IPR,910)COUNT,NBIS,U,V,TOL
 910  FORMAT('0TEST- COUNT=',I5,' NBIS=',I5,' U=',F15.5,' V=',F15.5,
     *' TOL=',F15.10)
      IF(COUNT.LE.0)STOP
      CALL A16GS(U,V,NBIS,FUNS,COUNT,TOL,IPR,T,VAL)
      WRITE(IPR,963)U,V,T,COUNT
 963  FORMAT('0FINAL INTERVAL=(',1PE16.8,',',E16.8,')  LOWEST VALUE=',
     #E16.8,' COUNT=',I4)
      GOTO 1
      END 
      REAL FUNCTION FUNS(X,NOCOM)
      LOGICAL NOCOM 
      REAL X
      NOCOM=.FALSE. 
      IF(X.LT.-5.0.OR.X.GT.10.0)NOCOM=.TRUE.
C  DEFINE RANGE OF INTEREST
      IF(NOCOM)RETURN
      FUNS=X*(X*X-2.0)-5.0
      RETURN
      END 
      SUBROUTINE A16GS(U,V,N,FNS,IFN,TOL,IPR,T,VAL)
C  ALGORITHM 16  GRID SEARCH
C   J.C.NASH  FEBRUARY 1980
C  U,V DEFINE THE INTERVAL OF INTEREST
C  N  GIVES THE NUMBER OF DIVISIONS (N+1) POINTS
C  FNS IS THE NAME OF THE FUNCTION    VAL=FNS(B,NOCOM)
C     NOCOM SET .TRUE. IF NOT COMPUTABLE. PROGRAM HALTS IN THIS CASE
C  IFN  IS LIMIT ON FUNCTION EVALUATIONS ALLOWED. RETURNS ACTUAL USED 
C  TOL =CONVERGENCE TOLERANCE ON ABS(V-U)*2/N
C  IPR = PRINT CHANNEL   IPR.GT.0 FOR PRINTING.
C  T  = LOWEST VALUE FOUND
C  VAL    =  WORKING VECTOR OF VALUES AT GRID POINTS
C  STEP 0 
      LOGICAL NOCOM 
      INTEGER N,K,J,LIM,N1
      REAL H,U,V,T,TOL,P,SV,X,VAL(N)
C  N.LE.2  CAN'T REDUCE INTERVAL
      IF(N.LT.3)STOP
      LIM=IFN
      IFN=0
      NOCOM = .FALSE.
      T=FNS(U,NOCOM)
      IF(NOCOM)STOP 
      IFN=IFN+1
      IF(IPR.GT.0)WRITE(IPR,956)IFN,U,T 
      VAL(1)=T
      SV=FNS(V,NOCOM)
      IF(NOCOM)STOP 
      IFN=IFN+1
      IF(IPR.GT.0)WRITE(IPR,956)IFN,V,SV
C  STEP 1 
  10  K=0 
      IF(SV.GE.T)GOTO 15
      K=N 
      T=SV
  15  H=(V-U)/N
C  STEP 2 
C  S(U) ALREAD IN T 
      N1=N-1
      DO 60 J=1,N1
C  STEP 3 
        X=U+J*H
        P=FNS(X,NOCOM)
        IF(NOCOM)STOP
        IFN=IFN+1
        IF(IFN.GE.LIM)RETURN
      IF(IPR.GT.0)WRITE(IPR,956)IFN,X,P 
 956  FORMAT( 8H EVALN #,I4,4H  F(,1PE16.8,2H)=,E16.8)
C  SAVE VALUE
      VAL(J+1)=P
C  STEP 4 
        IF(P.GE.T)GOTO 60
C  STEP 5 
        T=P
        K=J
C  STEP 6 
  60  CONTINUE
C  STEP 7 
      IF(ABS(H).LT.0.5*TOL)RETURN
C  STEP 8 
      V=U+(K+1)*H
      U=V-2*H
      IF(K.EQ.0)GOTO 82
C  S(U) IS IN VAL(K)
      T=VAL(K)
      GOTO 84
  82  T=FNS(U,NOCOM)
      IF(NOCOM)STOP 
      IFN=IFN+1
      IF(IPR.GT.0)WRITE(IPR,956)IFN,U,T 
      IF(IFN.GE.LIM)RETURN
  84  IF(K.GT.N-2)GOTO 86
      SV=VAL(K+2)
      GOTO 10
  86  IF(K.EQ.N1)GOTO 10
C  SV  ALREADY IN PLACE IF K=N-1
      SV=FNS(V,NOCOM)
      IF(NOCOM)STOP 
      IFN=IFN+1
      IF(IPR.GT.0)WRITE(IPR,956)IFN,V,SV
      IF(IFN.GE.LIM)RETURN
      GOTO 10
      END 
   80    5   0.0        3.0     0.00001 
   80    5   0.0        3.0     0.001
   80    5   0.0        3.0     0.01
   80    5   0.0        3.0     0.0
   80   10   0.0        3.0     0.00001 
   80   10   0.0        3.0     0.001
   80   10   0.0        3.0     0.01
   80   10   0.0        3.0     0.0
    5   10   0.0        3.0     0.0
    0    0