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

program mat_mat

      implicit none
      include 'mpif.h'

      integer 				:: i
      integer, parameter		:: col = 3
      integer 				:: ierr,rank,nproc,root
      real, dimension(1:col,1:col)	:: A,B,C,A_exact
      real, dimension(1:col,1) 		:: Apart,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 column size.

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


      if (rank == root) then

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

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

      endif   

      ! Put up a barrier until I/O is complete

      call MPI_BARRIER(mpi_comm_world,ierr)

      ! Broadcast matrix B.

      call MPI_BCAST(B,9,mpi_real,root,mpi_comm_world,ierr)

      ! Scatter matrix C by columns.

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

      ! Do the matrix multiplication.

      do i = 1,3
      
         Apart(i,1) = B(i,1)*Cpart(1,1) + B(i,2)*Cpart(2,1) + &
			B(i,3)*Cpart(3,1)
      
      enddo

       ! Gather matrix A.

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

	if (rank == 0) then

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

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


	endif

      call MPI_FINALIZE(ierr)
     

end program mat_mat
