c     
c     Second Edition
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1996
c     ISBN 0-534-33892-5
c     COPYRIGHT (c) 1996
c
c     Section 6.13
c
c     Example of Adaptive approximation
c 
c
c     file: adapta.f
c
      parameter (m=100) 
      dimension t(0:m),y(0:m),c(1:m),d(0:m)
      data a,b /0.0,1.0/    
      external f
c
      print *
      print *,' Adaptive approximation example'
      print *,' Section 6.13, Kincaid-Cheney'
      print *
c
      e = 1.0
      d(0) = 0.0
      do 7 icase=1,3
         e = e/10.   
         print *,' e =',e
         t(0) = a
         t(1) = b
         y(0) = f(t(0))
         y(1) = f(t(1))
         call maximum(f,t(0),t(1),c(1),d(1)) 
         do 4 n=1,m-1
            v = 0.0 
            do 2 j=1,n
               if (d(j) .gt. v) then
                  v = d(j)  
                  i = j 
               endif
 2          continue
            if (v .le. e) goto 5 
            do 3 j=n,i+1,-1 
               t(j+1) = t(j)
               y(j+1) = y(j)
               d(j+1) = d(j)
               c(j+1) = c(j)
 3          continue  
            t(i+1) = t(i)
            y(i+1) = y(i)
            t(i) = c(i) 
            y(i) = f(t(i))
            call maximum(f,t(i-1),t(i),c(i),d(i))       
            call maximum(f,t(i),t(i+1),c(i+1),d(i+1)) 
 4    continue
c
 5    print *,' n =',n
      print *
      print *,'    knots        value of f     deviations'
c
      do 6 i = 0,n 
         print 8,t(i),y(i),d(i) 
 6    continue  
      print *
 7    continue
c
 8    format (2x,3(e13.6,2x))
      stop
      end 
c  
      function f(x) 
      f = sqrt(x)
      return
      end 
c  
      subroutine maximum(f,a,b,c,d) 
c 
c     determine the maximum
c
      external f
      parameter (k=10) 
      dimension z(0:k)
c
      h = (b-a)/real(k) 
      do 2 i=0,k
         z(i) = f(a+h*real(i)) 
 2    continue  
      do 3 i=1,k-1
        z(i) = abs(z(i)-(z(k)*real(i)+z(0)*real(k-i))/real(k))
 3    continue  
      d = 0.0     
      do 4 i = 1,k-1
        if (z(i) .gt. d) then   
           d = z(i)
           c = a + h*real(i) 
        end if
 4    continue
c
      return
      end 
