# This is a shell archive. Remove anything before this line, # then unpack it by saving it in a file and typing "sh file". # # Wrapped by truffles!kincaid on Thu Mar 14 14:00:25 CST 1991 # Contents: elimit.f sqrt2.f nest.f epsi.f depsi.f ex2s22.f unstab1.f # unstab2.f instab.f ex1s31.f ex1s32.f ex2s32.f ex3s32.f ex1s33.f # ex3s34.f ex3s35.f ex6s35.f ex7s35.f laguerre.f forsub.f bacsub.f # pforsub.f pbacsub.f genlu.f doolt.f cholsky.f bgauss.f pbgauss.f # gauss.f paxeb.f yaec.f tri.f ex1s45.f ex2s45.f ex1s46.f ex2s46.f # jacobi.f ex3s46.f ex6s46.f steepd.f cg.f pcg.f ex1s51.f poweracc.f # ex2s51.f ipoweracc.f ex1s52.f qrshif.f ex1s53.f ex2s55.f ex3s55.f # coef.f fft.f adapta.f ex1s71.f ex2s71.f ex5s71.f ex6s71.f gauss5.f # romberg.f adapt.f taylor.f rk4.f rkfelberg.f taysys.f exs91.f exs92.f # exs93.f ex3s96.f mgrid1.f exs98.f mgrid2.f code-info.tex code-info.tty echo x - elimit.f sed 's/^@//' > "elimit.f" <<'@//E*O*F elimit.f//' 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 1.2 c c Example of slowly converging sequence for irrational number e c c c file: elimit.f c double precision x dimension x(1000) c print * print *,' Slowly converging sequence for irrational number e' print *,' Section 1.2, Kincaid-Cheney' print * c do 2 n=1,1000 x(n) = (1.0d0 + 1.0d0/real(n))**n 2 continue c print *,' 1, x(1) =',x(1) print *,' 10, x(10) =',x(10) print *,' 30, x(30) =',x(30) print *,' 50, x(50) =',x(50) print *,' 1000, x(1000) =',x(1000) print *,' exp(1.0) =',exp(1.0d0) print * print *,' error =',abs(x(1000) - exp(1.0d0)) c stop end @//E*O*F elimit.f// chmod u=rwx,g=x,o=x elimit.f echo x - sqrt2.f sed 's/^@//' > "sqrt2.f" <<'@//E*O*F sqrt2.f//' 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 1.2 c c Example of rapidly convergent sequence c c c file: sqrt2.f c double precision x c print * print *,' Rapidly converging sequence' print *,' Section 1.2, Kincaid-Cheney' print * c x = 2.0d0 print *,'k x' print *,1,x c do 2 k=2,4 x = x/2.0d0 + 1.0d0/x print *,k,x 2 continue c print *,' ',sqrt(2.0d0),'= sqrt(2)' c stop end @//E*O*F sqrt2.f// chmod u=rwx,g=x,o=x sqrt2.f echo x - nest.f sed 's/^@//' > "nest.f" <<'@//E*O*F nest.f//' 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 1.2 c c Example of Nested Multiplication c c c file: nest.f c dimension a(0:4) data (a(i),i=0,4) /-2.0,-5.0,7.0,-4.0,1.0/ data n,x /4,3.0/ c print * print *,' Nested Multiplication example' print *,' Section 1.2, Kincaid-Cheney' print * c p = a(n) do 2 k=n-1,0,-1 p = x*p + a(k) 2 continue c print *,' The value of the polynomial at x=3 is ',p c stop end @//E*O*F nest.f// chmod u=rwx,g=x,o=x nest.f echo x - epsi.f sed 's/^@//' > "epsi.f" <<'@//E*O*F epsi.f//' 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 2.1 c c Computing an approximate value of machine precision c c c file: epsi.f c print * print *,' Approximate value of machine precision' print *,' Section 2.1, Kincaid-Cheney' print * print *,' n computed 2**(-n)' c s = 1.0 c do 2 k=1,100 s = 0.5*s t = s + 1.0 if (t .le. 1.0) then s = 2.0*s t = 1.0/2.0**(k-1) print 3,k-1,s,t stop endif 2 continue c 3 format (i3,2x,2(e13.6,2x)) stop end @//E*O*F epsi.f// chmod u=rwx,g=x,o=x epsi.f echo x - depsi.f sed 's/^@//' > "depsi.f" <<'@//E*O*F depsi.f//' 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 2.1 c c Computing an approximate value of c double precision machine precision c c c file: depsi.f c double precision s,t c print * print *,' Approximate value of double precision machine epsilon' print *,' Section 2.1, Kincaid-Cheney' print * print *,' n computed 2**(-n)' c s = 1.0 c do 2 k=1,100 s = 0.5d0*s t = s + 1.0d0 if (t .le. 1.0d0) then s = 2.0d0*s t = 1.0d0/2.0d0**(k-1) print 3,k-1,s,t stop endif 2 continue c 3 format (i3,2(2x,d22.15)) stop end @//E*O*F depsi.f// chmod u=rwx,g=x,o=x depsi.f echo x - ex2s22.f sed 's/^@//' > "ex2s22.f" <<'@//E*O*F ex2s22.f//' 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 2.2 c c Example 2 c loss of significance in y = x - sin(x) c c c file: ex2s22.f c double precision t,s,yd c print * print *,' Using double precision to avoid loss of significance' print *,' Section 2.2, Kincaid-Cheney' print * c M = 7 x = 1.9 y = x - sin(x) yd = dble(x) - sin(dble(x)) p = (x**3/6.)*(1. - (x*x/20.)*(1. - (x*x/42.)*(1. - (x*x/72.)))) print *,' x =',x,' p =',p print *,' y =',y,' yd =',yd print * c t = x**3/6.0 s = t error = abs((s - yd)/yd) print *,' n s(n) rel. error' print 3,1,s,error c do 2 n=1,M t = -t*x*x/real((2.*real(n)+2.)*(2.*real(n)+3.)) s = s + t error = abs((s - yd)/yd) print 3,n+1,s,error 2 continue c 3 format(i3,4x,2(d22.15,4x)) stop end @//E*O*F ex2s22.f// chmod u=rwx,g=x,o=x ex2s22.f echo x - unstab1.f sed 's/^@//' > "unstab1.f" <<'@//E*O*F unstab1.f//' 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 2.3 c c Example of unstable sequence c c c file: unstab1.f c dimension x(0:15),y(0:15),e(0:15),f(0:15) double precision y c print * print *,' Unstable sequence' print *,' Section 2.3, Kincaid-Cheney' print * print 4,'n','x(n)','true','error','rel. error' c x(0) = 1.0 x(1) = 1.0/3.0 y(0) = 1.d0 y(1) = 1.d0/3.d0 c do 2 n=1,14 x(n+1) = (13./3.)*x(n) - (4./3.)*x(n-1) y(n+1) = (1.d0/3.d0)**(n+1) 2 continue c do 3 n=0,15 e(n) = real(dble(x(n)) - y(n)) f(n) = abs(e(n)/real(y(n))) print 5,n,x(n),real(y(n)),e(n),f(n) 3 continue c 4 format (a6,a10,a15,a16,a19) 5 format(1x,i5,4(2x,e13.6)) stop end @//E*O*F unstab1.f// chmod u=rwx,g=x,o=x unstab1.f echo x - unstab2.f sed 's/^@//' > "unstab2.f" <<'@//E*O*F unstab2.f//' 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 2.3 c c Example of unstable sequence c c c file: unstab2.f c dimension x(0:15),y(0:15),e(0:15),f(0:15) double precision y c print * print *,' Unstable sequence' print *,' Section 2.3, Kincaid-Cheney' print * print 4,'n','x(n)','true','error','rel. error' c x(0) = 1.0 x(1) = 4.0 y(0) = 1.d0 y(1) = 4.d0 c do 2 n=1,14 x(n+1) = (13./3.)*x(n) - (4./3.)*x(n-1) y(n+1) = (4.d0)**(n+1) 2 continue c do 3 n=0,15 e(n) = real(dble(x(n)) - y(n)) f(n) = abs(e(n)/real(y(n))) print 5,n,x(n),real(y(n)),e(n),f(n) 3 continue c 4 format (a6,a10,a15,a16,a19) 5 format(1x,i5,4(2x,e13.6)) stop end @//E*O*F unstab2.f// chmod u=rwx,g=x,o=x unstab2.f echo x - instab.f sed 's/^@//' > "instab.f" <<'@//E*O*F instab.f//' 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 2.3 c c Example of numerical instability c c c file: instab.f c print * print *,' Numerical instability' print *,' Section 2.3, Kincaid-Cheney' print * print *,' n y(n)' c e = exp(1.0) y = 1.0 print 3,1,y c do 2 n=2,20 y = e - y*real(n) print 3,n,y 2 continue c 3 format (i3,2x,e13.6) stop end @//E*O*F instab.f// chmod u=rwx,g=x,o=x instab.f echo x - ex1s31.f sed 's/^@//' > "ex1s31.f" <<'@//E*O*F ex1s31.f//' 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 3.1 c c Example of Bisection Method c c c file: ex1s31.f c external f data a,b,epsi,M,delta/-4.0,-3.0,1.0e-5,25,1.0e-5/ c print * print *,' Bisection Method to find the root of exp(x)=sin x' print *,' Section 3.1, Kincaid-Cheney' print * c call bisect(f,a,b,epsi,M) stop end c function f(x) f = exp(x) - sin(x) return end c subroutine bisect(f,a,b,epsi,M) c c bisecting the interval c external f fa = f(a) fb = f(b) error = b - a print *,' a =',a,' b =',b print *,' f(a) =',fa,' f(b) =',fb print * if (sign(1.0,fa) .eq. sign(1.0,fb)) then print 5 return endif print 3 c do 2 k=1,M error = error/2.0 c = a + error fc = f(c) print 4,k,c,fc,error if ((abs(fc) .lt. epsi) .or. (abs(error) .lt. delta)) return if (sign(1.0,fa) .ne. sign(1.0,fc)) then b = c fb = fc else a = c fa = fc endif 2 continue c 3 format (2x,' step',6x,'c',13x,'f(c)',10x,'error') 4 format (2x,i3,3(2x,e13.6)) 5 format (2x,'function does not change sign over interval',2e13.6) return end @//E*O*F ex1s31.f// chmod u=rwx,g=x,o=x ex1s31.f echo x - ex1s32.f sed 's/^@//' > "ex1s32.f" <<'@//E*O*F ex1s32.f//' 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 3.2 c c Example of Newton's Method c c c file: ex1s32.f c double precision x,x0,x1,f,fx,fp,epsi,delta data x0/-7.0d0/, M/8/, epsi/1.0e-12/, delta/1.0e-12/ f(x) = exp(x) - 1.5d0 - atan(x) fp(x) = exp(x) - 1.0d0/(1.0d0 + x*x) c print * print *,' Newton method example' print *,' Section 3.2, Kincaid-Cheney' print * c fx = f(x0) print 3 print 4,0,x0,fx if (abs(fx) .lt. epsi) stop do 2 k = 1,M x1 = x0 - fx/fp(x0) fx = f(x1) print 4,k,x1,fx if((abs(x1 - x0) .lt. delta) .or. (abs(fx) .lt. epsi)) stop x0 = x1 2 continue c 3 format(3x,'k',12x,'x',21x,'f(x)') 4 format(1x,i3,2x,d22.15,2x,d22.15) stop end @//E*O*F ex1s32.f// chmod u=rwx,g=x,o=x ex1s32.f echo x - ex2s32.f sed 's/^@//' > "ex2s32.f" <<'@//E*O*F ex2s32.f//' 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 3.2 c c Example of simple Newton's method c c c file: ex2s32.f c double precision x,f,fp data x/4.0d0/, M/4/ f(x) = x*x - 17.0d0 fp(x) = 2.0d0*x c print * print *,' Simple Newton method example' print *,' Section 3.2, Kincaid-Cheney' print * c fx = f(x) print 3 print 4,0,x,fx do 2 k = 1,M x = x - fx/fp(x) fx = f(x) print 4,k,x,fx 2 continue c 3 format (3x,' k',12x,'x',21x,'f(x)') 4 format(2x,i3,2x,d22.15,2x,d22.15) stop end @//E*O*F ex2s32.f// chmod u=rwx,g=x,o=x ex2s32.f echo x - ex3s32.f sed 's/^@//' > "ex3s32.f" <<'@//E*O*F ex3s32.f//' 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 3.2 c c Example of implicit function c c c file: ex3s32.f c double precision x,y,g,gp,gy,h data x/0.0d0/, y/1.0d0/, h/0.1/, M/100/, N/4/ g(x,y) = (3.0d0*x**4 - 1.0d0)*x**3 + (2.0d0*y**2 + 1.0d0)*y**3 - 3.0d0 gp(x,y) = (10.0d0*y**2 + 3.0d0)*y**2 c print * print *,' Implicit Function Example' print *,' Section 3.2, Kincaid-Cheney' print * c print 4 print 7,0,x,y,g(x,y) do 3 i=1,M x = x + h gy = g(x,y) print 5 do 2 j = 1,N y = y - gy/gp(x,y) gy = g(x,y) print 6,j,y,gy 2 continue print * print 7,i,x,y,gy 3 continue c 4 format(4x,'i',14x,'x',23x,'y',19x,'g(x,y)') 5 format(7x,'j',13x,'y',19x,'g(x,y)') 6 format(5x,i3,2x,d21.15,2x,d22.15) 7 format(2x,i3,2x,3(d22.15,2x)) stop end @//E*O*F ex3s32.f// chmod u=rwx,g=x,o=x ex3s32.f echo x - ex1s33.f sed 's/^@//' > "ex1s33.f" <<'@//E*O*F ex1s33.f//' 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 3.3 c c Example of Secant Method c c c file: ex1s33.f c data a/7.0/, b/8.0/, M/20/ data epsi/1.0e-6/, delta/1.0e-6/ fa = f(a) fb = f(b) c print * print *,' Secant Method Example' print *,' Section 3.3, Kincaid-Cheney' print * c print *,' n x(n) f(x(n))' print 3,0,a,fa print 3,1,b,fb c do 2 k = 2,M if (abs(fa) .le. abs(fb)) then tmp = a a = b b = tmp tmp = fa fa = fb fb = tmp endif s = (b - a)/(fb - fa) a = b fa = fb b = b - fb*s fb = f(b) print 3,k,b,fb if ((abs(fb) .lt. epsi) .or. (abs(a - b) .lt. delta)) stop 2 continue print 3,M,b,f(b) c 3 format(2x,i3,2(2x,e13.6)) stop end c function f(x) f = (((x + 4.0)*x + 6.0)*x + 9.0) - sinh(x) return end @//E*O*F ex1s33.f// chmod u=rwx,g=x,o=x ex1s33.f echo x - ex3s34.f sed 's/^@//' > "ex3s34.f" <<'@//E*O*F ex3s34.f//' 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 3.4 c c Example of Contractive Mapping Theorem c c c file: ex3s34.f c double precision x data x,M/4.0d0,20/ c print * print *,' Contractive Mapping Example' print *,' Section 3.4, Kincaid-Cheney' print * c print *,' k x' do 2 k=1,M x = 4.0d0 + (1.0d0/3.0d0)*sin(2.0d0*x) print 3,k,x 2 continue c 3 format (3x,i3,2x,d22.15) stop end @//E*O*F ex3s34.f// chmod u=rwx,g=x,o=x ex3s34.f echo x - ex3s35.f sed 's/^@//' > "ex3s35.f" <<'@//E*O*F ex3s35.f//' 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 3.5 c c Example of Horner's ALgorithm c c c file: ex3s35.f c dimension a(0:4),b(-1:3) data (a(i),i=0,4) /-2.0,-5.0,7.0,-4.0,1.0/ data n,z0 /4,3.0/ c print * print *,' Horner Method example' print *,' Section 3.5, Kincaid-Cheney' print * c b(n-1) = a(n) do 2 k=n-1,0,-1 b(k-1) = a(k) + z0*b(k) 2 continue c print *,' The value of the polynomial at x=3 is',b(-1) c stop end @//E*O*F ex3s35.f// chmod u=rwx,g=x,o=x ex3s35.f echo x - ex6s35.f sed 's/^@//' > "ex6s35.f" <<'@//E*O*F ex6s35.f//' 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 3.5 c c Example of Newton's Method on a polynomial of degree n c c Given starting point z c Coefficients in the polynomial: c a(0) is the constant term and a(n) is the coefficient of z**n c c c file: ex6s35.f c dimension a(0:4) data a/-2.0,-5.0,7.0,-4.0,1.0/ data n/4/ c print * print *,' Newton''s method on a given polynomial' print *,' Section 3.5, Kincaid-Cheney' print * c print *,' j alpha beta z' z = 0.0 call horner(n,a,z) stop end c subroutine horner(n,a,z) dimension a(0:n) data epsi/.5e-5/ do 3 j=1,8 alpha = a(n) beta = 0.0 do 2 k=n-1,0,-1 beta = alpha + z*beta alpha = a(k) + z*alpha 2 continue z1 = z - alpha/beta if (abs(z1-z) .lt. epsi) stop z = z1 print 4,j,alpha,beta,z 3 continue c 4 format(1x,i3,2x,3(e15.8,2x)) return end @//E*O*F ex6s35.f// chmod u=rwx,g=x,o=x ex6s35.f echo x - ex7s35.f sed 's/^@//' > "ex7s35.f" <<'@//E*O*F ex7s35.f//' 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 3.5 c c Example of Bairstow's method c applied to a polynomial of degree n c c c file: ex7s35.f c parameter (n=4,m=10) double precision a(0:n),b(0:n),c(0:n),u,v,xj data (a(j),j=0,n)/-2.0,-5.0,7.0,-4.0,1.0/ data u/3.0/ data v/-4.0/ c print * print *,' Bairstow''s Method example' print *,' Section 3.5, Kincaid-Cheney' print * print 10 c b(n) = a(n) c(n) = 0.0 c(n-1) = a(n) c do 3 j=1,M b(n-1) = a(n-1) + u*b(n) c do 2 k=n-2,0,-1 b(k) = a(k) + u*b(k+1) + v*b(k+2) c(k) = b(k+1) + u*c(k+1) + v*c(k+2) 2 continue c xj = c(0)*c(2) - c(1)**2 u = u+(c(1)*b(1) - c(2)*b(0))/xj v = v + (c(1)*b(0) - c(0)*b(1))/xj print 11,j,u,v,b(0),b(1) 3 continue c 10 format(1x,'n',12x,'u',23x,'v',17x,'b(0)',11x,'b(1)') 11 format(i2,1x,2(d22.15,2x),2(e13.6,2x)) stop end @//E*O*F ex7s35.f// chmod u=rwx,g=x,o=x ex7s35.f echo x - laguerre.f sed 's/^@//' > "laguerre.f" <<'@//E*O*F laguerre.f//' 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 3.5 c c Example of Laguerre's Method c c c file: laguerre.f c parameter (n=4, M=7) dimension a(0:n) implicit complex (a-h,o-z) data a/(-72.,68.),(1140.,636.),(-216.,190.),(-10.,-26.),(1.,0.)/ data z/(1.,1.)/ c print * print *,' Laguerre''s Method Example' print *,' Section 3.5, Kincaid-Cheney' print * print *,' Value of n and values of coefficients:' print *,n,(a(j),j=0,n) print * print *,' Starting point for the iteration:',z print * c sn = real(n) c do 3 k=1,m alpha = a(n) beta = (0.0,0.0) gamma = (0.0,0.0) c do 2 j=n-1,0,-1 gamma = z*gamma + beta beta = z*beta + alpha alpha = z*alpha + a(j) 2 continue c ca = -beta/alpha cb = ca*ca - 2.0*gamma/alpha c1 = sqrt((sn-1.0)*(sn*cb - ca*ca)) cc1 = (ca + c1)/sn cc2 = (ca - c1)/sn cc = cc2 if (abs(cc1) .gt. abs(cc2)) then cc = cc1 c2 = 1.0/cc z = z + c2 d = (ca - cc)/(sn - 1.0) rho1 = sqrt(sn)*c2 rho2 = sqrt(sn/(cc*cc + (sn-1)*d*d)) print *,'values of z, c2, rho1, and rho2:' print *,z,c2,rho1,rho2 print * 3 continue c stop end @//E*O*F laguerre.f// chmod u=rwx,g=x,o=x laguerre.f echo x - forsub.f sed 's/^@//' > "forsub.f" <<'@//E*O*F forsub.f//' 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.2 c c Example of Forward Substitution c c c file: forsub.f c parameter (n=3) dimension a(n,n),x(n),b(n) data (a(1,j),j=1,n) /6.0,0.0,0.0/ data (a(2,j),j=1,n) /3.0,2.0,0.0/ data (a(3,j),j=1,n) /4.0,2.0,1.0/ data (b(i),i=1,n) /6.0,5.0,7.0/ c print * print *,' Forward Substitution example' print *,' Section 4.2, Kincaid-Cheney' print * c do 3 i=1,n sum = b(i) do 2 j=1,i-1 sum = sum - a(i,j)*x(j) 2 continue x(i) = sum/a(i,i) 3 continue c do 4 i=1,n print 5,i,x(i) 4 continue c 5 format (1x,'x(',i2,') =',2x,e13.6) stop end @//E*O*F forsub.f// chmod u=rwx,g=x,o=x forsub.f echo x - bacsub.f sed 's/^@//' > "bacsub.f" <<'@//E*O*F bacsub.f//' 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.2 c c Example of Backward Substitution c c c file: bacsub.f c parameter (n=3) dimension a(n,n),x(n),b(n) data (a(1,j),j=1,n) /4.0,2.0,1.0/ data (a(2,j),j=1,n) /0.0,3.0,2.0/ data (a(3,j),j=1,n) /0.0,0.0,6.0/ data (b(i),i=1,n) /7.0,5.0,6.0/ c print * print *,' Backward Substitution example' print *,' Section 4.2, Kincaid-Cheney' print * c do 3 i=n,1,-1 sum = b(i) do 2 j=i+1,n sum = sum - a(i,j)*x(j) 2 continue x(i) = sum/a(i,i) 3 continue c do 4 i=1,n print 5,i,x(i) 4 continue c 5 format (1x,'x(',i2,') =',2x,e13.6) stop end @//E*O*F bacsub.f// chmod u=rwx,g=x,o=x bacsub.f echo x - pforsub.f sed 's/^@//' > "pforsub.f" <<'@//E*O*F pforsub.f//' 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.2 c c Example of Forward Substitution for a permuted system c c c file: pforsub.f c parameter (n=3) dimension a(n,n),x(n),b(n) integer p(n) data (a(1,j),j=1,n) /3.0,2.0,0.0/ data (a(2,j),j=1,n) /4.0,2.0,1.0/ data (a(3,j),j=1,n) /6.0,0.0,0.0/ data (b(i),i=1,n) /5.0,7.0,6.0/ data (p(i),i=1,n) /3,1,2/ c print * print *,' Forward Substitution for a permuted system' print *,' Section 4.2, Kincaid-Cheney' print * c do 3 i=1,n sum = b(p(i)) do 2 j=1,i-1 sum = sum - a(p(i),j)*x(j) 2 continue x(i) = sum/a(p(i),i) 3 continue c do 4 i=1,n print 5,i,x(i) 4 continue c 5 format (1x,'x(',i2,') =',2x,e13.6) stop end @//E*O*F pforsub.f// chmod u=rwx,g=x,o=x pforsub.f echo x - pbacsub.f sed 's/^@//' > "pbacsub.f" <<'@//E*O*F pbacsub.f//' 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.2 c c Example of Backward Substitution for a permuted system c c c file: pbacsub.f c parameter (n=3) dimension a(n,n),x(n),b(n) integer p(n) data (a(1,j),j=1,n) /0.0,3.0,2.0/ data (a(2,j),j=1,n) /0.0,0.0,6.0/ data (a(3,j),j=1,n) /4.0,2.0,1.0/ data (b(i),i=1,n) /5.0,6.0,7.0/ data (p(i),i=1,n) /3,1,2/ c print * print *,' Backward Substitution for a permuted system' print *,' Section 4.2, Kincaid-Cheney' print * c do 3 i=n,1,-1 sum = b(p(i)) do 2 j=i+1,n sum = sum - a(p(i),j)*x(j) 2 continue x(i) = sum/a(p(i),i) 3 continue c do 4 i=1,n print 5,i,x(i) 4 continue c 5 format (1x,'x(',i2,') =',2x,e13.6) stop end @//E*O*F pbacsub.f// chmod u=rwx,g=x,o=x pbacsub.f echo x - genlu.f sed 's/^@//' > "genlu.f" <<'@//E*O*F genlu.f//' 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.2 c c Example of general LU-factorization c c c file: genlu.f c parameter (n=3) dimension a(n,n),l(n,n),u(n,n) real l logical t(n) data (a(1,j),j=1,n) / 60.0,30.0,20.0/ data (a(2,j),j=1,n) / 30.0,20.0,15.0/ data (a(3,j),j=1,n) / 20.0,15.0,12.0/ data l(1,1),u(2,2),l(3,3) /6.0,2.0,4.0/ data (t(i),i=1,n) /.true.,.false.,.true./ c print * print *,' General LU-factorization example' print *,' Section 4.2, Kincaid-Cheney' print * c do 3 k=1,n sum = a(k,k) do 2 m=1,k-1 sum = sum - l(k,m)*u(m,k) 2 continue if (t(k)) then u(k,k) = sum/l(k,k) else l(k,k) = sum/u(k,k) endif 3 continue c do 9 k=1,n sum1 = a(k,k) do 4 m=1,k-1 sum1 = sum1 - l(k,m)*u(m,k) 4 continue u(k,k) = sum1/l(k,k) do 6 j=k+1,n sum2 = a(k,j) do 5 m=1,k-1 sum2 = sum2 - l(k,m)*u(m,j) 5 continue u(k,j) = sum2/l(k,k) 6 continue do 8 i=k+1,n sum3 = a(i,k) do 7 m=1,k-1 sum3 = sum3 - l(i,m)*u(m,k) 7 continue l(i,k) = sum3/u(k,k) 8 continue 9 continue c do 11 i=1,n do 10 j=1,n print 12,i,j,l(i,j),i,j,u(i,j) 10 continue 11 continue c 12 format (1x,'l(',i2,',',i2,') =',e13.6,7x, + 'u(',i2,',',i2,') =',e13.6) stop end @//E*O*F genlu.f// chmod u=rwx,g=x,o=x genlu.f echo x - doolt.f sed 's/^@//' > "doolt.f" <<'@//E*O*F doolt.f//' 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.2 c c Example of Doolittle's-factorization c c c file: doolt.f c parameter (n=3) dimension a(n,n),l(n,n),u(n,n) real l data (a(1,j),j=1,n) / 60.0,30.0,20.0/ data (a(2,j),j=1,n) / 30.0,20.0,15.0/ data (a(3,j),j=1,n) / 20.0,15.0,12.0/ c print * print *,' Doolittle-factorization example' print *,' Section 4.2, Kincaid-Cheney' print * c do 6 k=1,n l(k,k) = 1.0 do 3 j=k,n sum1 = a(k,j) do 2 m=1,k-1 sum1 = sum1 - l(k,m)*u(m,j) 2 continue u(k,j) = sum1 3 continue do 5 i=k+1,n sum2 = a(i,k) do 4 m=1,k-1 sum2 = sum2 - l(i,m)*u(m,k) 4 continue l(i,k) = sum2/u(k,k) 5 continue 6 continue c do 8 i=1,n do 7 j=1,n print 9,i,j,l(i,j),i,j,u(i,j) 7 continue 8 continue c 9 format (1x,'l(',i2,',',i2,') =',e13.6,7x, + 'u(',i2,',',i2,') =',e13.6) stop end @//E*O*F doolt.f// chmod u=rwx,g=x,o=x doolt.f echo x - cholsky.f sed 's/^@//' > "cholsky.f" <<'@//E*O*F cholsky.f//' 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.2 c c Example of Cholesky-factorization c c c file: cholsky.f c parameter (n=3) dimension a(n,n),l(n,n) real l data (a(1,j),j=1,n) / 60.0,30.0,20.0/ data (a(2,j),j=1,n) / 30.0,20.0,15.0/ data (a(3,j),j=1,n) / 20.0,15.0,12.0/ c print * print *,' Cholesky-factorization' print *,' Section 4.2, Kincaid-Cheney' print * c do 5 k=1,n sum1 = a(k,k) do 2 m=1,k-1 sum1 = sum1 - l(k,m)**2.0 2 continue l(k,k) = sqrt(sum1) do 4 i=k+1,n sum2 = a(i,k) do 3 m=1,k-1 sum2 = sum2 - l(i,m)*l(k,m) 3 continue l(i,k) = sum2/l(k,k) 4 continue 5 continue c do 7 i=1,n do 6 j=1,n print 8,i,j,l(i,j) 6 continue 7 continue c 8 format (1x,'l(',i2,',',i2,') =',e13.6) stop end @//E*O*F cholsky.f// chmod u=rwx,g=x,o=x cholsky.f echo x - bgauss.f sed 's/^@//' > "bgauss.f" <<'@//E*O*F bgauss.f//' 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.3 c c Example of basic Gaussian elimination c c c file: bgauss.f c parameter (n = 4) dimension a(n,n) data (a(1,j),j=1,n) /6.0,-2.0,2.0,4.0/ data (a(2,j),j=1,n) /12.0,-8.0,6.0,10.0/ data (a(3,j),j=1,n) /3.0,-13.0,9.0,3.0/ data (a(4,j),j=1,n) /-6.0,4.0,1.0,-18.0/ c print * print *,' Basic gaussian elimination' print *,' Section 4.3, Kincaid-Cheney' print * c do 4 k=1,n-1 do 3 i=k+1,n xmult = a(i,k)/a(k,k) a(i,k) = 0.0 do 2 j=k+1,n a(i,j) = a(i,j) - xmult*a(k,j) 2 continue 3 continue 4 continue c do 6 i=1,n do 5 j=1,n print 7,i,j,a(i,j) 5 continue 6 continue c 7 format (1x,'a(',i2,',',i2,') =',e13.6) stop end @//E*O*F bgauss.f// chmod u=rwx,g=x,o=x bgauss.f echo x - pbgauss.f sed 's/^@//' > "pbgauss.f" <<'@//E*O*F pbgauss.f//' 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.3 c c Example of basic Gaussian elimination with pivoting c c c file: pbgauss.f c parameter (n = 2) dimension a(n,n) integer p(n) data (a(1,j),j=1,n) /10.0e-8,1.0/ data (a(2,j),j=1,n) /1.0,1.0/ data (p(i),i=1,n) /2,1/ c print * print *,' Naive Gaussian elimination with pivoting' print *,' Section 4.3, Kincaid-Cheney' print * c do 4 k=1,n-1 do 3 i=k+1,n z = a(p(i),k)/a(p(k),k) a(p(i),k) = 0.0 do 2 j=k+1,n a(p(i),j) = a(p(i),j) - z*a(p(k),j) 2 continue 3 continue 4 continue c do 6 i=1,n do 5 j=1,n print 7,i,j,a(i,j) 5 continue 6 continue c 7 format (3x,'a(',i2,',',i2,') =',e13.6) stop end @//E*O*F pbgauss.f// chmod u=rwx,g=x,o=x pbgauss.f echo x - gauss.f sed 's/^@//' > "gauss.f" <<'@//E*O*F gauss.f//' 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.3 c c Example of Gaussian elimination with scaled row pivoting c c c file: gauss.f c parameter (n = 3) dimension a(n,n),s(n),b(n),x(n) integer p(n) real max data (a(1,j),j=1,n) /2.0,3.0,-6.0/ data (a(2,j),j=1,n) /1.0,-6.0,8.0/ data (a(3,j),j=1,n) /3.0,-2.0,1.0/ data (b(i),i=1,n) /-1.0,3.0,2.0/ c print * print *,' Gaussian elimination with scaled row pivoting' print *,' Section 4.3, Kincaid-Cheney' print * c do 3 i=1,n p(i) = i smax = 0.0 do 2 j=1,n smax = max(smax,abs(a(i,j))) 2 continue s(i) = smax 3 continue c do 7 k=1,n-1 rmax = 0.0 do 4 i=k,n r = abs(a(p(i),k))/s(p(i)) if (r .gt. rmax) then j = i rmax = r endif 4 continue c pk = p(j) p(j) = p(k) p(k) = pk c do 6 i=k+1,n z = a(p(i),k)/a(p(k),k) a(p(i),k) = z do 5 j=k+1,n a(p(i),j) = a(p(i),j) - z*a(p(k),j) 5 continue 6 continue 7 continue c do 9 k=1,n-1 do 8 i=k+1,n b(p(i)) = b(p(i)) - a(p(i),k)*b(p(k)) 8 continue 9 continue do 11 i=n,1,-1 sum = b(p(i)) do 10 j=i+1,n sum = sum - a(p(i),j)*x(j) 10 continue x(i) = sum/a(p(i),i) 11 continue c do 12 i=1,n print 13,i,x(i) 12 continue c 13 format (3x,'x(',i2,') =',e13.6) stop end @//E*O*F gauss.f// chmod u=rwx,g=x,o=x gauss.f echo x - paxeb.f sed 's/^@//' > "paxeb.f" <<'@//E*O*F paxeb.f//' 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.3 c c Example of Gaussian elimination with scaled row pivoting c Solving Lz=Pb and Ux=z c c c file: paxeb.f c parameter (n = 3) dimension a(n,n),s(n),b(n),x(n),z(n) integer p(n) data (a(1,j),j=1,n) /2.0,3.0,-6.0/ data (a(2,j),j=1,n) /1.0,-6.0,8.0/ data (a(3,j),j=1,n) /3.0,-2.0,1.0/ data (b(i),i=1,n) /-1.0,3.0,2.0/ c print * print *,' Gaussian elimination with scaled row pivoting' print *,' Section 4.3, Kincaid-Cheney' print * c call gauss(n,a,s,p) c do 3 i=1,n sum = b(p(i)) do 2 j=1,i-1 sum = sum - a(p(i),j)*z(j) 2 continue z(i) = sum 3 continue c do 5 i=n,1,-1 sum = z(i) do 4 j=i+1,n sum = sum - a(p(i),j)*x(j) 4 continue x(i) = sum/a(p(i),i) 5 continue c do 6 i=1,n print 7,i,x(i) 6 continue c 7 format (3x,'x(',i2,') =',e13.6) stop end c subroutine gauss(n,a,s,p) c c gaussian elimination c dimension a(n,n),s(n) integer p(n) real max c do 9 i=1,n p(i) = i smax = 0.0 do 8 j=1,n smax = max(smax,abs(a(i,j))) 8 continue s(i) = smax 9 continue c do 13 k=1,n-1 rmax = 0.0 do 10 i=k,n r = abs(a(p(i),k))/s(p(i)) if (r .gt. rmax) then j = i rmax = r endif 10 continue c pk = p(j) p(j) = p(k) p(k) = pk c do 12 i=k+1,n z = a(p(i),k)/a(p(k),k) a(p(i),k) = z do 11 j=k+1,n a(p(i),j) = a(p(i),j) - z*a(p(k),j) 11 continue 12 continue 13 continue c return end @//E*O*F paxeb.f// chmod u=rwx,g=x,o=x paxeb.f echo x - yaec.f sed 's/^@//' > "yaec.f" <<'@//E*O*F yaec.f//' 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.3 c c Example of gaussian elimination with scaled row pivoting c Solving U(trans)z=c and L(trans)Py=z c c c file: yaec.f c parameter (n = 3) dimension a(n,n),s(n),c(n),y(n),z(n) integer p(n) data (a(1,j),j=1,n) /2.0,3.0,-6.0/ data (a(2,j),j=1,n) /1.0,-6.0,8.0/ data (a(3,j),j=1,n) /3.0,-2.0,1.0/ data (c(i),i=1,n) /6.0,-5.0,3.0/ c print * print *,' Gaussian elimination with scaled row pivoting' print *,' Section 4.3, Kincaid-Cheney' print * c call gauss(n,a,s,p) c do 3 j=1,n sum = c(j) do 2 i=1,j-1 sum = sum - a(p(i),j)*z(i) 2 continue z(j) = sum/a(p(j),j) 3 continue do 5 j=n,1,-1 sum = z(j) do 4 i=j+1,n sum = sum - a(p(i),j)*y(p(i)) 4 continue y(p(j)) = sum 5 continue c do 6 i=1,n print 7,i,y(i) 6 continue c 7 format (3x,'y(',i2,') =',e13.6) stop end c subroutine gauss(n,a,s,p) c c gaussian elimination c dimension a(n,n),s(n) integer p(n) real max c do 9 i=1,n p(i) = i smax = 0.0 do 8 j=1,n smax = max(smax,abs(a(i,j))) 8 continue s(i) = smax 9 continue c do 13 k=1,n-1 rmax = 0.0 do 10 i=k,n r = abs(a(p(i),k))/s(p(i)) if (r .gt. rmax) then j = i rmax = r endif 10 continue c pk = p(j) p(j) = p(k) p(k) = pk c do 12 i=k+1,n z = a(p(i),k)/a(p(k),k) a(p(i),k) = z do 11 j=k+1,n a(p(i),j) = a(p(i),j) - z*a(p(k),j) 11 continue 12 continue 13 continue c return end @//E*O*F yaec.f// chmod u=rwx,g=x,o=x yaec.f echo x - tri.f sed 's/^@//' > "tri.f" <<'@//E*O*F tri.f//' 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.3 c c Example of tridiagonal system c c c file: tri.f c parameter (n=10) real a(n),d(n),c(n),b(n),x(n) c do 2 i=1,n d(i) = 2.0 a(i) = 0.5 c(i) = 0.5 b(i) = 3.0 2 continue b(1) = 2.5 b(n) = 2.5 c print * print *,' Tridiagonal system example' print *,' Section 4.3, Kincaid-Cheney' print * c call tri(n,a,d,c,b,x) c do 3 i=1,n print 4,i,x(i) 3 continue c 4 format (3x,'x(',i2,') =',e13.6) stop end c subroutine tri(n,a,d,c,b,x) c c tridiagonalize the matrix c real a(n),d(n),c(n),b(n),x(n) real xmult integer i c do 5 i = 2,n xmult = a(i-1)/d(i-1) d(i) = d(i) - xmult*c(i-1) b(i) = b(i) - xmult*b(i-1) 5 continue x(n) = b(n)/d(n) do 6 i = n-1,1,-1 x(i) = (b(i) - c(i)*x(i+1))/d(i) 6 continue c return end @//E*O*F tri.f// chmod u=rwx,g=x,o=x tri.f echo x - ex1s45.f sed 's/^@//' > "ex1s45.f" <<'@//E*O*F ex1s45.f//' 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.5 c c Example of Neumann series to compute the inverse of a matrix c c c file: ex1s45.f c parameter (n=3) dimension a(n,n),c(n,n),s(n,n),t(n,n) data (a(1,i),i=1,n)/0.1,0.2,0.3/ data (a(2,i),i=1,n)/-0.1,0.0,0.1/ data (a(3,i),i=1,n)/-0.3,-0.2,-0.1/ c print * print *,' Neumann series example' print *,' Section 4.5, Kincaid-Cheney' print * c call setI(n,s) call setI(n,t) do 2 k=1,20 call mult(n,a,t,c) call store(n,c,t) call add(n,s,t) print *,' k=',k call prt(n,s) print * 2 continue stop end c subroutine setI(n,s) c c initializing the identity matrix c dimension s(n,n) c do 3 i=1,n do 2 j=1,n s(i,j) = 0.0 2 continue s(i,i) = 1.0 3 continue c return end c subroutine mult(n,a,b,c) c c computing the kth power of the given matrix c dimension a(n,n),c(n,n),b(n,n) c do 3 i=1,n do 3 j=1,n sum = 0.0 do 2 k=1,n sum = sum + a(i,k)*b(k,j) 2 continue c(i,j) = sum 3 continue c return end c subroutine store(n,a,b) c c storing the kth power of the matrix c dimension a(n,n),b(n,n) c do 2 i=1,n do 2 j=1,n b(i,j) = a(i,j) 2 continue c return end c subroutine add(n,s,t) c c computing the partial sums c dimension s(n,n),t(n,n) c do 2 i=1,n do 2 j=1,n s(i,j) = s(i,j) + t(i,j) 2 continue c return end c subroutine prt(n,A) c c printing the inverse of the matrix at the kth iteration c dimension A(n,n) c do 2 i=1,n print 3,(A(i,j),j=1,n) 2 continue c 3 format (3x,3(e13.6,2x)) c return end @//E*O*F ex1s45.f// chmod u=rwx,g=x,o=x ex1s45.f echo x - ex2s45.f sed 's/^@//' > "ex2s45.f" <<'@//E*O*F ex2s45.f//' 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.5 c c Example of Gaussian elimination followed by iterative improvement c c c file: ex2s45.f c parameter (n=4,ia=4) dimension a(ia,n),b(n),r(n),l(n),s(n),x(n) dimension aa(ia,n),e(n) data (a(1,j),j=1,n) /420.,210.,140.,105./ data (a(2,j),j=1,n) /210.,140.,105.,84./ data (a(3,j),j=1,n) /140.,105.,84.,70./ data (a(4,j),j=1,n) /105.,84.,70.,60./ data (b(i),i=1,n) /875.,539.,399.,319./ c print * print *,' Gaussian elimination followed by iterative improvement' print *,' Section 3.5, Kincaid-Cheney' print * c c save a c do 2 i=1,n do 2 j=1,n aa(i,j) = a(i,j) 2 continue c call gauss(a,ia,n,l,s) call solve(a,ia,n,l,b,x) print *,' Initial solution is:' print 5,x print * c do 4 k=1,5 call residual(aa,ia,n,b,x,r) print *,' At iteration',k print *,' Residual vector is:' print 5,r call solve(a,ia,n,l,r,e) print *,' Error vector is:' print 5,e c do 3 i=1,n x(i) = x(i) + e(i) 3 continue c print *,' Solution is:' print 5,x print * 4 continue c 5 format (3x,4(e13.6,2x)) stop end c subroutine residual(a,ia,n,b,x,r) c c computes residual vector r = b-Ax c dimension a(ia,n),b(n),r(n),x(n) double precision sum do 3 i=1,n sum = 0.0 do 2 j=1,n sum = sum + a(i,j)*x(j) 2 continue r(i) = b(i) - sum 3 continue c return end c subroutine gauss(a,ia,n,l,s) c c gaussian elimination with scaled row pivoting c dimension a(ia,n),l(n),s(n) real max do 2 i=1,n l(i) = i s(i) = 0.0 do 2 j=1,n s(i) = max(s(i),abs(a(i,j))) 2 continue do 4 k=1,n-1 rmax = 0.0 do 3 i=k,n r = abs(a(l(i),k))/s(l(i)) if (r .gt. rmax) then j = i rmax = r endif 3 continue c lk = l(j) l(j) = l(k) l(k) = lk do 4 i = k+1,n xmult = a(l(i),k)/a(lk,k) a(l(i),k) = xmult do 4 j = k+1,n a(l(i),j) = a(l(i),j) - xmult*a(lk,j) 4 continue c return end c subroutine solve(a,ia,n,l,b,x) c c back substitution LUx = b c dimension a(ia,n),l(n),b(n),x(n),z(4) c do 3 i=1,n sum = b(l(i)) do 2 j=1,i-1 sum = sum - a(l(i),j)*z(j) 2 continue z(i) = sum 3 continue do 5 i=n,1,-1 sum = z(i) do 4 j=i+1,n sum = sum - a(l(i),j)*x(j) 4 continue x(i) = sum/a(l(i),i) 5 continue c return end @//E*O*F ex2s45.f// chmod u=rwx,g=x,o=x ex2s45.f echo x - ex1s46.f sed 's/^@//' > "ex1s46.f" <<'@//E*O*F ex1s46.f//' 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.6 c c Example of Jacobi method and Gauss-Seidel method c c c file: ex1s46.f c data M/50/ c print * print *,' Jacobi method example' print *,' Section 4.6, Kincaid-Cheney' print * c x1 = 0. x2 = 0. print *,' k x1 x2' print 4,0,x1,x2 do 2 k=1,M y1 = (6./7.)*x2 + (3./7.) y2 = (8./9.)*x1 - (4./9.) x1 = y1 x2 = y2 if (0 .eq. mod(k,10)) print 4,k,x1,x2 2 continue c c Gauss-Seidel method c print * print * print *,' Gauss-Seidel method example' print *,' Section 4.6, Kincaid-Cheney' print * c x1 = 0. x2 = 0. print *,' k x1 x2' print 4,0,x1,x2 do 3 k=1,M x1 = (6./7.)*x2 + (3./7.) x2 = (8./9.)*x1 - (4./9.) if (0 .eq. mod(k,10)) print 4,k,x1,x2 3 continue c 4 format (3x,i2,2x,2(e13.6,2x)) stop end @//E*O*F ex1s46.f// chmod u=rwx,g=x,o=x ex1s46.f echo x - ex2s46.f sed 's/^@//' > "ex2s46.f" <<'@//E*O*F ex2s46.f//' 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.6 c c Example of Richardson's method c c c file: ex2s46.f c parameter (n=3) dimension r(n),x(n) data (x(i),i=1,n) /0.0,0.0,0.0/ c print * print *,' Richardson method example' print *,' Section 4.6, Kincaid-Cheney' print * c do 3 k=1,100 r(1) = (11./18.) - x(1) - .5*x(2) - (1./3.)*x(3) r(2) = (11./18.) - x(2) - .5*x(3) - (1./3.)*x(1) r(3) = (11./18.) - x(3) - .5*x(1) - (1./3.)*x(2) do 2 j=1,n x(j) = x(j) + r(j) 2 continue if (0 .eq. mod(k,10)) then print *,' At iteration',k print 4,r print 5,x print * endif 3 continue c 4 format (3x,'r is:',3(e13.6,2x)) 5 format (3x,'x is:',3(e13.6,2x)) stop end @//E*O*F ex2s46.f// chmod u=rwx,g=x,o=x ex2s46.f echo x - jacobi.f sed 's/^@//' > "jacobi.f" <<'@//E*O*F jacobi.f//' 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.6 c c Example of Jacobi's method c c c file: jacobi.f c parameter (n=3) dimension a(n,n),b(n),u(n),x(n) data (a(1,j),j=1,n) /1.0,0.5,0.3333/ data (a(2,j),j=1,n) /0.3333,1.0,0.5/ data (a(3,j),j=1,n) /0.5,0.3333,1.0/ data (b(i),i=1,n) /.6111,.6111,.6111/ data (x(i),i=1,n) /0.0,0.0,0.0/ c print * print *,' Jacobi method example' print *,' Section 4.6, Kincaid-Cheney' print * c do 7 k=1,100 do 3 i=1,n d = 1.0/a(i,i) b(i) = d*b(i) do 2 j=1,n a(i,j) = d*a(i,j) 2 continue 3 continue do 5 i=1,n sum = 0.0 do 4 j=1,n if (j .eq. i) goto 4 sum = sum + a(i,j)*x(j) 4 continue u(i) = b(i) - sum 5 continue do 6 i=1,n x(i) = u(i) 6 continue if (0 .eq. mod(k,10)) then print *,' At iteration',k print 8,x print * endif 7 continue c 8 format (3x,'x is:',3(e13.6,2x)) stop end @//E*O*F jacobi.f// chmod u=rw,g=,o= jacobi.f echo x - ex3s46.f sed 's/^@//' > "ex3s46.f" <<'@//E*O*F ex3s46.f//' 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.6 c c Example of Gauss-Seidel method c c c file: ex3s46.f c dimension x(3) data x/0.,0.,0./ c print * print *,' Gauss-Seidel example' print *,' Section 4.6, Kincaid-Cheney' print * c print *,' k x1 x2 x3' print 3,0,x do 2 k=1,20 x(1) = 0.5*x(2) + 1.0 x(2) = -(1.0/6.0)*x(1) + (1.0/3.0)*x(3) - (2.0/3.0) x(3) = -(0.5)*x(1) + (3.0/8.0)*x(2) + (5.0/8.0) print 3,k,x 2 continue c 3 format (3x,i2,2x,3(e13.6,2x)) stop end @//E*O*F ex3s46.f// chmod u=rwx,g=x,o=x ex3s46.f echo x - ex6s46.f sed 's/^@//' > "ex6s46.f" <<'@//E*O*F ex6s46.f//' 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.6 c c Example of Chebyshev acceleration on Jacobi method c c c file: ex6s46.f c parameter (n=4,m=50) dimension u(n),v(n),g(n,n),c(n) data (g(1,j),j=1,n) /0.0,0.25,0.25,0.0/ data (g(2,j),j=1,n) /0.25,0.0,0.0,0.25/ data (g(3,j),j=1,n) /0.25,0.0,0.0,0.25/ data (g(4,j),j=1,n) /0.0,0.25,0.25,0.0/ data c /-1.0,0.0,1.0,-1.0/ data u /0.0,0.0,0.0,0.0/ data delta /1.0e-4/ c print * print *,' Chebyshev acceleration example' print *,' Section 4.6, Kincaid-Cheney' print * c b = 0.5 a = -b gamma = 2.0/(2.0 - (b + a)) alpha = (0.5*(b - a)/(2.0 - (b + a)))**2 print 3,0,u call extrap(gamma,n,G,c,u,v) print 3,1,v rho = 1.0/(1.0 - 2.0*alpha) call cheb(rho,gamma,n,g,c,u,v) print 3,2,u do 2 k=3,m,2 rho = 1.0/(1.0 - rho*alpha) call cheb(rho,gamma,n,G,c,v,u) print 3,k,v rho = 1.0/(1.0 - rho*alpha) call cheb(rho,gamma,n,G,c,u,v) print 3,k+1,u if (vnorm(n,u,v) .lt. delta) stop 2 continue c 3 format (3x,i2,2x,4(e13.6,2x)) stop end c subroutine extrap(gamma,n,G,c,u,v) c c extrapolation procedure c dimension u(n),v(n),G(n,n),c(n) c do 2 i=1,n v(i) = gamma*c(i) + (1.0 - gamma)*u(i) do 2 j=1,n v(i) = v(i) + gamma*g(i,j)*u(j) 2 continue c return end c subroutine cheb(rho,gamma,n,G,c,u,v) c c chebyshev acceleration c dimension u(n),v(n),G(n,n),c(n) c do 2 i=1,n u(i) = gamma*rho*c(i)+rho*(1.0-gamma)*v(i)+(1.0-rho)*u(i) 2 continue c do 3 i=1,n do 3 j=1,n u(i) = u(i) + rho*gamma*G(i,j)*v(j) 3 continue c return end c function vnorm(n,x,y) real x(n),y(n) c c compute max-norm of x - y c temp = 0.0 do 2 i=1,n temp = max(temp,abs(x(i) - y(i))) 2 continue vnorm = temp c return end @//E*O*F ex6s46.f// chmod u=rwx,g=x,o=x ex6s46.f echo x - steepd.f sed 's/^@//' > "steepd.f" <<'@//E*O*F steepd.f//' 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 Steepest Descent method c c c file: steepd.f c parameter (n=4) dimension a(n,n),b(n),x(n),v(n),y(n) data (a(1,j),j=1,n) /420.,210.,140.,105./ data (a(2,j),j=1,n) /210.,140.,105.,84./ data (a(3,j),j=1,n) /140.,105.,84.,70./ data (a(4,j),j=1,n) /105.,84.,70.,60./ data (b(i),i=1,n) /875.,539.,399.,319./ data m/200/ data (x(i),i=1,n) /0.,0.,0.,0./ c print * print *,' Steepest Descent method example' print *,' Section 4.7, Kincaid-Cheney' print * print *,' At iteration ',0,' initial x is:' print 5,x print * c do 4 k=1,m call mult(n,a,x,y) do 2 i=1,n v(i) = b(i) - y(i) 2 continue c = prod(n,v,v) call mult(n,a,v,y) d = prod(n,v,y) t = c/d do 3 i=1,n x(i) = x(i) + t*v(i) 3 continue if (0 .eq. mod(k,10)) then print *,' At iteration ',k,' solution is:' print 5,x print * endif 4 continue c 5 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 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) 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 return end @//E*O*F steepd.f// chmod u=rw,g=,o= steepd.f echo x - cg.f sed 's/^@//' > "cg.f" <<'@//E*O*F cg.f//' 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 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 *,' Iterative refinement example' print *,' Section 4.5, 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 @//E*O*F cg.f// chmod u=rwx,g=x,o=x cg.f echo x - pcg.f sed 's/^@//' > "pcg.f" <<'@//E*O*F pcg.f//' 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 preconditioned conjugate gradient method c (conditioned with Jacobi method) c c c file: pcg.f c parameter (n=4) dimension a(n,n),r(n),b(n),v(n),z(n),x(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 *,' Conjugate Gradient preconditioned with Jacobi matrix' print *,' Section 4.5, 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 z(i) = r(i)/a(i,i) v(i) = z(i) 2 continue c c = prod(n,z,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) z(i) = r(i)/a(i,i) 4 continue c d = prod(n,z,r) c do 5 i=1,n v(i) = z(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) c 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 vector 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 @//E*O*F pcg.f// chmod u=rwx,g=x,o=x pcg.f echo x - ex1s51.f sed 's/^@//' > "ex1s51.f" <<'@//E*O*F ex1s51.f//' 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 5.1 c c Example of power method c c c file: ex1s51.f c parameter (n=3,m=30) dimension a(n,n),x(n),y(n) data (a(1,j),j=1,n) /6.,5.,-5./ data (a(2,j),j=1,n) /2.,6.,-2./ data (a(3,j),j=1,n) /2.,5.,-1./ data (x(j),j=1,n) /-1.,1.,1./ c print * print *,' Power method example' print *,' Section 5.1, Kincaid-Cheney' print * print 3 print 4,0,x c do 2 k=1,m call prod(n,a,x,y) r = y(2)/x(2) call dot(n,x,y,t) call dot(n,x,x,s) s = t/s call normal(n,y) call store(n,y,x) print 4,k,x,r,s 2 continue c 3 format (1x,'k',7x,'x(1)',10x,'x(2)',10x,'x(3)',9x,'ratio', + 5x,'rayleigh quot.') 4 format (1x,i2,1x,5(e13.6,1x)) stop end c subroutine dot(n,x,y,t) c c compute dot product t= c dimension x(n),y(n) c t = 0.0 do 2 j=1,n t = t + x(j)*y(j) 2 continue c return end c subroutine prod(n,a,x,y) c c compute y=ax c dimension a(n,n),x(n),y(n) c do 2 i=1,n y(i) = 0.0 do 2 j=1,n y(i) = y(i) + a(i,j)*x(j) 2 continue c return end c subroutine store (n,x,y) c c put x into y c dimension x(n),y(n) c do 2 j=1,n y(j) = x(j) 2 continue c return end c subroutine norm (n,x,t) c c compute t = max-norm of x c real x(n),max c t = 0.0 do 2 j=1,n t = max(t,abs(x(j))) 2 continue c return end c subroutine normal(n,x) c c normalize x with max norm c dimension x(n) c call norm (n,x,t) do 2 j=1,n x(j) = x(j)/t 2 continue c return end @//E*O*F ex1s51.f// chmod u=rw,g=,o= ex1s51.f echo x - poweracc.f sed 's/^@//' > "poweracc.f" <<'@//E*O*F poweracc.f//' 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 5.1 c c Example to test power method for eigenvalues c c c file: poweracc.f c parameter (n=3,m=30) dimension a(n,n),x(n),y(n),at(m),r(m) data (a(1,j),j=1,n)/6.,5.,-5./ data (a(2,j),j=1,n)/2.,6.,-2./ data (a(3,j),j=1,n)/2.,5.,-1./ data (x(j),j=1,n)/-1.,1.,1./ c print * print *,' Power method with Aitken Acceleration' print *,' Section 5.1, Kincaid-Cheney' print * print 3 print 4,0,x c do 2 k=1,m call prod(n,a,x,y) r(k) = y(2)/x(2) call dot(n,x,y,t) call dot(n,x,x,s) if ((k .ge. 3) .and. (k .le. m-2)) then at(k) = (r(k-2)*r(k) - r(k-1)**2)/(r(k)-2.*r(k-1)+r(k-2)) endif s = t/s call normal(n,y) call store(n,y,x) print 4,k,x,r(k),s,at(k) 2 continue c 3 format (1x,'k',7x,'x(1)',10x,'x(2)',10x,'x(3)',9x,'ratio', + 5x,'rayleigh quot.',3x,'acc. quot.') 4 format (1x,i2,1x,6(e13.6,1x)) stop end c subroutine dot(n,x,y,t) c c compute dot product t= c dimension x(n),y(n) c t = 0.0 do 2 j=1,n t = t + x(j)*y(j) 2 continue c return end c subroutine prod(n,a,x,y) c c compute y=ax c dimension a(n,n),x(n),y(n) c do 2 i=1,n y(i) = 0.0 do 2 j=1,n y(i)=y(i)+a(i,j)*x(j) 2 continue c return end c subroutine store(n,x,y) c c put x into y c dimension x(n),y(n) c do 2 j=1,n y(j)=x(j) 2 continue c return end c subroutine norm(n,x,t) c c compute t= max-norm of x c real x(n),max c t = 0.0 do 2 j=1,n t = max(t,abs(x(j))) 2 continue c return end c subroutine normal(n,x) c c normalize x with max norm c dimension x(n) c call norm (n,x,t) do 2 j=1,n x(j) = x(j)/t 2 continue c return end @//E*O*F poweracc.f// chmod u=rw,g=,o= poweracc.f echo x - ex2s51.f sed 's/^@//' > "ex2s51.f" <<'@//E*O*F ex2s51.f//' 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 5.1 c c Example of Inverse Power method c c c file: ex2s51.f c parameter (n=3,m=25) dimension a(n,n),x(n),xx(n),y(n),s(n) integer p(n) data (a(1,j),j=1,n) /6.,5.,-5./ data (a(2,j),j=1,n) /2.,6.,-2./ data (a(3,j),j=1,n) /2.,5.,-1./ data (x(j),j=1,n) /3.,7.,-13./ c print * print *,' Inverse Power method example' print *,' Section 5.1, Kincaid-Cheney' print * print 3 print 4,0,x c call gauss (n,n,a,p,s) c do 2 k=1,m call store (n,x,xx) call solve (n,n,a,p,xx,y) r = y(1)/x(1) call dot(n,x,y,t) call dot(n,x,x,ss) ss = t/ss call normal (n,y) call store (n,y,x) print 4,k,x,r,ss 2 continue c 3 format (1x,'k',7x,'x(1)',10x,'x(2)',10x,'x(3)',9x,'ratio', + 5x,'rayleigh quot.') 4 format (1x,i2,1x,5(e13.6,1x)) stop end c subroutine gauss(n,ia,a,p,s) c c gaussian elimination c dimension a(ia,n),s(n) integer p(n) real max c do 3 i=1,n p(i) = i smax = 0.0 do 2 j=1,n smax = max(smax,abs(a(i,j))) 2 continue s(i) = smax 3 continue c do 7 k=1,n-1 rmax = 0.0 do 4 i=k,n r = abs(a(p(i),k))/s(p(i)) if (r .gt. rmax) then j = i rmax = r endif 4 continue pk = p(j) p(j) = p(k) p(k) = pk do 6 i=k+1,n z = a(p(i),k)/a(p(k),k) a(p(i),k) = z do 5 j=k+1,n a(p(i),j) = a(p(i),j) - z*a(p(k),j) 5 continue 6 continue 7 continue c return end c subroutine dot(n,x,y,t) c c compute dot product t= c dimension x(n),y(n) c t = 0.0 do 2 j=1,n t = t + x(j)*y(j) 2 continue c return end c subroutine prod(n,a,x,y) c c compute y = ax c dimension a(n,n),x(n),y(n) c do 2 i=1,n y(i) = 0.0 do 2 j=1,n y(i) = y(i)+a(i,j)*x(j) 2 continue c return end c subroutine store (n,x,y) c c put x into y c dimension x(n),y(n) c do 2 j=1,n y(j) = x(j) 2 continue c return end c subroutine norm (n,x,t) c c compute t = max-norm of x c real x(n),max c t = 0.0 do 2 j=1,n t = max(t,abs(x(j))) 2 continue c return end c subroutine normal(n,x) c c normalize x with max norm c dimension x(n) c call norm (n,x,t) do 2 j=1,n x(j) = x(j)/t 2 continue c return end c subroutine solve(ia,n,a,p,b,x) c c solve for x c dimension a(ia,n),b(n),x(n) integer p(n) c do 2 k=1,n-1 do 2 i=k+1,n b(p(i)) = b(p(i)) - a(p(i),k)*b(p(k)) 2 continue c do 4 i=n,1,-1 sum = b(p(i)) do 3 j=i+1,n sum = sum - a(p(i),j)*x(j) 3 continue x(i) = sum/a(p(i),i) 4 continue c return end @//E*O*F ex2s51.f// chmod u=rw,g=,o= ex2s51.f echo x - ipoweracc.f sed 's/^@//' > "ipoweracc.f" <<'@//E*O*F ipoweracc.f//' c c Numerical Analysis: c The Mathematics of Scientific Computing c D.R. Kincaid & E.W. Cheney c Brooks/Cole., 1990 c c Section 5.1 c c Example to test inverse power method for eigenvalues c c c file: ipoweracc.f c parameter (n=3,m=12) dimension a(n,n),x(n),xx(n),y(n),s(n),at(m),r(m) integer p(n) data (a(1,j),j=1,n) /6.,5.,-5./ data (a(2,j),j=1,n) /2.,6.,-2./ data (a(3,j),j=1,n) /2.,5.,-1./ data (x(j),j=1,n) /3.,7.,-13./ c print * print *,' Inverse power method with Aitken Acceleration' print *,' Section 5.1, Kincaid-Cheney' print * print 3 print 4,0,x c call gauss(n,n,a,p,s) c do 2 k=1,m call store (n,x,xx) call solve (n,n,a,p,xx,y) r(k) = y(1)/x(1) call dot(n,x,y,t) call dot(n,x,x,ss) if ((k .ge. 3) .and. (k .le. m-2)) then at(k) = (r(k-2)*r(k) - r(k-1)**2)/(r(k)-2.*r(k-1)+r(k-2)) endif ss = t/ss call normal(n,y) call store(n,y,x) print 4,k,x,r(k),ss,at(k) 2 continue c 3 format (1x,'k',7x,'x(1)',10x,'x(2)',10x,'x(3)',9x,'ratio', + 5x,'rayleigh quot.',3x,'acc. quot.') 4 format (1x,i2,1x,6(e13.6,1x)) stop end c subroutine gauss(n,ia,a,p,s) c c gaussian elimination c dimension a(ia,n),s(n) integer p(n) c do 3 i=1,n p(i) = i smax = 0.0 do 2 j=1,n smax = amax1(smax,abs(a(i,j))) 2 continue s(i) = smax 3 continue c do 7 k=1,n-1 rmax = 0.0 do 4 i=k,n r = abs(a(p(i),k))/s(p(i)) if (r .gt. rmax) then j = i rmax = r endif 4 continue pk = p(j) p(j) = p(k) p(k) = pk do 6 i=k+1,n z = a(p(i),k)/a(p(k),k) a(p(i),k) = z do 5 j=k+1,n a(p(i),j) = a(p(i),j) - z*a(p(k),j) 5 continue 6 continue 7 continue c return end c subroutine dot(n,x,y,t) c c compute dot product t= c dimension x(n),y(n) c t = 0.0 do 2 j=1,n t = t + x(j)*y(j) 2 continue c return end c subroutine prod(n,a,x,y) c c compute y=ax c dimension a(n,n),x(n),y(n) c do 2 i=1,n y(i) = 0.0 do 2 j=1,n y(i) = y(i)+a(i,j)*x(j) 2 continue c return end c subroutine store(n,x,y) c c put x into y c dimension x(n),y(n) c do 2 j=1,n y(j) = x(j) 2 continue c return end c subroutine norm(n,x,t) c c compute t= max-norm of x c real x(n),max c t = 0.0 do 2 j=1,n t = max(t,abs(x(j))) 2 continue c return end c subroutine normal(n,x) c c normalize x with max norm c dimension x(n) c call norm(n,x,t) do 2 j=1,n x(j) = x(j)/t 2 continue c return end c subroutine solve(ia,n,a,p,b,x) c c solve for x c dimension a(ia,n),b(n),x(n) integer p(n) c do 2 k=1,n-1 do 2 i=k+1,n b(p(i)) = b(p(i)) - a(p(i),k)*b(p(k)) 2 continue c do 4 i=n,1,-1 sum = b(p(i)) do 3 j=i+1,n sum = sum - a(p(i),j)*x(j) 3 continue x(i) = sum/a(p(i),i) 4 continue c return end @//E*O*F ipoweracc.f// chmod u=rw,g=,o= ipoweracc.f echo x - ex1s52.f sed 's/^@//' > "ex1s52.f" <<'@//E*O*F ex1s52.f//' 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 5.2 c c Example of Schur factorization c c c file: ex1s52.f c parameter (n=3) dimension a(n,n),u(n,n),t(n,n),v(n,n) data (a(1,j),j=1,n) /361.,123.,-180./ data (a(2,j),j=1,n) /148.,414.,-240./ data (a(3,j),j=1,n) /-92.,169.,65./ data (u(1,j),j=1,n) /0.36,0.48,0.80/ data (u(2,j),j=1,n) /0.48,0.64,-0.60/ data (u(3,j),j=1,n) /0.80,-0.60,0.0/ c print * print *,' Schur factorization example' print *,' Section 5.2, Kincaid-Cheney' print * c print *,' Matrix A' call prtmtx(n,n,a) c print *,' Matrix U' call prtmtx(n,n,u) c call mult(n,u,u,v) print *,' Matrix (U*)U' call prtmtx(n,n,v) c call mult(n,u,a,t) print *,' Matrix UA' call prtmtx(n,n,t) c call mult(n,t,u,v) print *,' Matrix UA(U*)' call prtmtx(n,n,v) c stop end c subroutine prtmtx(n,n,a) c c print the matrix c dimension a(n,n) c do 2 i=1,n print 3,(a(i,j),j=1,n) 2 continue print * c return 3 format (2x,3(e13.6,2x)) end c subroutine mult(n,a,b,c) c c compute matrix product c=ab c dimension a(n,n),b(n,n),c(n,n) c do 3 i=1,n do 3 j=1,n x = 0.0 do 2 k=1,n x = x + a(i,k)*b(k,j) 2 continue c(i,j)=x 3 continue c return end @//E*O*F ex1s52.f// chmod u=rw,g=,o= ex1s52.f echo x - qrshif.f sed 's/^@//' > "qrshif.f" <<'@//E*O*F qrshif.f//' 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 5.3 c c Example of modified Gram-Schmidt algorithm c c c file: qrshif.f c parameter (n=4,m=10) dimension a(n,n),q(n,n),r(n,n) data (a(1,j),j=1,n)/1.,2.,3.,4./ data (a(2,j),j=1,n)/5.,6.,7.,8./ data (a(3,j),j=1,n)/0.,9.,10.,11./ data (a(4,j),j=1,n)/0.,0.,12.,13./ c print * print *,' Modified Gram-Schimdt example' print *,' Section 5.3, Kincaid-Cheney' print * c do 4 k=1,m print *,' Matrix A' call prtmtx(n,n,n,a) z = a(n,n) do 2 i=1,n a(i,i)= a(i,i) - z 2 continue call mgs(a,q,r,n,n) print *,' Matrix Q' call prtmtx(n,n,n,q) print *,' Matrix R' call prtmtx(n,n,n,r) call mult(r,q,a,n) do 3 i=1,n a(i,i)=a(i,i) + z 3 continue 4 continue c stop end c subroutine prtmtx(n,n,n,a) c c print array a c dimension a(n,n) c do 2 i=1,n print 3,(a(i,j),j=1,n) 2 continue print * c return 3 format(2x,4(e13.6,2x)) end c subroutine mgs(a,q,t,m,n) c c modified Gram-Schimdt c dimension a(m,n),q(m,n),t(n,n) c do 3 j=1,n do 2 i=1,n t(i,j)=0. 2 continue do 3 i=1,m q(i,j)=a(i,j) 3 continue do 7 k=1,n z=0. do 4 i=1,m z=z+q(i,k)**2 4 continue t(k,k)=sqrt(z) do 5 i=1,m q(i,k)=q(i,k)/t(k,k) 5 continue do 7 j=k+1,n z=0. do 6 i=1,m z=z+q(i,j)*q(i,k) 6 continue t(k,j)=z do 7 i=1,m q(i,j)=q(i,j)-t(k,j)*q(i,k) 7 continue c return end c subroutine mult (a,b,c,n) c c compute matrix product c=ab c dimension a(n,n),b(n,n),c(n,n) c do 3 i=1,n do 3 j=1,n x=0.0 do 2 k=1,n x=x+a(i,k)*b(k,j) 2 continue 3 c(i,j)=x c return end @//E*O*F qrshif.f// chmod u=rw,g=,o= qrshif.f echo x - ex1s53.f sed 's/^@//' > "ex1s53.f" <<'@//E*O*F ex1s53.f//' 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 5.3 c c Example of QR-factorization using Householder transformations c c c file: ex1s53.f c parameter (n=3,ia=4,m=4) dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) data (A(i,1),i=1,m) /63.,42.,0.,126./ data (A(i,2),i=1,m) /41.,60.,-28.,82./ data (A(i,3),i=1,m) /-88.,51.,56.,-71./ c print * print *,' QR-factorization using Householder transformations' print *,' Section 5.3, Kincaid-Cheney' print * print *,' Matrix A' call prtmtx(n,ia,m,A) c call QRfac(n,ia,m,A,V,U,Q,W) c print *,' Matrix Q' call prtmtx(m,ia,m,Q) call mult(n,ia,m,Q,A) print *,' Check A = QR' call prtmtx(n,ia,m,A) c stop end c subroutine QRfac(n,ia,m,A,V,U,Q,W) c c QR-factorization c dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) c call setoI(m,ia,m,Q) do 2 k=1,n mm = m - k + 1 nn = n - k + 1 call setoI(m,ia,m,W) call UtimesA(nn,ia,mm,A(k,k),V(k),U(k,k),W(k,k)) print *,' Matrix A' call prtmtx(n,ia,m,A) print *,' Matrix W' call prtmtx(m,ia,m,W) call mult(m,ia,m,W,Q) print *,' Matrix Q' call prtmtx(m,ia,m,Q) 2 continue call trans(m,ia,m,Q) c return end c subroutine setoI(n,ia,m,Q) c c set Q to I c dimension Q(ia,n) c do 3 i=1,m do 2 j=1,n Q(i,j) = 0.0 2 continue Q(i,i) = 1.0 3 continue c return end c subroutine prtmtx(n,ia,m,A) c c print array A c dimension A(ia,n) c do 2 i=1,m print 3,(A(i,j),j=1,n) 2 continue print * c return 3 format(2x,4(e13.6,2x)) end c subroutine UtimesA(n,ia,m,A,V,U,W) c c compute UA c dimension A(ia,n),U(ia,m),V(m),W(ia,m) c call findV(n,ia,m,A,V) print *,' V vector' call prtmtx(1,ia,m,V) call formU(m,ia,V,U) print *,' Matrix U' call prtmtx(m,ia,m,U) call mult(n,ia,m,U,A) call formW(m,ia,U,W) c return end c subroutine findV(n,ia,m,A,V) c c determine vector V c dimension A(ia,n),V(m) c beta = - sqrt(prod(m,A(1,1),A(1,1))) V(1) = A(1,1) - beta dem = V(1)*V(1) + prod(m-1,A(2,1),A(2,1)) alpha = sqrt(2.0/dem) V(1) = alpha*(A(1,1) - beta) do 2 i=2,m V(i) = alpha*A(i,1) 2 continue c return end c function prod(n,x,y) dimension x(n),y(n) c sum = 0.0 do 2 i=1,n sum = sum + x(i)*y(i) 2 continue prod = sum c return end c subroutine formU(m,ia,V,U) c c compute the unitary factor U c dimension U(ia,m),V(m) c do 3 i=1,m do 2 j=1,m U(i,j) = - V(i)*V(j) 2 continue U(i,i) = 1.0 + U(i,i) 3 continue c return end c subroutine formW(m,ia,U,W) c c store U in W c dimension U(ia,m),W(ia,m) c do 2 i=1,m do 2 j=1,m W(i,j)=U(i,j) 2 continue c return end c subroutine trans(n,ia,m,A) c c compute transpose of the matrix A c dimension A(ia,n),C(4,3) c do 2 i=1,m do 2 j=1,n C(i,j) = A(j,i) 2 continue c do 3 i=1,m do 3 j=1,n A(i,j) = C(i,j) 3 continue c return end c subroutine mult(n,ia,m,B,A) c c compute the matrix-matrix product c dimension B(ia,m),A(ia,n),C(4,3) c do 3 j=1,n do 3 i=1,m sum = 0.0 do 2 k=1,m sum = sum + B(i,k)*A(k,j) 2 continue C(i,j) = sum 3 continue c do 4 i=1,m do 4 j=1,n A(i,j) = C(i,j) 4 continue c return end @//E*O*F ex1s53.f// chmod u=rw,g=,o= ex1s53.f echo x - ex2s55.f sed 's/^@//' > "ex2s55.f" <<'@//E*O*F ex2s55.f//' c c Numerical Analysis: c The Mathematics of Scientific Computing c D.R. Kincaid & E.W. Cheney c Brooks/Cole., 1990 c c Section 5.5 c c Example of QR-factorization to reduce given matrix to upper c Hessenberg form c c c file: ex2s55.f c parameter (n=4,ia=4,m=4,itmax=10) dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) data (A(i,1),i=1,m) /1.,-6.,0.,0./ data (A(i,2),i=1,m) /-5.,8.5556,-5.9182,0./ data (A(i,3),i=1,m) /-1.5395,0.10848,-2.1689,-0.14276/ data (A(i,4),i=1,m) /-1.2767,3.6986,-1.1428,3.6133/ c print * print *,' QR-factorization example' print *,' Section 5.5, Kincaid-Cheney' print * print *,' Matrix A' call prtmtx(n,ia,m,A) c do 2 k=1,itmax call QRfac(n,ia,m,A,V,U,Q,W) call scale(n,ia,m,Q,A) print *,' Matrix Q' call prtmtx(m,ia,m,Q) print *,' Matrix R' call prtmtx(n,ia,m,A) call copy(m,ia,m,Q,W) call mult(n,ia,m,A,W) c call mult(n,ia,m,Q,A) print *,' Check A = QR' call prtmtx(n,ia,m,A) c call copy(n,ia,m,W,A) print *,' Matrix A(k+1)' call prtmtx(n,ia,m,A) 2 continue stop end c subroutine QRfac(n,ia,m,A,V,U,Q,W) c c QR-factorization c dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) c call setoI(m,ia,m,Q) do 2 k=1,3 mm = m - k + 1 nn = n - k + 1 call setoI(m,ia,m,W) call UtimesA(nn,ia,mm,A(k,k),V(k),U(k,k),W(k,k)) print *,' Matrix A' call prtmtx(n,ia,m,A) print *,' Matrix W' call prtmtx(m,ia,m,W) call mult(m,ia,m,W,Q) print *,' Matrix Q' call prtmtx(m,ia,m,Q) 2 continue call trans(m,ia,m,Q) c return end c subroutine setoI(n,ia,m,Q) c c set Q to I c dimension Q(ia,n) c do 3 i=1,m do 2 j=1,n Q(i,j) = 0.0 2 continue Q(i,i) = 1.0 3 continue c return end c subroutine prtmtx(n,ia,m,A) c c print array A c dimension A(ia,n) c do 2 i=1,m print 3,(A(i,j),j=1,n) 2 continue print * c return 3 format(2x,4(e13.6,2x)) end c subroutine UtimesA(n,ia,m,A,V,U,W) c c compute UA c dimension A(ia,n),U(ia,m),V(m),W(ia,m) c call findV(n,ia,m,A,V) print *,' V vector' call prtmtx(1,ia,m,V) call formU(m,ia,V,U) print *,' Matrix U' call prtmtx(m,ia,m,U) call mult(n,ia,m,U,A) call formW(m,ia,U,W) c return end c subroutine findV(n,ia,m,A,V) c c determine vector V c dimension A(ia,n),V(m) c beta = - sqrt(prod(m,A(1,1),A(1,1))) V(1) = A(1,1) - beta dem = V(1)*V(1) + prod(m-1,A(2,1),A(2,1)) alpha = sqrt(2.0/dem) V(1) = alpha*(A(1,1) - beta) do 2 i=2,m V(i) = alpha*A(i,1) 2 continue c return end c function prod(n,x,y) 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 formU(m,ia,V,U) c c compute the unitary factor U c dimension U(ia,m),V(m) c do 3 i=1,m do 2 j=1,m U(i,j) = - V(i)*V(j) 2 continue U(i,i) = 1.0 + U(i,i) 3 continue c return end c subroutine formW(m,ia,U,W) c c store U in W c dimension U(ia,m),W(ia,m) c do 2 i=1,m do 2 j=1,m W(i,j)=U(i,j) 2 continue c return end c subroutine trans(n,ia,m,A) c c compute transpose of the matrix A c dimension A(ia,n),C(4,4) c do 2 i=1,m do 2 j=1,n C(i,j) = A(j,i) 2 continue c do 3 i=1,m do 3 j=1,n A(i,j) = C(i,j) 3 continue c return end c subroutine mult(n,ia,m,B,A) c c compute the matrix-matrix product c dimension B(ia,m),A(ia,n),C(4,4) c do 3 j=1,n do 3 i=1,m sum = 0.0 do 2 k=1,m sum = sum + B(i,k)*A(k,j) 2 continue C(i,j) = sum 3 continue c do 4 i=1,m do 4 j=1,n A(i,j) = C(i,j) 4 continue c return end c subroutine copy(n,ia,m,A,B) c c store A in B c dimension A(ia,n),B(ia,n) c do 2 i=1,m do 2 j=1,n B(i,j) = A(i,j) 2 continue c return end c subroutine scale(n,ia,m,Q,R) c c scale the matrix Q & R c dimension Q(ia,m),R(ia,n),D(4) c do 2 k=1,m D(k) = sign(1.0,R(k,k)) 2 continue do 3 j=1,n do 3 i=1,m Q(i,j) = Q(i,j)*D(i) 3 continue do 4 i=1,m do 4 j=1,n R(i,j) = D(j)*R(i,j) 4 continue c return end @//E*O*F ex2s55.f// chmod u=rw,g=,o= ex2s55.f echo x - ex3s55.f sed 's/^@//' > "ex3s55.f" <<'@//E*O*F ex3s55.f//' 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 5.5 c c Example of shifted QR-factorization c c c file: ex3s55.f c parameter (ia=4,m=4,n=4,itmax=4) dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) data (A(1,j),j=1,m) /1.,2.,3.,4./ data (A(2,j),j=1,m) /4.,5.,6.,7./ data (A(3,j),j=1,m) /2.,1.,5.,0./ data (A(4,j),j=1,m) /4.,2.,1.,0./ c print * print *,' Shifted QR-factorization example' print *,' Section 5.5, Kincaid-Cheney' print * print *,' Matrix A' call prtmtx(n,ia,m,A) c call hess(n,ia,A,V,U) print *,' After Hess' c do 2 nn=n,2,-1 print *,' Matrix A' call prtmtx(nn,ia,nn,A) call submtx(nn,ia,nn,A,itmax,Q,U,W,V) print *,' Eigenvalue=',A(nn,nn) 2 continue print *,' Eigenvalue=',A(1,1) stop end c subroutine submtx(n,ia,m,A,itmax,Q,U,W,V) c c shifted QR on deflated matrix c dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) c do 2 k=1,itmax print *,' k =',k z = A(n,n) call shiftA(n,ia,m,A,z) print *,' Shifted A-zI' call prtmtx(n,ia,m,A) call QRfac(n,ia,m,A,V,U,Q,W) c call scale(n,ia,m,Q,A) print *,' Matrix Q' call prtmtx(m,ia,m,Q) print *,' Matrix R' call prtmtx(n,ia,m,A) c call copy(n,ia,m,A,W) call mult(n,ia,m,Q,W) print *,' Check A = QR + zI' call unshiftA(n,ia,m,W,z) call prtmtx(n,ia,m,W) c c A = R*Q + zI c call mult(n,ia,m,A,Q) call copy(n,ia,m,Q,A) call unshiftA(n,ia,m,A,z) print *,' Matrix A(k+1)' call prtmtx(n,ia,m,A) 2 continue c return end c subroutine shiftA(n,ia,m,A,z) c c form shifted matrix c dimension A(ia,n) c do 2 i=1,n a(i,i) = a(i,i) - z 2 continue c return end c subroutine unshiftA(n,ia,m,A,z) c c form unshifted matrix c dimension A(ia,n) c do 2 i=1,n a(i,i) = a(i,i) + z 2 continue c return end c subroutine hess(n,ia,A,V,U) c c compute Hessenberg matrix of A c dimension A(ia,n),U(ia,n),V(n) c do 2 k=1,n-2 nn = n - k call findU(nn,ia,A(k+1,k),V,U) call mult(nn+1,ia,nn,U,A(k+1,k)) print *,' Matrix UA' call prtmtx(n,ia,n,A) call trans(nn,ia,nn,U) call mult(nn,ia,n,A(1,k+1),U) call copy(nn,ia,n,U,A(1,k+1)) print *,' Matrix UAU*' call prtmtx(n,ia,n,A) 2 continue c return end c subroutine QRfac(n,ia,m,A,V,U,Q,W) c c compute QR factorization c dimension A(ia,n),Q(ia,m),U(ia,m),W(ia,m),V(m) c call setoI(m,ia,m,Q) do 2 k=1,n-1 mm = m - k + 1 nn = n - k + 1 call setoI(m,ia,m,W) call UtimesA(nn,ia,mm,A(k,k),V(k),U(k,k),W(k,k)) print *,' Matrix A' call prtmtx(n,ia,m,A) print *,' Matrix W' call prtmtx(m,ia,m,W) call mult(m,ia,m,W,Q) print *,' Matrix Q*' call prtmtx(m,ia,m,Q) 2 continue call trans(m,ia,m,Q) c return end c subroutine setoI(n,ia,m,Q) c c set Q to I c dimension Q(ia,n) c do 3 i=1,m do 2 j=1,n Q(i,j) = 0.0 2 continue Q(i,i) = 1.0 3 continue c return end c subroutine prtmtx(n,ia,m,A) c c print array A c dimension A(ia,n) c do 2 i=1,m print 3,(A(i,j),j=1,n) 2 continue print * c return 3 format(2x,4(e13.6,2x)) end c subroutine UtimesA(n,ia,m,A,V,U,W) c c compute U times A, store into A c dimension A(ia,n),U(ia,m),V(m),W(ia,m) c call findV(n,ia,m,A,V) print *,' V vector' call prtmtx(1,ia,m,V) call formU(m,ia,V,U) print *,' Matrix U' call prtmtx(m,ia,m,U) call mult(n,ia,m,U,A) call formW(m,ia,U,W) c return end c subroutine findU(n,ia,A,V,U) c c find matrix U c dimension A(ia,n),U(ia,n),V(n) c call findV(n,ia,n,A,V) print *,' V vector' call prtmtx(1,ia,n,V) call formU(n,ia,V,U) print *,' Matrix U' call prtmtx(n,ia,n,U) c return end c subroutine findV(n,ia,m,A,V) c c determine vector V c dimension A(ia,n),V(m) c beta = - sqrt(prod(m,A(1,1),A(1,1))) V(1) = A(1,1) - beta dem = V(1)*V(1) + prod(m-1,A(2,1),A(2,1)) alpha = sqrt(2.0/dem) V(1) = alpha*(A(1,1) - beta) do 2 i=2,m V(i) = alpha*A(i,1) 2 continue print *,' beta=',beta,' alpha=',alpha c return end c function prod(n,x,y) 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 formU(m,ia,V,U) c c form matrix U c dimension U(ia,m),V(m) c do 3 i=1,m do 2 j=1,m U(i,j) = - V(i)*V(j) 2 continue U(i,i) = 1.0 + U(i,i) 3 continue c return end c subroutine formW(m,ia,U,W) c c form matrix W c dimension U(ia,m),W(ia,m) c do 2 i=1,m do 2 j=1,m W(i,j)=U(i,j) 2 continue c return end c subroutine trans(n,ia,m,A) c c compute transpose of the matrix A c dimension A(ia,n),C(10,10) c do 2 i=1,m do 2 j=1,n C(i,j) = A(j,i) 2 continue do 3 i=1,m do 3 j=1,n A(i,j) = C(i,j) 3 continue c return end c subroutine mult(n,ia,m,B,A) c c compute the matrix-matrix product c dimension B(ia,m),A(ia,n),C(10,10) c do 3 j=1,n do 3 i=1,m sum = 0.0 do 2 k=1,m sum = sum + B(i,k)*A(k,j) 2 continue C(i,j) = sum 3 continue do 4 i=1,m do 4 j=1,n A(i,j) = C(i,j) 4 continue c return end c subroutine copy(n,ia,m,A,B) c c copy A into B c dimension A(ia,n),B(ia,n) c do 2 i=1,m do 2 j=1,n B(i,j) = A(i,j) 2 continue c return end c subroutine scale(n,ia,m,Q,R) c c scale Q & R (positive diagonals in R) c dimension Q(ia,m),R(ia,n),D(4) c do 2 k=1,m D(k) = sign(1.0,R(k,k)) 2 continue do 3 j=1,n do 3 i=1,m Q(i,j) = Q(i,j)*D(i) 3 continue do 4 i=1,m do 4 j=1,n R(i,j) = D(j)*R(i,j) 4 continue c return end @//E*O*F ex3s55.f// chmod u=rw,g=,o= ex3s55.f echo x - coef.f sed 's/^@//' > "coef.f" <<'@//E*O*F coef.f//' 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.1 c c Computing coefficients in the Newton Form of a polynomial c c c file: coef.f c parameter (n=3) dimension x(0:n),y(0:n),c(0:n) data (x(i),i=0,n) /5.0,-7.0,-6.0,0.0/ data (y(i),i=0,n) /1.0,-23.0,-54.0,-954.0/ c print * print *,' Coefficients in the Newton form of a polynomial' print *,' Section 6.1, Kincaid-Cheney' print * c c(0) = y(0) print 4,0,c(0) do 3 k=1,n d = x(k) - x(k-1) u = c(k-1) do 2 i=k-2,0,-1 u = u*(x(k) - x(i)) + c(i) d = d*(x(k) - x(i)) 2 continue c(k) = (y(k) - u)/d print 4,k,c(k) 3 continue c 4 format (1x,'c(',i1,') =',2x,e13.6) stop end @//E*O*F coef.f// chmod u=rw,g=,o= coef.f echo x - fft.f sed 's/^@//' > "fft.f" <<'@//E*O*F fft.f//' 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.12 c c Example of Fast Fourier Transform method c c c file: fft.f c parameter (m=4,n3=15) complex z(0:n3),c(0:n3),d(0:n3),w,f,u,v c print * print *,' Fast Fourier Transform example' print *,' Section 6.12, Kincaid-Cheney' print * c nn=2**m m2=nn/2 pi=4.0*atan(1.0) y=2.0*pi/float(nn) w=cmplx(0.0 , -y) w=cexp(w) z(0)=cmplx(1.0,0.0) c(0)=f(0.0) do 2 k=1,nn-1 z(k)=w*z(k-1) c(k)=f(y*float(k)) 2 continue print *,' Initial c-vector' print 5,c print * do 4 n=0,m-1 print *,' n = ',n print 5,c print * mn=2**(m-n-1) n2=2**(n-1) do 3 k=0,mn-1 do 3 j=0,n2 u=c(2*n2*k+j) v=z(j*mn)*c(2*n2*k+j+m2) d(4*n2*k+j)=(u+v)*cmplx(0.5,0.0) d(4*n2*k+j+2*n2)=(u-v)*cmplx(0.5,0.0) 3 continue do 4 j=0,nn-1 c(j)=d(j) 4 continue print *,' Final c-vector' print 5,c c 5 format (2x,16( '(',e13.6,',',e13.6,')',4x)) stop end c function f(x) complex f,z,d u=cos(x) v=sin(x) z=cmplx(u,v) d=cmplx(16.0,0.0) do 2 k=1,15 y=float(k) d=d*z + cmplx(16.0-y,0.0) 2 continue f=d c return end @//E*O*F fft.f// chmod u=rw,g=,o= fft.f echo x - adapta.f sed 's/^@//' > "adapta.f" <<'@//E*O*F adapta.f//' 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 max(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 max(f,t(i-1),t(i),c(i),d(i)) call max(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 max(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 @//E*O*F adapta.f// chmod u=rw,g=,o= adapta.f echo x - ex1s71.f sed 's/^@//' > "ex1s71.f" <<'@//E*O*F ex1s71.f//' 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 7.1 c c Computing the derivative of a function c using forward difference formula c c c file: ex1s71.f c parameter (M=26) f(x) = atan(x) c print * print *,' Derivative approximations: forward difference formula' print *,' Section 7.1, Kincaid-Cheney' print * print 3,'k','h','F2','F1','d','r' c h = 1.0 a = sqrt(2.0) f1 = f(a) c do 2 k=0,M F2 = f(a+h) d = F2 - F1 r = d/h print 4,k,h,F2,F1,d,r h = 0.5*h 2 continue c 3 format(a4,a9,a16,a15,a14,a15) 4 format(1x,i3,5(2x,e13.6)) stop end @//E*O*F ex1s71.f// chmod u=rwx,g=x,o=x ex1s71.f echo x - ex2s71.f sed 's/^@//' > "ex2s71.f" <<'@//E*O*F ex2s71.f//' 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 7.1 c c Derivative approximation c using central difference formula c c c file: ex2s71.f c parameter (M=26) f(x) = atan(x) c print * print *,' Derivative approximation: central difference' print *,' Section 7.1, Kincaid-Cheney' print * print 3,'k','h','F2','F1','d','r' c a = sqrt(2.0) h = 1.0 c do 2 k=0,M F2 = f(a+h) F1 = f(a-h) d = F2 - F1 r = d/(2.0*h) print 4,k,h,F2,F1,d,r h = 0.5*h 2 continue c 3 format(a4,a9,a16,a15,a14,a15) 4 format(1x,i3,5(2x,e13.6)) stop end @//E*O*F ex2s71.f// chmod u=rwx,g=x,o=x ex2s71.f echo x - ex5s71.f sed 's/^@//' > "ex5s71.f" <<'@//E*O*F ex5s71.f//' 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 7.1 c c Computing the derivative using Richardson extrapolation c c c file: ex5s71.f c parameter (M=30) dimension d(0:M),r(0:M) f(x) = atan(x) c print * print *,' Derivative approximation: Richardson extrapolation' print *,' Section 7.1, Kincaid-Cheney' print * print 5,'k','d(k)','r(k)' c a = sqrt(2.0) h = 1.0 c do 2 k=1,M h = 0.5*h d(k) = (f(a+h) - f(a-h))/(2.0*h) 2 continue c do 3 k=1,M r(k) = d(k) + (d(k) - d(k-1))/3.0 3 continue c do 4 k=1,M print 6,k,d(k),r(k) 4 continue c 5 format(a6,a12,a14) 6 format(1x,i5,2x,2(e13.6,2x)) stop end @//E*O*F ex5s71.f// chmod u=rwx,g=x,o=x ex5s71.f echo x - ex6s71.f sed 's/^@//' > "ex6s71.f" <<'@//E*O*F ex6s71.f//' 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 7.1 c c Generate a triangular array of derivatives c using Richardson interpolation c c c file: ex6s71.f c parameter (M=4) dimension d(0:M,0:M) f(x) = atan(x) phi(x,h) = (f(x+h) - f(x-h))/(2.0*h) c print * print *,' Complete Richardson extrapolation' print *,' Section 7.1, Kincaid-Cheney' print * print 6,'n','D(n,0)','D(n,1)','D(n,2)','D(n,3)','D(n,4)' c h = 1.0 a = sqrt(2.0) c do 2 n=0,M d(n,0) = phi(a,h/(2**n)) 2 continue c do 4 k=1,M do 3 n=k,M d(n,k) = d(n,k-1) + (d(n,k-1) - d(n-1,k-1))/((4**k) - 1.) 3 continue 4 continue c do 5 n=0,M print 7,n,(d(n,k),k=0,n) 5 continue c 6 format(a6,a13,4(a14)) 7 format(1x,i5,2x,5(e13.6,2x)) stop end @//E*O*F ex6s71.f// chmod u=rwx,g=x,o=x ex6s71.f echo x - gauss5.f sed 's/^@//' > "gauss5.f" <<'@//E*O*F gauss5.f//' 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 7.3 c c Example of 5-point Gaussian integration formula c to integrate sqrt(x) over [0,1] c c c file: gauss5.f c dimension w(0:3),x(0:3) real a,b,c,d double precision x,w,S,f,sqrt,tmp,tmn data x(0) /0.0d0/ data x(1) /0.53846 93101 05683d0/ data x(2) /0.90617 98459 38664d0/ data w(0) /0.56888 88888 88889d0/ data w(1) /0.47862 86704 99366d0/ data w(2) /0.23692 68850 56189d0/ data a,b,c,d/0.,1.,-1.,1./ c f(x) = sqrt(x) c print * print *,' 5-point Gaussian integration example' print *,' Section 7.3, Kincaid-Cheney' print * c S = 0.0d0 do 2 i=0,3 tmp = ( (b - a)*x(i) + a*d - b*c)/(d - c) tmn = (-(b - a)*x(i) + a*d - b*c)/(d - c) S = S + w(i)*(sqrt(tmp) + sqrt(tmn)) 2 continue S = (b - a)*S/(d - c) print *,' S = ',S c stop end @//E*O*F gauss5.f// chmod u=rwx,g=x,o=x gauss5.f echo x - romberg.f sed 's/^@//' > "romberg.f" <<'@//E*O*F romberg.f//' 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 7.4 c c Computing the integral of a function using Romberg algorithm c c c file: romberg.f c parameter (MM = 4) dimension r(0:MM,0:MM) data a,b/1.0,3.0/ f(x) = 1.0/x c print * print *,' Romberg extrapolation' print *,' Section 7.4, Kincaid-Cheney' print * print 6,'n','R(n,0)','R(n,1)','R(n,2)','R(n,3)','R(n,4)' c h = b-a r(0,0) = 0.5*h*(f(a) + f(b)) c do 4 n = 1,MM h = 0.5*h sum = 0.0 do 2 i = 1,2**(n-1) sum = sum + f(a + real(2*i-1)*h) 2 continue r(n,0) = 0.5*r(n-1,0) + h*sum c do 3 m = 1,MM r(n,m) = r(n,m-1) + (r(n,m-1) - r(n-1,m-1))/((4.**m) - 1.) 3 continue 4 continue c do 5 i=0,MM print 7,i,(r(i,j),j = 0,i) 5 continue c 6 format(a6,a13,4(a14)) 7 format(1x,i5,2x,5(e13.6,2x)) stop end @//E*O*F romberg.f// chmod u=rwx,g=x,o=x romberg.f echo x - adapt.f sed 's/^@//' > "adapt.f" <<'@//E*O*F adapt.f//' 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 7.5 c c adaptive quadrature by simpson rule c c c file: adapt.f c parameter (n = 5) dimension v(6,n) logical iprt data a,b,epsi/0.0,1.0,5.0e-5/ data iprt/.true./ f(x) = sqrt(x) c print * print *,' Adaptive Quadrature' print *,' Section 7.5, Kincaid-Cheney' print * print 4,'k','v(1,n)','v(2,n)','v(3,n)','v(4,n)','v(5,n)','v(6,n)' c c initialize everything, c particularly the first column vector in the stack. c delta = b-a sigma = 0.0 h = delta/2.0 c = 0.5*(a + b) k = 1 abar = f(a) bbar = f(b) cbar = f(c) S = (abar + 4.0*cbar + bbar)*h/3.0 v(1,1) = a v(2,1) = h v(3,1) = abar v(4,1) = cbar v(5,1) = bbar v(6,1) = S if (iprt) print 5,k,(v(i,1),i=1,6) c 2 continue if ( (k .lt. 1) .or. (k .gt. n) ) go to 3 c c take the last column off the stack and process it. c h = 0.5*v(2,k) y = v(1,k) + h ybar = f(y) Star = ( v(3,k) + 4.0*ybar + v(4,k) )*h/3.0 z = v(1,k) + 3.0*h zbar = f(z) SStar = ( v(4,k) + 4.0*zbar + v(5,k) )*h/3.0 c tmp = Star + SStar - v(6,k) if ( abs(tmp) .le. 30.*epsi*h/delta ) then c if the tolerance is being met, add partial integral c and take a new vector from the bottom of the stack. c sigma = sigma + Star + SStar + tmp/15.0 if (iprt) print *,' k=',k,' sigma=',sigma k = k-1 c if (k .le. 0) then print *,' value of integral is ',sigma stop endif else if (k .ge. n) then c if the stack has reached it maximum number of columns (n), c then stop and report a failure. c print *,' method fails' stop endif c c if the tolerance is not being met, subdivide the interval c and create two new vectors in the stack, c one of which overwrites the vector just processed. if (iprt) print *,' split' vbar = v(5,k) v(2,k) = h v(5,k) = v(4,k) v(4,k) = ybar v(6,k) = Star if (iprt) print 5,k,(v(i,k),i=1,6) c k = k+1 v(1,k) = v(1,k-1) + 2.0*h v(2,k) = h v(3,k) = v(5,k-1) v(4,k) = zbar v(5,k) = vbar v(6,k) = SStar if (iprt) print 5,k,(v(i,k),i=1,6) endif c c having created two new vectors in the stack, c pick off the last one and start the proces anew. c go to 2 c 3 continue c 4 format (a2,a12,a14,a14,a14,a14,a14) 5 format (i2,1x,6(e13.6,1x)) stop end @//E*O*F adapt.f// chmod u=rwx,g=x,o=x adapt.f echo x - taylor.f sed 's/^@//' > "taylor.f" <<'@//E*O*F taylor.f//' 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 8.2 c c Solving the initial value problem using Taylor Series c c c file: taylor.f c data M,h,t,x/200,0.01,-1.0,3.0/ c print * print *,' Taylor series method (order 4) ' print *,' Section 8.1, Kincaid-Cheney' print * print 3,'k','t','x' print 4,0,t,x c do 2 k=1,M x1 = cos(t) - sin(x) + t**2.0 x2 = -sin(t) - x1*cos(x) + 2.0*t x3 = -cos(t) - x2*cos(x) + (x1**2.0)*sin(x) + 2.0 x4 = sin(t) + ((x3)**3.0 -x3)*cos(x) + 3.0*x1*x2*sin(x) x = x + h*(x1 + (h/2.)*(x2 + (h/6.)*(x3 + (h/24.)*x4))) t = t + h print 4,k,t,x 2 continue c 3 format(a6,a9,a15) 4 format(1x,i5,2x,e13.6,2x,e13.6) stop end @//E*O*F taylor.f// chmod u=rwx,g=x,o=x taylor.f echo x - rk4.f sed 's/^@//' > "rk4.f" <<'@//E*O*F rk4.f//' 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 8.3 c c Solving an initial value problem using Runge-Kutta method of order 4 c c c file: rk4.f c data M/256/, t/1.0/, x/2.0/, h/7.8125e-3/ c print * print *,' Runge-Kutta method of order 4' print *,' Section 8.3, Kincaid-Cheney' print * print 3,'k','t','x','e' c e = abs(u(t) - x) print 4,0,t,x,e h2 = 0.5*h c do 2 k = 1,M f1 = h*f(t,x) f2 = h*f(t + h2,x + 0.5*f1) f3 = h*f(t + h2,x + 0.5*f2) f4 = h*f(t + h,x + f3) x = x + (f1 + 2.0*(f2 + f3) + f4)/6.0 t = t + h e = abs(u(t) - x) print 4,k,t,x,e 2 continue c 3 format(a6,a9,2a16) 4 format(1x,i5,2x,3(e14.7,2x)) stop end c function f(t,x) f = (t*x -x**2.0)/t**2.0 return end c function u(t) u = t/(0.5 + alog(t)) return end @//E*O*F rk4.f// chmod u=rwx,g=x,o=x rk4.f echo x - rkfelberg.f sed 's/^@//' > "rkfelberg.f" <<'@//E*O*F rkfelberg.f//' 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 8.3 c c Runge-Kutta-Fehlberg method c for solving an initial value problem c c c file: rkfelbergg.f c external f dimension a(6),e(6),c(6),d(6,6),p(6) data M/500/, t/1.0/, x/2.0/, h/0.1/, delta/5.0e-5/, b/3.0/ c data (a(i),i=1,6)/ + .118518518518518518518518519, 0., .518986354775828460038986355, + .50613149034201665780613149, -.18, .0363636363636363636363636364/ data (c(i),i=1,6)/ + 0., .25, .375, .923076923076923076923076923, 1., .5/ data (e(i),i=1,6)/ + .00277777777777777777777777778,0.,-.0299415204678362573099415205, + -.0291998936735778841041998937,.02,.0363636363636363636363636364/ data ((d(i,j),j=1,6),i=1,6)/ + 0., 0., 0., 0., 0., 0., + .25, 0., 0., 0., 0., 0., + .09375, 0., .28125, 0., 0., 0., + .879380974055530268548020027,0.,-3.27719617660446062812926718, + 3.32089212562585343650432408, 0., 0., + 2.03240740740740740740740741,0.,-8.,7.17348927875243664717348928, + -.205896686159844054580896686, .45297270955165692007797271, + -.296296296296296296296296296,0.,2., + -1.38167641325536062378167641, 0., -.275/ c print * print *,' Runge-Kutta-Fehlberg method' print *,' Section 8.3, Kincaid-Cheney' print * print 7,'k','t','x','e','h' iflag = 1 c do 6 k=1,M xd = b-t if (abs(xd) .le. h) then iflag = 0 h = xd endif y = x s = t c do 3 i=1,6 sum = 0.0 do 2 j=1,i-1 sum = sum + d(i,j)*p(j) 2 continue p(i) = h*f(t+(c(i)*h),x+sum) 3 continue c sum = 0.0 do 4 i=1,6 sum = sum + a(i)*p(i) 4 continue x = x + sum c est = 0.0 do 5 i=1,6 est = est + e(i)*p(i) 5 continue c t = t + h print 8,k,t,y,est,h c if (iflag .eq. 0) then stop endif c if (abs(est) .ge. delta) then print *,' half h' h = h/2.0 t = s x = y k = k-1 else if (abs(est) .lt. delta/128.0) then print *,' double h' h = 2*h endif endif 6 continue c 7 format(a6,a9,2a15,a15) 8 format(1x,i5,2x,4(e13.6,2x)) stop end c function f(t,x) f = (t*x - x**2.0)/t**2.0 return end @//E*O*F rkfelberg.f// chmod u=rwx,g=x,o=x rkfelberg.f echo x - taysys.f sed 's/^@//' > "taysys.f" <<'@//E*O*F taysys.f//' 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 8.6 c c Taylor series method (order 3) for solving a c system of ordinary differential equations c c c file: taysys.f c data t/1.0/, x/3.0/, y/1.0/, h/-0.1/, m/30/ c print * print *,' Taylor series method (order 3) for a system' print *,' Section 8.6, Kincaid-Cheney' print * print 3,'k','t','x','y' print 4,0,t,x,y c do 2 k = 1,m x1 = x + y**2.0 - t**3.0 y1 = y + x**3.0 + cos(t) x2 = x1 + 2.0*y*y1 - 3.0*(t**2.0) y2 = y1 + 3.0*x1*(x**2.0) - sin(t) x3 = x2 + 2.0*y*y2 + 2.0*(y1**2.0) - 6.0*t y3 = y2 + 6.0*x*(x1**2.0) + 3.0*(x**2.0)*x2 -cos(t) x = x + h*(x1 + (h/2.0)*(x2 + (h/3.0)*(x3))) y = y + h*(y1 + (h/2.0)*(y2 + (h/3.0)*(y3))) t = t + h print 4,k,t,x,y 2 continue c 3 format(a6,a9,2a15) 4 format(1x,i5,2x,3(e13.6,2x)) stop end @//E*O*F taysys.f// chmod u=rwx,g=x,o=x taysys.f echo x - exs91.f sed 's/^@//' > "exs91.f" <<'@//E*O*F exs91.f//' 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 9.1 c c Example of solving a P.D.E using Explicit Method c c c file: exs91.f c parameter (n=10, np1=11) dimension w(0:np1),v(0:np1),error(0:np1) real a,b data M,hk/200,0.005125/ c print * print *,' Explicit method example' print *,' Section 9.1, Kincaid-Cheney' print * print 7 c pi = 4.0*atan(1.0) h = 1.0/real(np1) s = hk/h**2 c do 2 i=0,np1 w(i) = sin(pi*real(i)*h) 2 continue c t = 0.0 print 8,0,t,w(0),w(1),w(10),w(11) print * c do 6 j=1,M v(0) = a(real(j)*hk) v(np1) = b(real(j)*hk) do 3 i = 1,n v(i) = s*w(i-1) + (1.0 - 2.0*s)*w(i) + s*w(i+1) 3 continue t = real(j)*hk if (0 .eq. mod(j,10)) then print 8,j,t,v(0),v(1),v(10),v(11) print * endif c do 4 i=1,n w(i) = v(i) 4 continue c do 5 i=0,np1 x = real(i)*h error(i) = abs(sin(pi*x)/exp(pi*pi*t) - v(i)) 5 continue c if (0 .eq. mod(j,10)) then print *,' norm of error =', vnorm(n,error) print * endif 6 continue c 7 format (3x,'j',8x,'t',13x,'v(0)',11x,'v(1)',6x,'.....', + 5x,'v(10)',10x,'v(11)') 8 format (1x,i3,2x,3(e13.6,2x),5x,2(e13.6,2x)) stop end c real function a(t) a = 0.0 return end c real function b(t) b = 0.0 return end c function vnorm(n,v) dimension v(0:n) c c computes infinity norm of n component vector v c real max c temp = 0.0 do 2 i=0,n temp = max(temp,abs(v(i))) 2 continue vnorm = temp c return end @//E*O*F exs91.f// chmod u=rw,g=,o= exs91.f echo x - exs92.f sed 's/^@//' > "exs92.f" <<'@//E*O*F exs92.f//' 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 9.2 c c Example of solving a P.D.E using implicit method c c c file: exs92.f c parameter (n=10) dimension v(n),a(n),d(n),c(n),error(n) data M,hk/200,0.005125/ c print * print *,' Implicit method example' print *,' Section 9.2, Kincaid-Cheney' print * print 7 c pi = 4.0*atan(1.0) h = 1.0/real(n+1) s = hk/h**2 c do 2 i=1,n v(i) = sin(pi*real(i)*h) d(i) = 1.0 + 2.0*s 2 continue c t = 0.0 print 8,0,t,v(0),v(1),v(10),v(11) print * c do 3 i=1,n-1 c(i) = -s a(i) = -s 3 continue c do 6 j=1,M call tri(n,a,d,c,v,v) t = real(j)*hk if (0 .eq. mod(j,10)) then print 8,j,t,v(0),v(1),v(10),v(11) print * endif c do 5 i=1,n x = real(i)*h error(i) = abs(sin(pi*x)/exp(pi*pi*t) - v(i)) 5 continue c if (0 .eq. mod(j,10)) then print *,' norm of error =', unorm(n,error) print * endif 6 continue c 7 format (3x,'j',8x,'t',13x,'v(0)',11x,'v(1)',6x,'.....', + 5x,'v(10)',10x,'v(11)') 8 format (1x,i3,2x,3(e13.6,2x),5x,2(e13.6,2x)) stop end c subroutine tri(n,a,d,c,b,x) c c solve a tridiagonal system c dimension a(n),d(n),c(n),b(n),x(n) c do 2 i = 2,n xmult = a(i-1)/d(i-1) d(i) = d(i) - xmult*c(i-1) b(i) = b(i) - xmult*b(i-1) 2 continue c x(n) = b(n)/d(n) do 3 i = n-1,1,-1 x(i) = (b(i) - c(i)*x(i+1))/d(i) 3 continue c return end c function unorm(n,u) c c computes infinity norm of n component vector v c dimension u(n) real max c temp = 0.0 do 2 i=1,n temp = max(temp,abs(u(i))) 2 continue unorm = temp c return end @//E*O*F exs92.f// chmod u=rw,g=,o= exs92.f echo x - exs93.f sed 's/^@//' > "exs93.f" <<'@//E*O*F exs93.f//' 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 9.3 c c Example of Finite Difference Method c c c file: exs93.f c parameter (n=18) dimension v(0:n,0:n),x(0:n),y(0:n),u(0:n,0:n) c print * print *,' Finite Difference Method' print *,' Section 9.3, Kincaid-Cheney' print * print *,' error' c m = 200 nn = n+1 h = 1.0/nn do 2 i=0,n+1 x(i) = i*h y(i) = x(i) 2 continue do 3 i=0,n+1 v(i,0)=g(x(i),0) v(i,n+1)=g(x(i),1) v(0,i)=g(0,y(i)) v(n+1,i)=g(1,y(i)) 3 continue do 4 k=1,m do 4 j=1,n do 4 i=1,n v(i,j) = (v(i-1,j) + v(i+1,j) + v(i,j-1) + v(i,j+1))/4.0 4 continue k = 0 do 5 i=0,n do 5 j=0,n u(i,j) = g(x(i),y(j)) err = abs(v(i,j) - u(i,j)) if (0 .eq. mod(k,20)) print 6,err k = k+1 5 continue c 6 format (e13.6) stop end c real function g(x,y) g = (10.0e-4)*sinh(3.0*22.0*x/7.0)*sin(3.0*22.0*y/7.0) return end @//E*O*F exs93.f// chmod u=rw,g=,o= exs93.f echo x - ex3s96.f sed 's/^@//' > "ex3s96.f" <<'@//E*O*F ex3s96.f//' 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 9.6 c c Example of use of method of characteristics c c c file: ex3s96.f c parameter(n=8) dimension x(n),y(n),u(n),p(n),q(n) c print * print *,' Method of Characteristics' print *,' Section 9.6, Kincaid-Cheney' print * print 5 c h = 1.0/(n-1) do 2 j=1,n x(j) = (j-1)*h y(j) = 0.0 u(j) = f(x(j)) p(j) = df(x(j)) q(j) = g(x(j)) print 6,x(j),y(j),p(j),q(j),u(j) 2 continue print * c do 4 i=1,n y(i) = i*h do 3 j=1,n-i x(j) = x(j) + h/2 pp = (p(j) + p(j+1))/2. + (1.0 - h/8.)*(q(j+1) - q(j)) qq = (p(j+1)-p(j) + (2.- h/4.)*(q(j+1) + q(j)))/(4.+ h/2.) u(j) = u(j) + h*(pp+p(j))/4. + h*(qq+q(j))/2. p(j) = pp q(j) = qq print 6,x(j),y(j),p(j),q(j),u(j) 3 continue print * 4 continue c 5 format (8x,'x',14x,'y',14x,'p',14x,'q',14x,'u') 6 format (2x,5(e13.6,2x)) stop end c real function f(x) f = exp(x/4.) return end c real function g(x) g = (-1 + sqrt(2.))*exp(x/4.)/8. return end c real function df(x) df = exp(x/4.0)/4. return end c real function tu(x,y) tu = exp(x/4.)*exp((-1+sqrt(2.))*y/8.) return end @//E*O*F ex3s96.f// chmod u=rw,g=,o= ex3s96.f echo x - mgrid1.f sed 's/^@//' > "mgrid1.f" <<'@//E*O*F mgrid1.f//' 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 9.8 c c multigrid program used on the equation u = f, in one variable. c c parameter m - the number of mesh refinements used in the method c n - 2**m c k - number of iterations in the Gauss-Seidel method c to be used on each grid. c c c file: mgrid1.f c parameter (m=6, k=3, np1=65) dimension v(0:np1),w(0:np1) c c the function f is the right-hand side of the differential equation. c the function g is the solution function. c real f,g f(x) = cos(x) g(x) = -cos(x) + x*(cos(1.0) - 1.0) + 1.0 c print * print *,' Multigrid method example' print *,' Section 9.8, Kincaid-Cheney' print * c c first we define the solution of the difference equations on the c coarsest grid. this grid has three points, two of which are on the c boundary. c h = 0.5 n = 1 v(0) = 0.0 v(2) = 0.0 v(1) = -f(0.5)/8.0 print *,' n =',n,' h =',h print * c c now the main iteration (outer loop) begins. c do 8 i=1,m-1 h = 0.5*h n = 2*n + 1 print *,' n =',n,' h =',h c c the next part is the interpolation phase in which an approximate c solution on a finer grid is manufactured from an approximation on c a coarse grid. c do 2 j=0,n-1 w(2*j) = v(j) 2 continue c do 3 j=1,n-1 w(2*j-1) = 0.5*(v(j-1) + v(j)) 3 continue c do 4 j=0,2*n-2 v(j) = w(j) 4 continue c c next is the Gauss-Seidel iteration on the finer grid (k steps) c do 6 p=1,k do 5 j=1,n v(j) = 0.5*(v(j-1) + v(j+1) - h*h*f(real(j)*h)) 5 continue 6 continue c do 7 j=0,n+1 w(j) = abs(g(real(j)*h) - v(j)) 7 continue c print *,' error norm =',vnorm(n+1,w) print * 8 continue c stop end c function vnorm(n,v) c c infinity norm of vector v c dimension v(0:n) real max c temp = 0.0 do 2 i=0,n temp = max(temp,abs(v(i))) 2 continue vnorm = temp c return end @//E*O*F mgrid1.f// chmod u=rw,g=,o= mgrid1.f echo x - exs98.f sed 's/^@//' > "exs98.f" <<'@//E*O*F exs98.f//' 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 9.8 c c Example of Damping of Errors c c c file: exs98.f c parameter (n=63) dimension v(0:n+1),rho(100) integer p(4) data (p(i),i=1,4) /1,4,7,16/ c print * print *,' Damping of Errors' print *,' Section 9.8, Kincaid-Cheney' print * print 6 c k = 100 do 5 l=1,4 do 2 j=0,n+1 v(j) = sin((j*p(l)*22.0/7.0)/(n+1)) 2 continue do 5 i=1,k do 3 j=1,n v(j) = (v(j-1) + v(j+1))/2.0 3 continue rho(i) = abs(v(1)) do 4 j=2,n if (rho(i) .ge. abs(v(j))) goto 4 rho(i) = abs(v(j)) 4 continue print 7,p(l),i,rho(i) 5 continue c 6 format (2x,'p',5x,'i',8x,'norm of v') 7 format (1x,i2,3x,i3,6x,e13.6) stop end @//E*O*F exs98.f// chmod u=rw,g=,o= exs98.f echo x - mgrid2.f sed 's/^@//' > "mgrid2.f" <<'@//E*O*F mgrid2.f//' 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 9.8 c c complete V-cycle in the multigrid algorithm to solve c uxx = g(x), with u(0)=u(l)=0. c parameter m gives the number of grids to be employed. c parameter k gives the number of iterations to be performed c in each gauss-seidel iteration. c parameter np1 is simply n+1 c c c file: mgrid2.f c parameter (m=7, np1=128) dimension v(m,0:np1), f(m,0:np1), w(0:np1), error(15) integer p real g,true g(x) = cos(x) true(x) = -g(x) + x*(g(1.0) - 1.0) + 1.0 c print * print *,' V-cycle in the multigrid method' print *,' Section 9.8, Kincaid-Cheney' print * c error(1) = 1.0 do 12 k=2,15 c c initialize arrays and c put best available guess into v**m c n = 2**m - 1 h = 1.0/real(n+1) print *,' n =',n,' h =',h c do 1 i=1,m do 1 j=0,np1 f(i,j) = 0.0 v(i,j) = 0.0 1 continue c c put data from the bvp into f**m c do 2 j=1,n f(m,j) = g(real(j)*h) 2 continue c do 6 i=m,2,-1 c c gauss-seidel iteration (k steps) c do 3 p=1,k do 3 j=1,n v(i,j) = 0.5*(v(i,j-1) + v(i,j+1) - (h**2)*f(i,j)) 3 continue c c residual vector c do 4 j=1,n w(j) = f(i,j) - (v(i,j-1) - 2.0*v(i,j) + v(i,j+1))/h**2 4 continue c c apply restriction operator c do 5 j=1,(n-1)/2 f(i-1,j) = w(2*j) 5 continue c h = 2.0*h n = (n-1)/2 print *,' n =',n,' h =',h 6 continue c c solve coarsest (smallest) system exactly c v(1,1) = -g(0.5)/8.0 print *,' Bottom of V-cycle' c c end downward phase of V-cycle c start upward phase of V-cycle c do 10 i=2,m h = 0.5*h n = 2*n + 1 print *,' n =',n,' h =',h c c apply extension operation, i.e. c interpolation up from coarse grid to fine grid c do 7 j=0,(n+1)/2 w(2*j) = v(i-1,j) 7 continue c do 8 j=1,(n+1)/2 w(2*j-1) = 0.5*(v(i-1,j-1) + v(i-1,j)) 8 continue c c add correction to extension operator c do 9 j=0,n+1 v(i,j) = v(i,j) + w(j) 9 continue c c gauss-seidel iteration (k steps) c do 10 p=1,k do 10 j=1,n v(i,j) = 0.5*(v(i,j-1) + v(i,j+1) - (h**2)*f(i,j)) 10 continue c c output phase c print *,' k =',k,' m =',m do 11 j=0,np1 w(j) = true(real(j)*h) - v(m,j) 11 continue error(k) = vnorm(np1,w) print *,' maximun error =', error(k), + ' ratio =', error(k)/error(k-1) c 12 continue c stop end c function vnorm(n,v) c c infinity norm of vectorv v c dimension v(0:n) real max c temp = 0.0 do 2 i=0,n temp = max(temp,abs(v(i))) 2 continue vnorm = temp c return end @//E*O*F mgrid2.f// chmod u=rw,g=,o= mgrid2.f echo x - code-info.tex sed 's/^@//' > "code-info.tex" <<'@//E*O*F code-info.tex//' %15 Mar 91 \documentstyle{article} \topmargin=-.5in \textwidth 6.0in \textheight 8.5in \oddsidemargin .25in \evensidemargin .25in \begin{document} \title{\bf NUMERICAL ANALYSIS \\ MATHEMATICS OF SCIENTIFIC COMPUTING} \smallskip \author{ David Kincaid and Ward Cheney \\ \copyright 1991 Brooks/Cole Publishing Co.\\ ISBN 0--534--13014--3} \date{March 15, 1991} \maketitle \section*{Introduction} The following table lists files containing sample programs based on the pseudocode given in the textbook cited above. They are intended primarily as a learning and teaching aid for use with this book. We believe that these computer routines are coded in a clear and easy-to-understand style. We have intentionally included few comment statements so that students will read the code and study the algorithms---they can add comments as they decipher them. These programs are usable on computer systems with Fortran 77 compilers, from small personal computers to large scientific computing machines. However, they do not contain all of the ``bells-and-whistles'' of robust state-of-the-art software such as may be found in general-purpose scientific libraries. Nevertheless, they are adequate for many small nonpathological problems. \noindent \section*{Installation and Usage} On a Unix system, unpack the file kincaid-cheney.shar as follows \indent {\tt sh kincaid-cheney.shar} \noindent These programs will run as is on any computer with a standard Fortran 77 compiler. However, the statement {\tt data epsi/1.0e-6/} in some routines should be changed to the machine epsilon (single precision roundoff error) for the computer that is to be used. Compile and execute the program on file romberg.f as follows \indent {\tt f77 romberg.f} \indent {\tt a.out} \section*{Availability} For information on the availability of this software, contact either the publisher of the textbook or the authors at the following addresses. \medskip \begin{center} \begin{tabular}{lcl} Brooks/Cole Publishing Co.& &Center for Numerical Analysis\\ 511 Forest Lodge Road& &University of Texas at Austin\\ Pacific Grove, CA 93950--5098& &Austin, TX 78713--8510\\[0.1in] (408) 373--0728& &(512) 471--1242\\ fax: (408) 375--6414& &fax: (512) 471--9038 \\ &&kincaid@cs.utexas.edu \end{tabular} \end{center} \newpage \begin{center} \begin{tabular}{lrl} {\bf File Name} & {\bf Pages} & {\bf Description of Code\quad (Subprogram Names)} \\[0.1in] {\bf \tt elimit.f} & 10 & Example of a slowly converging sequence \\ {\bf \tt sqrt2.f} & 10 & Example of a rapidly converging sequence \\ {\bf \tt nest.f} & 14 & Nested multiplication \\[0.1in] {\bf \tt epsi.f} & 36 & Approximate value of machine precision \\ {\bf \tt depsi.f} & 36 & Approximate value of double precision machine precision \\ {\bf \tt ex2s22.f} & 43--44 & Loss of significance \\ {\bf \tt unstab1.f} & 49 & Example of an unstable sequence \\ {\bf \tt unstab2.f} & 50 & Example of another unstable sequence \\ {\bf \tt instab.f} & 50 & Example of numerical instability \\[0.1in] {\bf \tt ex1s31.f} & 58--59 & Bisection method to find root of $\exp(x)=\sin x$ \quad({\bf bisect}) \\ {\bf \tt ex1s32.f} & 65 & Newton's method example \\ {\bf \tt ex2s32.f} & 68 & Simple Newton's method \\ {\bf \tt ex3s32.f} & 69--70 & Implicit function example \\ {\bf \tt ex1s33.f} & 76 & Secant method example \quad ({\bf f\,}) \\ {\bf \tt ex3s34.f} & 83--84 & Contractive mapping example \\ {\bf \tt ex3s35.f} & 93 & Horner's method example \\ {\bf \tt ex6s35.f} & 94 & Newton's method on a given polynomial \quad({\bf horner}) \\ {\bf \tt ex7s35.f} & 99 & Bairstow's method example \\ {\bf \tt laguerre.f} & 102 & Laguerre's method example \\[0.1in] {\bf \tt forsub.f} & 127 & Forward substitution example \\ {\bf \tt bacsub.f} & 127 & Backward substitution example \\ {\bf \tt pforsub.f} & 128 & Forward substitution for a permuted system \\ {\bf \tt pbacsub.f} & 128 & Backward substitution for a permuted system \\ {\bf \tt genlu.f} & 130 & General LU-factorization example \\ {\bf \tt doolt.f} & 131 & Doolittle's-factorization example \\ {\bf \tt cholsky.f} & 134 & Cholesky-factorization example\\ {\bf \tt bgauss.f} & 143 & Basic Gaussian elimination \\ {\bf \tt pbgauss.f} & 145 & Basic Gaussian elimination with pivoting \\ {\bf \tt gauss.f} & 148 & Gaussian elimination with scaled row pivoting \\ {\bf \tt paxeb.f} & 150 & Solves $Lz=Pb$ and then $Ux=z$ \quad ({\bf gauss}) \\ {\bf \tt yaec.f} & 151 & Solves $U^{T}z=c$ and then $L^{T}Py=z$ \quad ({\bf gauss}) \\ {\bf \tt tri.f} & 155 & Tridiagonal system solver \quad ({\bf tri}) \\ {\bf \tt ex1s45.f} & 173 & Neumann series example \quad ({\bf setI, mult, store, add, prt}) \\ {\bf \tt ex2s45.f} & 175 & Gaussian elimination followed by iterative improvement\\ && \qquad ({\bf residual, gauss, solve}) \\ {\bf \tt ex1s46.f} & 182 & Example of Jacobi and Gauss-Seidel methods \\ {\bf \tt ex2s46.f} & 184 & Richardson method example (with scaling) \\ {\bf \tt jacobi.f} & 186 & Jacobi method example (with scaling) \\ {\bf \tt ex3s46.f} & 190 & Gauss-Seidel method (with scaling)\\ {\bf \tt ex6s46.f} & 200 & Chebyshev acceleration example \quad ({\bf extrap, cheb, vnorm}) \\ {\bf \tt steepd.f} & 207 & Steepest descent method example \quad ({\bf prod, mult}) \\ {\bf \tt cg.f} & 211 & Conjugate gradient method \quad ({\bf prod, residual, mult}) \\ {\bf \tt pcg.f} & 217 & Jacobi preconditioned conjugate gradient method\\ & & \qquad ({\bf prod, residual, mult}) \end{tabular} \end{center} \newpage \begin{center} \begin{tabular}{lrl} {\bf File Name} & {\bf Pages} & {\bf Description of Code \quad (Subprogram Names)} \quad (continued)\\[0.1in] {\bf \tt ex1s51.f} & 231 & Power method example \quad ({\bf dot, prod, store, norm, normal}) \\ {\bf \tt poweracc.f} & 231 & Power method with Aitken acceleration\\ & &\qquad ({\bf dot, prod, store, norm, normal}) \\ {\bf \tt ex2s51.f} & 233 & Inverse power method example \\ & &\qquad ({\bf gauss, dot, prod, store, norm, normal, solve)} \\ {\bf \tt ipoweracc.f} & 233 & Inverse power method with Aitken acceleration \\ & & \qquad ({\bf gauss, dot, prod, store, norm, normal, solve}) \\ {\bf \tt ex1s52.f} & 239 & Schur factorization example \quad ({\bf prtmtx, mult}) \\ {\bf \tt qrshif.f} & 248 & Modified Gram-Schmidt example \quad ({\bf prtmtx, mgs, mult}) \\ {\bf \tt ex1s53.f} & 253--255 & QR-factorization using Householder transformations \\ & &\qquad ({\bf QRfac, setoI, prtmtx, UtimesA, findV, prod,}\\ & &\qquad {\bf formU, formW, trans, mult}) \\ {\bf \tt ex2s55.f} & 272--273 & QR-factorization example \\ & &\qquad ({\bf QRfac, setoI, prtmtx, UtimesA, findV, prod, } \\ & &\qquad {\bf formU, formW, trans, mult, copy, scale}) \\ {\bf \tt ex3s55.f} & 274 & Shifted QR-factorization example \\ & &\qquad ({\bf submtx, shiftA, unshiftA, hess, QRfac, setoI, }\\ & &\qquad {\bf prtmtx, UtimesA, findU, findV, prod, formU, }\\ & &\qquad {\bf formW, Trans, mult, copy, scale}) \\[0.1in] {\bf \tt coef.f} & 280--281 & Coefficients in the Newton form of a polynomial \\ {\bf \tt fft.f} & 419 & Fast Fourier transform example \quad ({\bf f\,}) \\ {\bf \tt adapta.f} & 426--428 & Adaptive approximation example \quad ({\bf f, max}) \\[0.1in] {\bf \tt ex1s71.f} & 431--432 & Derivative approximations: forward difference formula \\ {\bf \tt ex2s71.f} & 434 & Derivative approximation: central difference \\ {\bf \tt ex5s71.f} & 437--438 & Derivative approximation: Richardson extrapolation \\ {\bf \tt ex6s71.f} & 440--441 & Richardson extrapolation \\ {\bf \tt gauss5.f} & 459 & Gaussian five-point quadrature example \\ {\bf \tt romberg.f} & 468 & Romberg extrapolation \\ {\bf \tt adapt.f} & 475 & Adaptive quadrature \\[0.1in] {\bf \tt taylor.f} & 492 & Taylor-series method \\ {\bf \tt rk4.f} & 501--502 & Runge-Kutta method \quad ({\bf f, u}) \\ {\bf \tt rkfelberg.f} & 503--505 & Runge-Kutta-Fehlberg method \quad ({\bf f\,}) \\ {\bf \tt taysys.f} & 526--527 & Taylor series for systems \\[0.1in] {\bf \tt exs91.f} & 576 & Boundary value problem (BVP): Explicit method example\\ & & \qquad ({\bf a, b, vnorm}) \\ {\bf \tt exs92.f} & 582 & BVP: Implicit method example \quad ({\bf tri, unorm}) \\ {\bf \tt exs93.f} & 589 & Finite difference method \quad ({\bf g\,}) \\ {\bf \tt ex3s96.f} & 613 & BVP: Method of characteristics \quad ({\bf f, g, df, tu}) \\ {\bf \tt mgrid1.f} & 624 & Multigrid method example \quad ({\bf vnorm}) \\ {\bf \tt exs98.f} & 625 & Damping of errors \\ {\bf \tt mgrid2.f} & 630 & Multigrid method V-cycle \quad ({\bf vnorm}) \\[0.1in] {\bf \tt code-info.tex}&&This \LaTeX\ file\\ {\bf \tt code-info.tty}&&This tty file \end{tabular} \end{center} \end{document} @//E*O*F code-info.tex// chmod u=rw,g=,o= code-info.tex echo x - code-info.tty sed 's/^@//' > "code-info.tty" <<'@//E*O*F code-info.tty//' NUMERICAL ANALYSIS MATHEMATICS OF SCIENTIFIC COMPUTING David Kincaid and Ward Cheney (c)1991 Brooks/Cole Publishing Co. ISBN 0-534-13014-3 March 15, 1991 Introduction: The following table lists files containing sample programs based on the pseudocode given in the textbook cited above. They are intended primarily as a learning and teaching aid for use with this book. We believe that these computer routines are coded in a clear and easy-to-understand style. We have intentionally included few comment statements so that students will read the code and study the algorithms---they can add comments as they decipher them. These programs are usable on computer systems with Fortran 77 compilers, from small personal computers to large scientific computing machines. However, they do not contain all of the "bells-and-whistles" of robust state-of-the-art software such as may be found in general-purpose scientific libraries. Nevertheless, they are adequate for many small nonpathological problems. Installation and Usage: On a Unix system, unpack the file kincaid-cheney.shar as follows sh kincaid-cheney.shar These programs will run as is on any computer with a standard Fortran 77 compiler. However, the statement data epsi/1.0e-6/ in some routines should be changed to the machine epsilon (single precision roundoff error) for the computer that is to be used. Compile and execute the program on file romberg.f as follows f77 romberg.f a.out Availability: For information on the availability of this software, contact either the publisher or the textbook or the authors at the following addresses. Brooks/Cole Publishing Co. Center for Numerical Analysis 511 Forest Lodge Road University of Texas at Austin Pacific Grove, CA 93950-5098 Austin, TX 78713-8510 (408) 373-0728 (512) 471-1242 fax: (408) 375-6414 fax: (512) 471-9038 kincaid@cs.utexas.edu 1 File Name Pages Description of Code (Subprogram Names) elimit.f 10 Example of a slowly converging sequence sqrt2.f 10 Example of a rapidly converging sequence nest.f 14 Nested multiplication epsi.f 36 Approximate value of machine precision depsi.f 36 Approximate value of double precision machine precision ex2s22.f 43-44 Loss of significance unstab1.f 49 Example of an unstable sequence unstab2.f 50 Example of another unstable sequence instab.f 50 Example of numerical instability ex1s31.f 58-59 Bisection method to find root of exp(x) = sinx (bisect) ex1s32.f 65 Newton's method example ex2s32.f 68 Simple Newton's method ex3s32.f 69-70 Implicit function example ex1s33.f 76 Secant method example (f ) ex3s34.f 83-84 Contractive mapping example ex3s35.f 93 Horner's method example ex6s35.f 94 Newton's method on a given polynomial (horner) ex7s35.f 99 Bairstow's method example laguerre.f 102 Laguerre's method example forsub.f 127 Forward substitution example bacsub.f 127 Backward substitution example pforsub.f 128 Forward substitution for a permuted system pbacsub.f 128 Backward substitution for a permuted system genlu.f 130 General LU-factorization example doolt.f 131 Doolittle's-factorization example cholsky.f 134 Cholesky-factorization example bgauss.f 143 Basic Gaussian elimination pbgauss.f 145 Basic Gaussian elimination with pivoting gauss.f 148 Gaussian elimination with scaled row pivoting paxeb.f 150 Solves Lz = Pb and then Ux = z (gauss) yaec.f 151 Solves UT z = c and then LTPy = z (gauss) tri.f 155 Tridiagonal system solver (tri) ex1s45.f 173 Neumann series example (setI, mult, store, add, prt) ex2s45.f 175 Gaussian elimination followed by iterative improvement (residual, gauss, solve) ex1s46.f 182 Example of Jacobi and Gauss-Seidel methods ex2s46.f 184 Richardson method example (with scaling) jacobi.f 186 Jacobi method example (with scaling) ex3s46.f 190 Gauss-Seidel method (with scaling) ex6s46.f 200 Chebyshev acceleration example (extrap, cheb, vnorm) steepd.f 207 Steepest descent method example (prod, mult) cg.f 211 Conjugate gradient method (prod, residual, mult) pcg.f 217 Jacobi preconditioned conjugate gradient method (prod, residual, mult) 2 File Name Pages Description of Code (Subprogram Names) (continued) ex1s51.f 231 Power method example (dot, prod, store, norm, normal) poweracc.f 231 Power method with Aitken acceleration (dot, prod, store, norm, normal) ex2s51.f 233 Inverse power method example (gauss, dot, prod, store, norm, normal, solve) ipoweracc.f 233 Inverse power method with Aitken acceleration (gauss, dot, prod, store, norm, normal, solve) ex1s52.f 239 Schur factorization example (prtmtx, mult) qrshif.f 248 Modified Gram-Schmidt example (prtmtx, mgs, mult) ex1s53.f 253-255 QR-factorization using Householder transformations (QRfac, setoI, prtmtx, UtimesA, findV, prod, formU, formW, trans, mult) ex2s55.f 272-273 QR-factorization example (QRfac, setoI, prtmtx, UtimesA, findV, prod, formU, formW, trans, mult, copy, scale) ex3s55.f 274 Shifted QR-factorization example (submtx, shiftA, unshiftA, hess, QRfac, setoI, prtmtx, UtimesA, findU, findV, prod, formU, formW, Trans, mult, copy, scale) coef.f 280-281 Coefficients in the Newton form of a polynomial fft.f 419 Fast Fourier transform example (f ) adapta.f 426-428 Adaptive approximation example (f, max) ex1s71.f 431-432 Derivative approximations: forward difference formula ex2s71.f 434 Derivative approximation: central difference ex5s71.f 437-438 Derivative approximation: Richardson extrapolation ex6s71.f 440-441 Richardson extrapolation gauss5.f 459 Gaussian five-point quadrature example romberg.f 468 Romberg extrapolation adapt.f 475 Adaptive quadrature taylor.f 492 Taylor-series method rk4.f 501-502 Runge-Kutta method (f, u) rkfelberg.f 503-505 Runge-Kutta-Fehlberg method (f ) taysys.f 526-527 Taylor series for systems exs91.f 576 Boundary value problem (BVP): Explicit method example (a, b, vnorm) exs92.f 582 BVP: Implicit method example (tri, unorm) exs93.f 589 Finite difference method (g ) ex3s96.f 613 BVP: Method of characteristics (f, g, df, tu) mgrid1.f 624 Multigrid method example (vnorm) exs98.f 625 Damping of errors mgrid2.f 630 Multigrid method V-cycle (vnorm) code-info.tex This LaTeX file code-info.tty This tty file 3 @//E*O*F code-info.tty// chmod u=rw,g=,o= code-info.tty echo Inspecting for damage in transit... temp=/tmp/shar$$; dtemp=/tmp/.shar$$ trap "rm -f $temp $dtemp; exit" 0 1 2 3 15 cat > $temp <<\!!! 36 117 867 elimit.f 33 85 594 sqrt2.f 31 88 602 nest.f 35 99 674 epsi.f 38 108 772 depsi.f 46 153 1043 ex2s22.f 42 109 906 unstab1.f 42 109 892 unstab2.f 31 81 555 instab.f 67 189 1476 ex1s31.f 39 123 920 ex1s32.f 36 101 722 ex2s32.f 44 128 1042 ex3s32.f 54 163 1115 ex1s33.f 30 80 596 ex3s34.f 31 87 618 ex3s35.f 51 157 1168 ex6s35.f 48 140 1114 ex7s35.f 61 198 1448 laguerre.f 40 104 842 forsub.f 40 104 842 bacsub.f 42 116 944 pforsub.f 42 116 939 pbacsub.f 72 194 1705 genlu.f 53 141 1214 doolt.f 49 129 1067 cholsky.f 44 117 1012 bgauss.f 44 120 991 pbgauss.f 79 216 1757 gauss.f 95 246 1990 paxeb.f 94 245 2003 yaec.f 60 151 1187 tri.f 113 261 2048 ex1s45.f 137 360 3079 ex2s45.f 53 159 1123 ex1s46.f 41 128 959 ex2s46.f 55 153 1278 jacobi.f 33 95 724 ex3s46.f 98 254 2185 ex6s46.f 81 216 1770 steepd.f 114 286 2281 cg.f 118 300 2421 pcg.f 113 249 1954 ex1s51.f 116 257 2172 poweracc.f 181 414 3311 ex2s51.f 183 423 3505 ipoweracc.f 80 182 1585 ex1s52.f 112 236 2213 qrshif.f 212 472 3989 ex1s53.f 258 577 4935 ex2s55.f 350 781 6792 ex3s55.f 39 121 866 coef.f 73 151 1460 fft.f 98 256 2115 adapta.f 39 106 786 ex1s71.f 39 101 774 ex2s71.f 43 112 852 ex5s71.f 47 122 979 ex6s71.f 43 155 1060 gauss5.f 48 144 1100 romberg.f 117 426 2789 adapt.f 36 126 916 taylor.f 50 159 1065 rk4.f 102 289 2645 rkfelberg.f 40 161 1099 taysys.f 94 238 1930 exs91.f 101 257 2096 exs92.f 57 152 1184 exs93.f 71 190 1529 ex3s96.f 102 362 2338 mgrid1.f 45 123 983 exs98.f 145 457 3324 mgrid2.f 191 1348 8691 code-info.tex 176 853 7392 code-info.tty 5793 16526 129914 total !!! wc elimit.f sqrt2.f nest.f epsi.f depsi.f ex2s22.f unstab1.f unstab2.f instab.f ex1s31.f ex1s32.f ex2s32.f ex3s32.f ex1s33.f ex3s34.f ex3s35.f ex6s35.f ex7s35.f laguerre.f forsub.f bacsub.f pforsub.f pbacsub.f genlu.f doolt.f cholsky.f bgauss.f pbgauss.f gauss.f paxeb.f yaec.f tri.f ex1s45.f ex2s45.f ex1s46.f ex2s46.f jacobi.f ex3s46.f ex6s46.f steepd.f cg.f pcg.f ex1s51.f poweracc.f ex2s51.f ipoweracc.f ex1s52.f qrshif.f ex1s53.f ex2s55.f ex3s55.f coef.f fft.f adapta.f ex1s71.f ex2s71.f ex5s71.f ex6s71.f gauss5.f romberg.f adapt.f taylor.f rk4.f rkfelberg.f taysys.f exs91.f exs92.f exs93.f ex3s96.f mgrid1.f exs98.f mgrid2.f code-info.tex code-info.tty | sed 's=[^ ]*/==' | diff -b $temp - >$dtemp if [ -s $dtemp ] then echo "Ouch [diff of wc output]:" ; cat $dtemp else echo "No problems found." fi exit 0