// mpi example for matrix vector product
//
//   compute the matrix-vector production
//   by cutting the matrix along rows
//
#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

#define P  4
#define N  1024
#define NL N/P

int main(int argc, char *argv[])
{
  	int myid, np, root;
  	int i, j, lr;
  	double t, (*pA)[N], *px, *pyl, *py;

  	MPI_Init(&argc, &argv);
  	MPI_Comm_rank(MPI_COMM_WORLD, &myid);
  	MPI_Comm_size(MPI_COMM_WORLD, &np);

  	if ( np != P & myid == 0) 
	{  printf("Error: P not equal to np!\n"); MPI_Abort(MPI_COMM_WORLD, 1); }
  
  	pA  = malloc(NL*N*sizeof(double));
  	px  = malloc(N*sizeof(double));
  	pyl = malloc(NL*sizeof(double));
  	py  = malloc(N*sizeof(double));
  	if (!pA | !px | !pyl | !py) 
	{  printf("Out of memory!\n"); MPI_Abort(MPI_COMM_WORLD, 2); }

  	t = MPI_Wtime(); // time

	// set initial value 
  	for(j=0; j<N; j++)  px[j] = 1.0;
  
  	lr = myid*NL;
  	for(i=0; i<NL; i++) 
  	{
    	pyl[i] = 0.0;
    	for(j=0; j<N; j++) 
      		pA[i][j] = 1.0/(double)(i+lr+j+1);
  	}

  	for(i=0; i<NL; i++) 
    	for(j=0; j<N; j++)
	  		pyl[i] = pyl[i] + pA[i][j]*px[j];

  	root = 0;
  	MPI_Gather(pyl, NL, MPI_DOUBLE, py, NL, MPI_DOUBLE, root, MPI_COMM_WORLD);
 
  	t = MPI_Wtime() - t; // time

  	if (myid == root)
  		printf("y[0]=%.6f, y[N-1]=%.6f, Time=%.4e\n", py[0], py[N-1], t);

  	free(pA); free(px); free(pyl); free(py);

  	MPI_Finalize();
  	return 0;
}