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