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