c	matrix.f
c	Written By Drew Weitz
c	Last Modified: June 9, 2003

c	This program is a demostration of some very basic MPI function calls.
c	For information about each specific function call, please check out
c	the web for information on MPI.  Unfortuately, man pages are only
c	installed for the C language MPI calls.  One good source for Fortran
c	information is: www.iitk.ac.in/cc/param/mpi_calls.html

c	Programmers Note:  I am not a Fortran programmer.  This program was
c	written as an exercise to attempt to learn both Fortran and MPI at the
c	same time.  If the code presented contains awful formatting or coding
c	style, then I apologize.


	program main

c	To include MPI functions
	include 'mpif.h'

c	Rank is the "processor number" assigned each processor
c	at program startup.  Nprocs is the total number of processors
c	involved in the program run.  Each node has a copy of all
c	variables declared.

	integer rc, rank, nprocs

c	MPI_INIT insures that the cluster nodes can talk to one another.
c	It is necessary for MPI code.
	call MPI_INIT(ierr)

c	MPI_COMM_RANK initializes each nodes' copy of the "rank" variable
c	and sets it equal to the process number it was assigned.
c       MPI_COMM_WORLD is a macro that MPI uses to address the correct
c       MPI network.  Since multiple MPI jobs could be running at once,
c       we want to have a way of addressing only the processors in our program
c       run (or our "world").

	call MPI_COMM_RANK(MPI_COMM_WORLD, rank, ierr)

c	MPI_COMM_SIZE sets each nodes' copy of the "nprocs" equal to the
c	number of processors avaiable (the number of processors dedicated to
c	an MPI program run is designated at runtime as an argument for the
c	MPIRUN script.
	call MPI_COMM_SIZE(MPI_COMM_WORLD, nprocs, ierr)

c	Matrix is the function that does the bulk of the work in the program.
c	The idea is that each node is given a different "grid" or 2-D array.
c	Each node the modifies it's own section in some well-defined manner,
c	and then returns the result to node 0.  Node 0 then stacks these
c	together to form a "cube" or 3-D array, and then collapses this
c	cube (or Matrix) into a two dimensional array by taking the sum of
c	all the elements in one direction.
	call matrix(rank, nprocs, MPI_COMM_WORLD, MPI_INTEGER)

c	MPI_FINALIZE closes all the connections opened during the MPI program
c	run.  It is neccessary for MPI Programs.
	call MPI_FINALIZE(rc)

	stop
	end





	subroutine matrix(rank, nprocs, MPI_COMM_WORLD, MPI_INTEGER)

c	variable declarations and initializations
	integer rank, nprocs, i, j, z, maxsize, ierr, MPI_COMM_WORLD
	integer MPI_INTEGER
	integer totalsize
	integer status(4)
	integer thematrix(nprocs-1, nprocs-1)
	integer thecube(nprocs-1, nprocs-1, nprocs-1)

	maxsize = nprocs-1
	if (rank .eq. 0) then
c	Initialize the matrix
		do 10 i=1, maxsize
			do 20 j=1, maxsize
				thematrix(i,j)=0
 20			continue
 10		continue

		write (*,*) "The matrix is: "

		call showmatrix(thematrix, maxsize)






c	Send the matrix out to the other ranks
		write(*,*) "Sending it out"
		do 30 i=1, maxsize

c	MPI_SEND is the function call that sends the values of one nodes'
c	variable out to another node.  In effect, it is the way in which
c	data can be passed to another node.  A node can either receieve this
c	"message" (via MPI_RECV), or it can be ignored.
			call MPI_SEND(thematrix, maxsize*maxsize,
     +				MPI_INTEGER, i, 100, MPI_COMM_WORLD,
     +				ierr)
 30		continue

c	Receive the matricies back from the ranks

		write(*,*) "waiting to receive..."
		do 50 z=1, maxsize
c       MPI_Recv is the compliament to MPI_Send.  This is one way in which
c       the processor will receive data.  In the case of MPI_RECV, the
c       processor will stop executing its code and wait for the expected
c       message.   If you do not want the processor to stop executing code
c       while waiting for a "message", you must use a "non-blocking" 
c       receive function.  See the website listed at the top of the program
c       for more information.

			call MPI_RECV(thematrix, maxsize*maxsize,
     +				MPI_INTEGER, z, 100+z, MPI_COMM_WORLD, 
     +				status, ierr)
			do 60 i=1, maxsize
				do 70 j=1, maxsize
					thecube(z,i,j)=thematrix(i,j)
 70				continue
 60			continue
 50		continue

		write (*,*) "Recombining data..."

		do i=1, maxsize
			do j=1, maxsize
				thematrix(i,j)=0
				do z=1, maxsize
					thematrix(i,j)=thematrix(i,j)+
     +							thecube(z,i,j)
				enddo
			enddo
		enddo

		call showmatrix(thematrix, maxsize)
	elseif (rank .gt. 0) then




c	Receive the matrix from the master

		call MPI_RECV(thematrix, maxsize*maxsize, MPI_INTEGER, 
     +				0, 100, MPI_COMM_WORLD, status, ierr)

c	Update the matrix
		do i=1, maxsize
			thematrix(rank,i)=rank*i
		enddo
c	Send the matrix back

		call MPI_SEND(thematrix, maxsize*maxsize, MPI_INTEGER,
     +				0, 100+rank, MPI_COMM_WORLD, ierr)

	endif
	

	return
	end






c	This subroutine takes a 2-D array and prints it to the screen.

	subroutine showmatrix(thematrix, maxsize)
	integer maxsize
	integer thematrix(maxsize, maxsize)
	integer i,j
	do 30 i=1, maxsize
		write (*, 200) (thematrix(i,j), j=1,maxsize)
 200	format (10000I8)
 30	continue
	return
	end

c	This subroutine takes a 3-D array and prints each "level" to
c	the screen.

	subroutine showcube(thecube, maxsize)
	integer maxsize
	integer thecube(maxsize, maxsize, maxsize)
	integer i,j,z
	do z=1, maxsize
		write (*,*) "Level ", z
		do i=1, maxsize
			write (*,*) (thecube(z,i,j), j=1, maxsize)
		enddo
	enddo
	return
	end