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