	PROGRAM main

	include "mpif.h"

	include 'matrix.inc'

C	Total data area for each processor 
C	(p_imin, p_jmin) ====> left most corner point
C	(p_imax, p_jmax) ====> right most corner point
	INTEGER p_imin(0:MAX_PROCS), p_imax(0:MAX_PROCS)
	INTEGER p_jmin(0:MAX_PROCS), p_jmax(0:MAX_PROCS)

	COMMON	/matrix/ p_imin, p_jmin, p_imax, p_jmax

C	Data area for each time processing in each processor
C	(m_start, n_start) ====> (m_end, n_end)
	INTEGER m_start, m_end
	INTEGER n_start, n_end


C	Maximun memory in each processor
	REAL Y(0:MAX_MEM)
	COMMON /memory/ Y


C	Matrix size (M_total * N_total) and Processor number P
	INTEGER M_total, N_total, P

C	Memory size of each processor
	INTEGER LONG

C	Split processor number into pm * pn
	INTEGER pm, pn

C	MM * NN is the total data amount each processor will process
C	m * n is the data amount each processor will processor each time
C	m, MM is number of rows, n, NN is number of columns
	INTEGER m, n, MM, NN

C	Total computing times will be mdiv * ndiv
C	MM will be split into mdiv parts, and NN will be split into ndiv parts
	INTEGER mdiv, ndiv

C	INTEGER	count
	INTEGER i, j, k

C	me is the processor ID number
	INTEGER me

C	comm(0:4) is the communication buffer
	INTEGER comm(0:4)
	
C	(imin(k), jmax(k)) ====> (imax(k), jmax(k))
C	is the data area each time process on processor k
	INTEGER imin(0:MAX_PROCS), imax(0:MAX_PROCS)
	INTEGER jmin(0:MAX_PROCS), jmax(0:MAX_PROCS)

C	MPI communication error buffer
	INTEGER ierr

C	Initialize
	mdiv = 1
	ndiv = 1
C	count = 0

C	Initilize MPI
	CALL MPI_INIT(ierr)

C	Get the rank of processor
	CALL MPI_COMM_RANK(MPI_COMM_WORLD, me, ierr)

C	in processor 0
	IF (me .eq. 0) THEN

C	Input the matrix size, processor number, and memory size 
	  print *, 'Input the Matrix size (M_rows x N_columns) : '
	  read (5, *) M_total, N_total
	  print *, 'How many processors: '
	  read (5, *) P
	  print *, 'Maximum memory of each processor : '
	  read (5, *) LONG
	  print *, M_total, N_total, P, LONG

C	initializing the matrix
	CALL initialize()

C	Mesh the processors
	CALL process_mesh(P, pm, pn)

C	Caculate the maximum size of matrix each processor will process (MM*NN)
	IF (mod(M_total, pm) .eq. 0) THEN
	  MM = M_total / pm
	ELSE
	  MM = M_total / pm + 1
	END IF

	IF (mod(N_total, pn) .eq. 0) THEN
	  NN = N_total / pn
	ELSE
	  NN = N_total / pn + 1
	END IF

C	partition the matrix
	CALL matrix_part(M_total, N_total, MM, NN, pm, pn)

C 	print out the mesh
	CALL print_mesh(pm, pn)

C	partition the computing area to fit into each processor
	CALL partition(MM, NN, LONG, m, n, mdiv, ndiv)
	print *, 'm = ', m, ' n = ', n, ' mdiv = ', mdiv, 
     &		' ndiv = ', ndiv

C	Broadcast the computing times to each processor
	comm(0) = mdiv
	comm(1) = ndiv
	CALL MPI_BCAST(comm, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)

C	open the result file
	open (UNIT=10, FILE = 'matrix.f.rst')

C	write the matrix size to file
	write(10, *) M_total, N_total

C	Computing repeatly
	DO i = 0, mdiv-1
	  DO j = 0, ndiv-1
	    DO k = 0, P-1

C	Caculate the computing area for each processor
		imin(k) = p_imin(k) + i * m
		IF (imin(k)+m-1 .le. p_imax(k)) THEN
		  imax(k) = imin(k) + m - 1
		ELSE
		  imax(k) = p_imax(k)
		END IF

		jmin(k) = p_jmin(k) + j * n
		IF (jmin(k)+n-1 .le. p_jmax(k)) THEN
		  jmax(k) = jmin(k) + n - 1
		ELSE
		  jmax(k) = p_jmax(k)
		END IF
		
C	Send the computing area and data to each processor
		IF (k .lt. P-1) THEN
		  CALL send_data(k, imin(k), jmin(k), imax(k), jmax(k))

