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