/* Mesh skeleton:

   containts: graphical routines setup/visualize, master(), 
   registercoord and registerworkers, plus stub worker() process
*/

#include "mesh.h"

int main(int argc, char** argv)
{
  int rank, size, i, j;
  int themaster = 0;
  MPI_Init(&argc,&argv);
  MPI_Comm_rank(MPI_COMM_WORLD,&rank);
  MPI_Comm_size(MPI_COMM_WORLD,&size);
  if (size != (P*P)+1) {printf("wrong np\n"); exit(1);}

  // seperate worker from master process:  (must be done by all procs)
  MPI_Group worldgroup, workersgroup;
  MPI_Comm_group(MPI_COMM_WORLD,&worldgroup);
  MPI_Group_excl(worldgroup,1,&themaster,&workersgroup);
  MPI_Comm_create(MPI_COMM_WORLD,workersgroup,&WORKERS);

  /***  Define MPI_Datatype for sub-matrix and sub-matrix column ***/
  MPI_Type_vector(M,M,N+2,MPI_DOUBLE,&QUAD);
  MPI_Type_vector(M,1,M,MPI_DOUBLE,&COLUMN);
  MPI_Type_commit(&QUAD);
  MPI_Type_commit(&COLUMN);


  if (rank==0)
      master();
  else  
      worker();

  MPI_Finalize();
}


void registerworkers()
{
  int i, y, x;
  int buf[2];  
  MPI_Status stat;
    // needs to know how to contact each proc in mesh
  for(i=0;i<P*P;i++)
    {
      MPI_Recv(buf,2,MPI_INT,MPI_ANY_SOURCE,1,MPI_COMM_WORLD,&stat);
      y = buf[0];  x = buf[1];
      locate[y][x] = stat.MPI_SOURCE;
      // printf("%d,%d has rank %d\n",y,x,locate[y][x]);  fflush(stdout);
    }
} // registerworkers

void registercoord(int y, int x) // args are your coords in mesh
{  int buf[2];
   buf[0] = y;  buf[1] = x;
   MPI_Send(buf,2,MPI_INT,0,1,MPI_COMM_WORLD);  // tag 1 used
}


/* 3D visualization routines:  - active only when GRAPHICS is 1*/
void setup()
{ 
  int M2;
  M2 = N+2;
  tiltR = toRadians(rotation);   // 140 for rear view, 220 for front
  ZR = 0.5; 
  XDIM=M2*3;  ZDIM=M2;  YDIM = 2*M2;  // x,z will be x,y in simulation
  xshift = (int)(ZDIM * cos(tiltR));
  yshift = (int)(ZDIM * sin(tiltR));

  MPE_Open_graphics(&graph,MPI_COMM_WORLD,(char*)0,0,0,
		    XDIM+(int)abs(xshift)+10,YDIM+(int)abs(yshift)+10, 
		    MPE_GRAPH_INDEPENDENT);
} //setup

void visualize(NROW *BB)
{
  int i, j, k, kr, kd, mrank, U; 
  double z, zr, zd;
  MPE_Color mc;
  // erase background
  MPE_Fill_rectangle(graph,0,0,XDIM+(int)abs(xshift)+10,
		    YDIM+(int)abs(yshift)+10,MPE_WHITE);
  drawbox(graph,0,0,0,XDIM-1,YDIM-1,ZDIM-1,MPE_GREEN);
  U = 4;  // gap between lines in grid
  // draw as a grid
  for(i=0;i<N+2;i+=U)   // rows  - z coord on screen
    for(j=0;j<N+2;j+=U) // cols  - x coord on screen
      {
	mrank = locate[j/(M-1)][i/(M-1)];
        mc = mcolor[mrank%9];
	z = -1*BB[i][j]*(YDIM/2.0);
	k = (int)(z + (YDIM/2.0)+0.5);
	if (k<0) k=0; if (k>=YDIM) k=YDIM-1;
        if (j+U<N+2) 
	  {
	    zr=-1*BB[i][j+U]*(YDIM/2.0);
	    kr = (int)(zr + (YDIM/2.0)+0.5);
	    if (kr<0) kr=0; if (kr>=YDIM) kr=YDIM-1;
	    drawline(graph,3*j,k,i,3*(j+U),kr,i,mc);
	  }
	if (i+U<N+2)
	  {
	    zd=-1*BB[i+U][j]*(YDIM/2.0);
	    kd = (int)(zd + (YDIM/2.0)+0.5);
	    if (kd<0) kd=0; if (kd>=YDIM) kd=YDIM-1;
	    drawline(graph,3*j,k,i,3*j,kd,i+U,mc);
	  }
      } // for for
  MPE_Update(graph);
  if (mousesynch) MPE_Get_mouse_press( graph, &i, &j, &k );
} //visualize