C	remain one part in node 0 computing
		ELSE
		  m_start = imin(k)
		  n_start = jmin(k)
		  m_end = imax(k)
		  n_end = jmax(k)
		END IF
	    END DO
	    
C	Computing on node 0
	    CALL compute(me, m_start, n_start, m_end, n_end)

C	write result in node 0 to file
	    CALL write_to_file(m_start, n_start, m_end, n_end)

C	receive result from other processors and write to file
	    DO k = 0, P-2
		IF (imin(k) .le. p_imax(k) .and. jmin(k) .le. p_jmax(k))
     &			THEN
		  CALL receive_result(k, imin(k), jmin(k), imax(k), jmax(k))
		  CALL write_to_file(imin(k), jmin(k), imax(k), jmax(k))
		END IF
	    END DO
	  END DO
	END DO
	
C	in other processors
	ELSE

C	receive the computing times from node 0
	CALL MPI_BCAST(comm, 2, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
	mdiv = comm(0)
	ndiv = comm(1)
C	print *, 'processor ', me, 'mdiv = ', mdiv, 'ndiv = ', ndiv

	DO i = 0, mdiv-1
	  DO j = 0, ndiv-1

C	receive computing area and data from node 0
	    CALL receive_data(m_start, n_start, m_end, n_end)

C	Computing
	    CALL compute(me, m_start, n_start, m_end, n_end)

C	send the result to node 0
	    CALL send_result(m_start, n_start, m_end, n_end)
	  END DO
	END DO

	END IF

C	MPI finalize
	CALL MPI_FINALIZE(ierr)

	END PROGRAM



C	initialize the matrix data
	SUBROUTINE initialize()
	END SUBROUTINE



C	partition the processor number P to m * n 
	SUBROUTINE process_mesh(P, m, n)
	INTEGER P
	INTEGER m, n

	INTEGER i

	i = 1
	m = 1
	n = P
300	IF (mod(P, i) .eq. 0) THEN
	  m = i
	  n = P / i
	END IF
	i = i + 1
	IF (i .le. n) GOTO 300
	END SUBROUTINE



C	partition the matrix to each processor
C	each processor will process the data between
C	(p_imin(k), p_jmin(k))  =====>  (p_imax(k), p_jmax(k))
	SUBROUTINE matrix_part(M, N, MM, NN, pm, pn)
	include 'matrix.inc'

	INTEGER M, N
	INTEGER MM, NN
	INTEGER pm, pn
	COMMON /matrix/ p_imin, p_jmin, p_imax, p_jmax
	INTEGER p_imin(0:MAX_PROCS), p_imax(0:MAX_PROCS)
	INTEGER p_jmin(0:MAX_PROCS), p_jmax(0:MAX_PROCS)

	INTEGER i, j
	INTEGER procs

	procs = 0
	DO i = 0, pm-1
	  DO j = 0, pn-1
	    p_imin(procs) = i * MM
	    IF (i .lt. pm-1) THEN
	      p_imax(procs) = p_imin(procs) + MM -1
	    ELSE
	      p_imax(procs) = M - 1
	    END IF

	    p_jmin(procs) = j * NN
	    IF (j .lt. pn-1) THEN
	      p_jmax(procs) = p_jmin(procs) + NN -1
	    ELSE
	      p_jmax(procs) = N - 1
	    END IF

	    procs = procs +1
	  END DO
	END DO

	END SUBROUTINE



C	print out the partitioned mesh of data on each processor
	SUBROUTINE print_mesh(pm, pn)
	include 'matrix.inc'

	INTEGER pm, pn
	COMMON /matrix/ p_imin, p_jmin, p_imax, p_jmax
	INTEGER p_imin(0:MAX_PROCS), p_imax(0:MAX_PROCS)
	INTEGER p_jmin(0:MAX_PROCS), p_jmax(0:MAX_PROCS)

	INTEGER procs

	procs = 0

	print *, 'pm = ', pm, ' pn = ', pn
	DO i = 0, pm-1
	  DO j = 0, pn-1
	    print *, '(', p_imin(procs), ',', p_jmin(procs), ') ---> (',
     &		p_imax(procs), ',', p_jmax(procs), ')'
	    procs = procs +1
	  END DO
	END DO

	END SUBROUTINE


	
C	partition the processor area of each processor to fit in
C	the maximum memory ( m * n < LONG )
C	compute times = mdiv * ndiv
	SUBROUTINE partition(MM, NN, LONG, m, n, mdiv, ndiv)
	INTEGER MM, NN
	INTEGER LONG
	INTEGER m, n
	INTEGER mdiv, ndiv

	m = MM
	n = NN
400	IF ((m*n) .le. LONG) GOTO 410

C	split the larger one of column or row
	  IF (n .gt. m) THEN
	    ndiv = ndiv + 1
	    IF (mod(NN, ndiv) .eq. 0) THEN
	      n = NN / ndiv
	    ELSE
	      n = NN / ndiv + 1
	    END IF
	  ELSE
	    mdiv = mdiv + 1
	    IF (mod(MM, mdiv) .eq. 0) THEN
	      m = MM / mdiv
	    ELSE
	      m = MM / mdiv + 1
	    END IF
	  END IF
	GOTO 400
410	END SUBROUTINE



C	send the processing data and area index from node 0 to 
C	other processors
	SUBROUTINE send_data(k, i_start, j_start, i_end, j_end)
	include "mpif.h"
	include 'matrix.inc'
	INTEGER k
	INTEGER i_start, j_start
	INTEGER i_end, j_end

	INTEGER comm(0:4)
	INTEGER TAG
	INTEGER ierr

	TAG = 99

	comm(0) = i_start
	comm(1) = j_start
	comm(2) = i_end
	comm(3) = j_end
	CALL MPI_SEND(comm, 4, MPI_INTEGER, k+1, TAG, 
     &		MPI_COMM_WORLD, ierr)

	END SUBROUTINE



C	processor receive processing data and area index from node 0
	SUBROUTINE receive_data(m_start, n_start, m_end, n_end)
	include "mpif.h"
	include 'matrix.inc'

	INTEGER m_start, n_start
	INTEGER m_end, n_end

	INTEGER comm(0:4)
        INTEGER status(MPI_STATUS_SIZE)
	INTEGER TAG
	INTEGER ierr

	TAG = 99

	CALL MPI_RECV(comm, 4, MPI_INTEGER, 0, TAG, MPI_COMM_WORLD, 
     &		status, ierr)
	m_start = comm(0)
	n_start = comm(1)
	m_end = comm(2)
	n_end = comm(3)

	END SUBROUTINE



C	each processor computing its own data
	SUBROUTINE compute(me, m_start, n_start, m_end, n_end)
	include 'matrix.inc'

	INTEGER me
	INTEGER m_start, n_start
	INTEGER m_end, n_end

	COMMON /memory/ Y
	REAL Y(0:MAX_MEM)

	INTEGER addr

	addr = 0

C	print *, 'processor', me, '(', m_start, ',', n_start, ') ====> (',
C     &		m_end, ',', n_end, ')'

	DO i = m_start, m_end
	  DO j = n_start, n_end
	    Y(addr) = i + j
	    addr = addr + 1
	  END DO
	END DO

	END SUBROUTINE



C	write the results to file
	SUBROUTINE write_to_file(m_start, n_start, m_end, n_end)
	include 'matrix.inc'

	INTEGER m_start, n_start
	INTEGER m_end, n_end

	COMMON /memory/ Y
	REAL Y(0:MAX_MEM)

C	INTEGER addr

	addr = 0

	write(10, *) m_start, n_start, m_end, n_end
	write(10, *) (Y(i), i=0,(m_end-m_start+1)*(n_end-n_start+1)-1)

C	DO i = m_start, m_end
C	  DO j = n_start, n_end
C	    write(10, *) Y(addr)
C	    addr = addr + 1
C	  END DO
C	END DO

	END SUBROUTINE


	
C	send the result on each processor to node 0
	SUBROUTINE send_result(m_start, n_start, m_end, n_end)
	include "mpif.h"
	include 'matrix.inc'

	INTEGER m_start, n_start
	INTEGER m_end, n_end

	COMMON /memory/ Y
	REAL Y(0:MAX_MEM)
	INTEGER TAG
	INTEGER ierr

	TAG = 99

	IF (m_start .le. m_end .and. n_start .le. n_end) THEN
	  CALL MPI_SEND(Y, (m_end-m_start+1)*(n_end-n_start+1), MPI_REAL,
     &		0, TAG, MPI_COMM_WORLD, ierr)
	END IF

	END SUBROUTINE



C	node 0 receive result from other processors
	SUBROUTINE receive_result(k, i_start, j_start, i_end, j_end)
	include "mpif.h"
	include 'matrix.inc'

	INTEGER k
	INTEGER i_start, j_start
	INTEGER i_end, j_end

	COMMON /memory/ Y
	REAL Y(0:MAX_MEM)
	INTEGER TAG
	INTEGER	ierr

        INTEGER status(MPI_STATUS_SIZE)
	
	TAG = 99

	CALL MPI_RECV(Y, (i_end-i_start+1)*(j_end-j_start+1), MPI_REAL,
     &		k+1, TAG, MPI_COMM_WORLD, status, ierr)

	END SUBROUTINE
