// Program 12.5 Release 1.3 on 30.08.2006
// The Parallel Gauss-Seidel Method (wavefront scheme)
  bool IsPrintEnabled;
  double *u;

int main(int argc, char *argv[])
{
  if (argc < 2)
  {
    printf("USAGE: GZ_ParW.exe N\n");
    return 1;
  }
  if (argc == 2)
    IsPrintEnabled = 0;
  else
    IsPrintEnabled = atoi(argv[2]);

  int N = atoi(argv[1]);

  GZ_Par(N);  

  return 0;
}

int GZ_Par(int N)
{
  int i, j, nx;
  clock_t begin,end;
  
  double temp;
  double dmax;
  double * dm; 
  double d; 
  int Step;
  omp_lock_t dmax_lock;
  omp_init_lock(&dmax_lock);
  
  printf("The Gauss - Seidel method(wavefront scheme). OpenMP Parallel version.\r\n");
  
  if ((u = CreateMatrix(N)) == NULL)
    return 1;

  dm = new double [N];
  
  Step = 0;
  
  printf("Matrix size: %d\r\n", N);

  begin = clock();
  
  do
  {
    dmax = 0;
    for (nx = 1; nx < N - 1; nx++)
    {
      dm[nx] = 0;
#pragma omp parallel for shared(u,nx,dm) private(i,j,temp,d)
      for (i = 1; i < nx + 1; i++ )
      {
        j = nx + 1 - i;
        temp = u[N * i + j];
        u[N * i + j] = 0.25 * (u[N * i + j + 1] + u[N * i + j - 1] + 
          u[N * (i + 1) + j] + u[N * (i - 1) + j]);
        d = fabs(u[N * i + j] - temp);
        if ( dm[i] < d ) 
          dm[i] = d;
      } 
    }
    
    for (nx = N - 2; nx > 1; nx--)
    {
#pragma omp parallel for shared(u,nx,dm) private(i,j,temp,d)
      for (j = N - 2; j > N - nx - 1; j--)
      {
        i = 2 * N - 2 - nx - j;
        temp = u[N * i + j];
        u[N * i + j] = 0.25 * (u[N * i + j + 1] + u[N * i + j - 1] + 
          u[N * (i + 1) + j] + u[N * (i - 1) + j]);
        d = fabs(u[N * i + j] - temp);
        if ( dm[i] < d ) 
          dm[i] = d;
      } 
    }
    
#pragma omp parallel for shared(N, dm, dmax) private(i)
    for (i = 1; i < N - 1; i++ ) 
    {
      omp_set_lock(&dmax_lock);
      if ( dmax < dm[i] ) 
        dmax = dm[i];
      omp_unset_lock(&dmax_lock);
    } 
    
    Step++;
  }
  while (dmax > EPS);
  
  end = clock();
  
  long tck = end - begin;
  
  printf("Results:\r\n- current eps = %f\r\n- required eps = %f\r\n- matrix Size = %d\r\n- iterations = %d\r\n- time = %f\r\n",
    dmax, EPS, N, Step, ((double)tck)/CLOCKS_PER_SEC);
  
  if (IsPrintEnabled)
    PrintMatrix(u, N);
  
  DeleteMatrix(u);
  
  return 0;
}