#define SQUARE(x) SYM_multiply(x,x)

void printSolution() {
  int i, j;

  update_seqpar();
#ifdef MPI_VERBOSE
  for (i = 0; i < nx; i++) {
    for (j = 0; j < ny; j++) {
      printf("u[%d][%d] = ", i, j);
      SYM_printExpr(now.u_seq[i].data[j]);
      printf("\n");
    }
    printf("\n");
  }
  printf("\n");
#endif
}

void init_seq() {
  int i, j;

  for (i = 0; i < nx; i++)
    for (j = 0; j < ny; j++)
      now.u_seq[i].data[j] =
	SYM_if
	(SYM_lessThan
	 (SYM_add
	  (SQUARE
	   (SYM_subtract
	    (SYM_multiply(dx, SYM_intConstant(i)),
	     SYM_multiply
	     (POINT5,
	      SYM_multiply(dx, SYM_intConstant(nx-1))))),
	   SQUARE
	   (SYM_subtract
	    (SYM_multiply(dx, SYM_intConstant(j)),
	     SYM_multiply
	     (POINT5,
	      SYM_multiply(dx, SYM_intConstant(ny-1)))))),
	  SQUARE(r)),
	 SYM_intConstant(100),
	 SYM_ZERO);
}

void update_seq() {
  int i, j;
  MPI_Symbolic u_new[nx][ny];
  MPI_Symbolic k = SYM_divide(SYM_multiply(D,dt),SQUARE(dx));

  for (i = 0; i < nx; i++)
    for (j = 0; j < ny; j++)
      u_new[i][j] =
	SYM_add
	(now.u_seq[i].data[j],
	 SYM_multiply
	 (k,
	  SYM_subtract
	  (SYM_add
	   (SYM_add
	    (SYM_add(now.u_seq[(i+1)%nx].data[j],
		     now.u_seq[(i+nx-1)%nx].data[j]),
	     now.u_seq[i].data[(j+1)%ny]),
	    now.u_seq[i].data[(j+ny-1)%ny]),
	   SYM_multiply(SYM_intConstant(4), now.u_seq[i].data[j]))));
  for (i = 0; i < nx; i++)
    for (j = 0; j < ny; j++)
      now.u_seq[i].data[j] = u_new[i][j];
}

void main_seq() {
  int iter;

  init_seq();
  for (iter = 1; iter <= nsteps; iter++) {
    update_seq();
  }
}


void init_par(struct ParallelColumn u_par[], uchar xcoord, uchar ycoord) {
  int i, j;

  for (i = 1; i <= nxl; i++)
    for (j = 1; j <= nyl; j++)
      u_par[i].data[j] =
	SYM_if
	(SYM_lessThan
	 (SYM_add
	  (SQUARE
	   (SYM_subtract
	    (SYM_multiply(dx, SYM_intConstant(nxl*xcoord+i-1)),
	     SYM_multiply
	     (POINT5,
	      SYM_multiply(dx, SYM_intConstant(nx-1))))),
	   SQUARE
	   (SYM_subtract
	    (SYM_multiply(dx, SYM_intConstant(nyl*ycoord+j-1)),
	     SYM_multiply
	     (POINT5,
	      SYM_multiply(dx, SYM_intConstant(ny-1)))))),
	  SQUARE(r)),
	 SYM_intConstant(100),
	 SYM_ZERO);
}

void update_par(struct ParallelColumn u_par[]) {
  int i, j;
  MPI_Symbolic u_new[nxl+2][nyl+2];
  MPI_Symbolic k = SYM_divide(SYM_multiply(D,dt),SQUARE(dx));

  for (i = 1; i <= nxl; i++)
    for (j = 1; j <= nyl; j++)
      u_new[i][j] =
	SYM_add
	(u_par[i].data[j],
	 SYM_multiply
	 (k,
	  SYM_subtract
	  (SYM_add
	   (SYM_add
	    (SYM_add(u_par[i+1].data[j],
		     u_par[i-1].data[j]),
	     u_par[i].data[j+1]),
	    u_par[i].data[j-1]),
	   SYM_multiply(SYM_intConstant(4), u_par[i].data[j]))));
  for (i = 1; i <= nxl; i++)
    for (j = 1; j <= nyl; j++)
      u_par[i].data[j] = u_new[i][j];
}


int agree(uchar iter, uchar m, uchar n, uchar j, MPI_Symbolic buf[],
	  struct SequentialColumn u_seq[]) {
  int i;

  if (iter != nsteps) return 1;
  for (i = 1; i <= nxl; i++) {
    int i_g = m*nxl+i-1;
    int j_g = n*nyl+j-1;
    if (buf[i-1] !=
	now.u_seq[i_g].data[j_g]) {
      printf("%s: m=%d, i=%d, n=%d, j=%d, (%d,%d):\n",
	     "Discrepancy", m, i, n,
	     j, i_g, j_g);
      printf("Sequential result: ");
      SYM_printExpr(now.u_seq[i_g].data[j_g]);
      printf("\nParallel result: ");
      SYM_printExpr(buf[i-1]);
      printf("\n");
      return 0;
    }
  }
  return 1;
}
