c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
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(1:m)
      data a,b /0.0,1.0/
c
      print *
      print *,' Adaptive approximation example'
      print *,' Section 6.13, Kincaid-Cheney'
      print *
c
      e = 1.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

From cs.utexas.edu!kincaid Thu May 27 02:02:40 0500 1993
From: kincaid@cs.utexas.edu (David R. Kincaid)
Date: Thu, 27 May 1993 02:02:40 -0500
To: ehg@research.att.com
Subject: cg.f
Cc: kincaid@cs.utexas.edu

c
c     Numerical Analysis:
c     The Mathematics of Scientific Computing
c     D.R. Kincaid & E.W. Cheney
c     Brooks/Cole Publ., 1990
c
c     Section 4.7
c
c     Example of Conjugate Gradient method
c
c     file: cg.f
c
      parameter (n=4)
      dimension a(n,n),r(n),b(n),v(n),x(n),z(n)
      data (a(1,i),i=1,n) /420.,210.,140.,105./
      data (a(2,i),i=1,n) /210.,140.,105.,84./
      data (a(3,i),i=1,n) /140.,105.,84.,70./
      data (a(4,i),i=1,n) /105.,84.,70.,60./
      data (b(i),i=1,n) /875.,539.,399.,319./
      data m/15/
      data (x(i),i=1,n) /0.,0.,0.,0./
      data delta /1.0e-6/
c
      print *
      print *,' Example of Conjugate Gradient method'
      print *,' Section 4.7, Kincaid-Cheney'
      print *
c
      call residual(n,a,x,b,r)
      print *,' The residual vector corresponding to x(0) is:'
      print 7,r
      print *
c
      do 2 i=1,n
         v(i) = r(i)
 2    continue
c
      c = prod(n,r,r)
      do 6 k=1,M
         if (sqrt(prod(n,v,v)) .lt. delta) stop
         call mult(n,a,v,z)
         t = c/prod(n,v,z)
         do 3 i=1,n
            x(i) = x(i) + t*v(i)
 3       continue
         do 4 i=1,n
            r(i) = r(i) - t*z(i)
 4       continue
c
         d = prod(n,r,r)
c
         do 5 i=1,n
            v(i) = r(i) + (d/c)*v(i)
 5       continue
c
         c = d
c
         print *,' At iteration',k
         print *,' Solution is:'
         print 7,x
         print *,' Residual vector is:'
         print 7,r
         print *
 6    continue
c
 7    format (3x,4(e13.6,2x))
      stop
      end
c
      function prod(n,x,y)
c
c     compute the vector product
c
      dimension x(n),y(n)
      sum = 0.0
      do 2 i=1,n
         sum = sum + x(i)*y(i)
 2    continue
      prod = sum
c
      return
      end
c
      subroutine residual(n,a,x,b,r)
c
c     compute the residual r = b-Ax
c
      dimension a(n,n),x(n),b(n),r(n)
      call mult(n,a,x,r)
c
      do 2 i=1,n
         r(i) = b(i) - r(i)
 2    continue
c
      return
      end
c
      subroutine mult(n,a,x,y)
c
c     compute the matrix-vector product
c
      dimension a(n,n),x(n),y(n)
c
      do 3 i=1,n
         sum = 0.0
         do 2 j=1,n
            sum = sum + a(i,j)*x(j)
 2       continue
         y(i) = sum
 3    continue
c
      return
      end