C====================== mpi_gather01.f ================================
C
C Calculates the matrix-vector product of a 4x4 matrix A and
C 4x1 vector x, using 4 processes. Originally, each process
C stores a row of A and a single entry of x.
C This method uses a set of 4 gather operations to place a full
C copy of x in each process, which then performs the row-column
C inner-product operation on its data.
C====================================================================

	PROGRAM mpi_gather01
	INCLUDE "mpif.h"

C	MyPE: My logical process number
C	NumPE: Number of processes
C	i: loop index variable
	INTEGER MyPE, NumPE, i, error_code

C	MyA: Single 1x4 row of 4x4 array A
C	MyX: Row entry of 4x1 column vector x
C	FullX: Full copy of x
C	ip: Inner product (result for row i)
C	product: Final vector result
	DOUBLE PRECISION MyA(4), MyX, FullX(4), ip, product(4)

	CALL MPI_INIT(error_code)
	CALL MPI_COMM_RANK(MPI_COMM_WORLD, MyPE, error_code)
	CALL MPI_COMM_SIZE(MPI_COMM_WORLD, NumPE, error_code)

C	Assign myself a row of A and x
	DO i=1, 4
		MyA(i) = MyPE + 1 + 4*(i-1)
	END DO

	MyX = MyPE + 17

C	We want every process to have a copy of x, so
C	we perform gather operations with each process as root
	DO i=1,4
		CALL MPI_GATHER(MyX, 1, MPI_DOUBLE_PRECISION, FullX, 1,
		1	MPI_DOUBLE_PRECISION, i-1, MPI_COMM_WORLD, error_code)
	END DO

C	Now, each process creates inner product for ith row
	ip = 0.0
	DO i=1, 4
		ip = ip + MyA(i)*FullX(i)
	END DO

C	Finally, we gather all process' ip into master process
	CALL MPI_GATHER(ip, 1, MPI_DOUBLE_PRECISION, product, 1,
	1	MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, error_code)
	IF (MyPE == 0) THEN
		WRITE(6,*) "Matrix-vector product is: "
		DO i=1, 4
			WRITE(6,*) product(i)
		END DO
	END IF
	CALL MPI_FINALIZE(error_code)
	END PROGRAM mpi_gather01