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