void master()
{
  int i, j,x, y, r, c;
  int resume = 1;
  int vote = 0;  // 0 means should not resume
  static MPI_Status stat[P*P];
  static MPI_Request req[P][P];

  registerworkers(); // sets locate[][] matrix

  // create Data matrix
  Data = (NROW *)malloc(sizeof(double)*(N+2)*(N+2));
  
 // GENERATE Data matrix
  // interior
  for(i=1;i<N+1;i++) // rows (y)
    for(j=1;j<N+1;j++) // cols (x)
    { 
      Data[i][j] = sin(PI/(1+(random() % 4)));
    }
  // set boundary values:
  for(j=0;j<N+2;j++) 
    { 
      Data[0][j] = sin((PI*j)/(N+1));
      Data[N+1][j] = sin((PI*j)/(N+1))*exp(-1*PI);
    }
  for(i=0;i<N+2;i++)
    { Data[i][0] = 0;  Data[i][N+1] = 0; }

  // setup graphics
  if (GRAPHICS)  
    { setup();   visualize(Data); }
  Q = (QROW *)malloc(sizeof(double)*M*M);  // for receiving from worker

  // send quadrants to workers
  for(i=0;i<P;i++) for(j=0;j<P;j++)
    MPI_Send(QUADRANT(Data,i,j),1,QUAD,locate[i][j],0,MPI_COMM_WORLD);
  ///////////////////////

  while (resume)
    {
        MPI_Allreduce(&vote,&resume,1,MPI_INT,MPI_LOR,MPI_COMM_WORLD);
	//printf("master:vote is %d\n",resume);
	if (GRAPHICS)  // if using graphics, need to receive from workers
	  {
	    for(i=0;i<P;i++) for(j=0;j<P;j++) 
	      { 

		//MPI_Recv(QUADRANT(Data,i,j),1,QUAD,locate[i][j],3,MPI_COMM_WORLD,stat+3);

    	        MPI_Recv(Q,M*M,MPI_DOUBLE,
			 locate[i][j],3,MPI_COMM_WORLD,stat); 
		// copy back to Data matrix, without affecting boundary 
		for(r=0;r<M-2;r++) 
		  for(c=0;c<M-2;c++)
		    { Data[r+(i*(N/2)+1)][c+(j*(N/2)+1)] = Q[r+1][c+1]; }
	      } // for i,j

	    visualize(Data);  
	  } // if show graphics
    } // while resume

  printf("Convergence!\n");

  // Receive final portions from worker
  if (!GRAPHICS)
  {
    for(i=0;i<P;i++) for(j=0;j<P;j++) 
      { 
	MPI_Recv(Q,M*M,MPI_DOUBLE,locate[i][j],3,MPI_COMM_WORLD,stat); 
	for(r=0;r<M-2;r++) 
	  for(c=0;c<M-2;c++)
	    { Data[r+(i*(N/2)+1)][c+(j*(N/2)+1)] = Q[r+1][c+1]; }
      } // for i,j
  }
  ///////////////////////
} //master



void worker()
{
  /* setup, create mesh ... */

  registercoord(YCORD,XCORD);

  // create work matrices:
  Q = (QROW *)malloc(sizeof(double)*M*M);
  R = (QROW *)malloc(sizeof(double)*M*M);      

  // receive matrix f
  MPI_Recv(Q,M*M,MPI_DOUBLE,0,0,MPI_COMM_WORLD,stat);

  //  ... 

} //worker