!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!								!
! mat_vec.f							!
!								!
! Created by John Grondalski 6/21/99				!
!								!
! This program computes, in parallel, the matrix-vector product	!
!								!
!	A = B*C							!
!								!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

program mat_vec

      implicit none
      include 'mpif.h'

      integer 			:: i
      integer, parameter	:: col = 4
      integer 			:: ierr,rank,nproc,root
      real, dimension(1:4,1) 	:: A,Apart,Bpart,C,A_exact
      real, dimension(1:4,1:4)	:: B
      real, dimension(1,1)	:: Cpart
      
      root = 0

      ! Initiate mpi.

      call MPI_INIT(ierr)
      call MPI_COMM_RANK(mpi_comm_world, rank, ierr)
      call MPI_COMM_SIZE(mpi_comm_world, nproc,ierr)

! Check that the number of processes started equals the number assumed.

      if (nproc - col /= 0) print*, "Error in number Processors."
      if (nproc - col /= 0) print*, 4, "processes declared in program."
      if (nproc - col /= 0) print*, nproc, "processes allocated."
      if (nproc - col /= 0) Stop



      ! Input B and C.

      if (rank == root) then

	B(1,1) = 1;	B(1,2) = 2;	B(1,3) = 3;	B(1,4) = 4;
	B(2,1) = 4;     B(2,2) = -5;     B(2,3) = 6;	B(2,4) = 4;
	B(3,1) = 7;     B(3,2) = 8;     B(3,3) = 9;	B(3,4) = 2;
        B(4,1) = 3;     B(4,2) = -1;     B(4,3) = 5;     B(4,4) = 0;

        C(1,1) = 1;
        C(2,1) = -4;
        C(3,1) = 7;
	C(4,1) = 3;

      endif   

      ! Put up a barrier until I/O is complete

      call MPI_BARRIER(mpi_comm_world,ierr)

      ! Scatter matrix B by columns.

     call MPI_SCATTER(B,col,mpi_real,Bpart,col,mpi_real, &
	root,mpi_comm_world,ierr)

      ! Scatter matrix C by rows.

       call MPI_SCATTER(C,1,mpi_real,Cpart,1,mpi_real, &
       root,mpi_comm_world,ierr)

      ! Do the vector-scalar multiplication.

        Apart = matmul(Bpart,Cpart)

      ! Reduce to  matrix A.

      call MPI_REDUCE(Apart,A,col,mpi_real,mpi_sum, &
      root,mpi_comm_world,ierr)

	if (rank == 0) then

      print*,''
      print*,'This is the result of the parallel computation:'
      print*,''
      print "(F6.1)", A(1,1)
      print "(F6.1)", A(2,1)
      print "(F6.1)", A(3,1)
      print "(F6.1)", A(4,1)
      
	A_exact = matmul(B,C)

      print*,''
      print*,'This is the result of the Fortran intrinsic subroutine:'
      print*,''
      print "(F6.1)", A_exact(1,1)
      print "(F6.1)", A_exact(2,1)
      print "(F6.1)", A_exact(3,1)
      print "(F6.1)", A_exact(4,1)

	endif

      call MPI_FINALIZE(ierr)
     

end program mat_vec
