c     examples/iter/iter_mpi
c
      include "mpif.h"
      integer my_rank,nproc,ierr
      integer n,iroot
      integer status(MPI_STATUS_SIZE)
      real*8 a(4,4),a1(4),b(4,4),x(4),c(4),tmp(4),y(4),tmp1
c     
c     INITIALIZE MPI AND DETERMINE BOTH INDIVIDUAL PROCESSOR #
C     AND THE TOTAL NUMBER OF PROCESSORS
c     
      call init(nproc,my_rank)
c
c     process 0 reads data
c
      if(my_rank.eq.0) then
         open(8,file='in.dat',status='unknown')
         read(8,*) n
         do i=1,n
            read(8,*) (a(i,j),j=1,n)
         enddo
         read(8,*) (c(i),i=1,n)
         write(*,*) (c(i),i=1,n)
      endif
      iroot=0
      call MPI_Bcast(n,1,MPI_INTEGER,iroot,MPI_COMM_WORLD,ierr)
      call MPI_Bcast(a,n*n,MPI_DOUBLE_PRECISION,iroot,
     *     MPI_COMM_WORLD,ierr)
      call MPI_Bcast(c,n,MPI_DOUBLE_PRECISION,iroot,
     *     MPI_COMM_WORLD,ierr)
      nn=n/nproc
      istart=1+my_rank*nn
      iend=(my_rank+1)*nn
      write(*,*) 'my_rank=',my_rank,'  istart=',istart,' iend=',iend
      do i=1,n
         x(i)=1.
      enddo
      diff=10.
      icount=1
      do while(diff.gt.1.e-8)
         if(my_rank.eq.0) then
            do i=1,n
               write(*,*) 'i=',i,'  x(i)=',x(i),'  y(i)=',y(i)
            enddo
         endif
         do i=1,n
            y(i)=x(i)
         enddo
         do i=istart,iend
            tmp1=0.0
            do j=1,n
               if(j.ne.i) then
                  tmp1=tmp1+a(i,j)*y(j)
               endif
            enddo
            x(i)=(c(i)-tmp1)/a(i,i)
            write(*,*) 'my_rank=',my_rank,'  x(i)=',x(i),'c=',c(i)
         enddo
         do i=1,nproc
            iroot=i-1
            call MPI_Bcast(x(i),1,MPI_DOUBLE_PRECISION,iroot,
     *           MPI_COMM_WORLD,ierr)
         enddo
         diff=0.
         do i=1,n
            diff=diff+dabs(y(i)-x(i))
         enddo
         icount=icount+1
         if(my_rank.eq.0) then
            write(*,*) 'iter=',icount,'  diff=',diff
         endif
      enddo
      if(my_rank.eq.0) then
         write(*,*) 'x(1)=',x(1),'  x(2)=',x(2),'  x(3)=',x(3),
     *        '  x(4)=',x(4)
         do i=1,n
            tmp(i)=0.0
            do j=1,n
               tmp(i)=tmp(i)+a(i,j)*x(j)
            enddo
         enddo
         write(*,*) tmp(1),tmp(2),tmp(3),tmp(4)
      endif
      call MPI_Finalize(ierr)
      stop
      end
c
c      
      subroutine init(nproc,my_rank)
      include "mpif.h"
c
c     initialize the MPI
c
      integer my_rank,nproc
      character*30 nodename
      integer nchar
      call MPI_init(ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,nproc,ierr)
      call MPI_Get_processor_name(nodename,nchar,ierr)
      write(*,*) 'my_rank=',my_rank,' processor=',nodename
      return
      end