The Higher Education and Research forge

Home My Page Projects Code Snippets Project Openings MMD
Summary Activity Forums Tracker Lists Tasks Docs Surveys News SCM Files

SCM Repository

authorAlain <alain.bastide@univ-reunion.fr>
Fri, 10 Jul 2020 08:44:02 +0000 (12:44 +0400)
committerAlain <alain.bastide@univ-reunion.fr>
Sat, 11 Jul 2020 13:54:11 +0000 (17:54 +0400)
CMakeLists.txt
TODO.md
src/CMakeLists.txt
src/main.c
src/mmd.h
test/CMakeLists.txt
test/test_hybrid.c [new file with mode: 0644]
test/test_mpi.c [new file with mode: 0644]
test/test_openmp.c [new file with mode: 0644]

index 40b019d..c084a03 100644 (file)
@@ -18,6 +18,7 @@ set(CMAKE_C_FLAGS "-Wall -Wextra")
 set(CMAKE_C_FLAGS_DEBUG "-g -DDEBUG -march=skylake     -mtune=skylake")
 set(CMAKE_C_FLAGS_RELEASE "-O3 -DNDEBUG -march=skylake         -mtune=skylake")
 
+
 #set(CMAKE_C_FLAGS_DEBUG "put your flags")
 #set(CMAKE_C_FLAGS_MINSIZEREL "put your flags")
 #set(CMAKE_C_FLAGS_RELWITHDEBINFO "put your flags")
@@ -55,8 +56,8 @@ set(DATA_FILE ${PROJECT_SOURCE_DIR}/README.md  ${PROJECT_SOURCE_DIR}/TODO.md )
 add_subdirectory(${SOURCES})
 add_subdirectory(${TESTS})
   
-enable_testing ()
-add_test (NAME MyTest COMMAND test01.out ${DATA_FILE})
+#enable_testing ()
+#add_test (NAME MyTest COMMAND test01.out ${DATA_FILE})
 
 
 
diff --git a/TODO.md b/TODO.md
index d9fd0c3..8f6f178 100644 (file)
--- a/TODO.md
+++ b/TODO.md
@@ -38,6 +38,9 @@ Reduce the execution time of the functions :
      * https://matteding.github.io/2019/04/25/sparse-matrices/
      * https://slepc.upv.es/
      
+    Scaling perf
+     * https://scalasca.org/scalasca/front_content.php?idart=1070
+     
 
 
 
index 9a062bf..9bccd26 100644 (file)
@@ -5,9 +5,15 @@ add_library (savedata STATIC savedata.c savedata.h mmd.h)
 add_library (memory STATIC memory.c memory.h mmd.h)
 add_library (tensor STATIC tensor.c tensor.h mmd.h)
 add_library (mmd STATIC mmd.c mmd.h)
-add_executable(${PROJECT_NAME} main.c ${DATA_FILE})
-#SET_TARGET_PROPERTIES(${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-DDEBUG")
+set(OTHER_FILES ../README.md ../LICENSE.md ../TODO.md)
 
+add_executable(${PROJECT_NAME} main.c ${DATA_FILE} ${OTHER_FILES})
+set_target_properties( ${PROJECT_NAME}
+    PROPERTIES
+    ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
+)
 target_link_libraries (${PROJECT_NAME}
     primal
     listop
@@ -20,6 +26,33 @@ target_link_libraries (${PROJECT_NAME}
     ${VTK_LIBRARIES}
     )
 
+
+add_executable(parallel_${PROJECT_NAME} main.c ${DATA_FILE})
+set_target_properties( parallel_${PROJECT_NAME}
+    PROPERTIES
+    ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
+)
+target_link_libraries (parallel_${PROJECT_NAME}
+    primal
+    listop
+    dual
+    savedata
+    tensor
+    memory
+    mmd
+    m
+    ${VTK_LIBRARIES}
+    )
+#set(FLAGS_GNU "-fopenmp -DMMD_PARALLEL")
+set(FLAGS_GNU "")
+set(CMAKE_C_FLAGS  "${CMAKE_C_FLAGS} ${FLAGS_GNU}")
+
+#set_target_properties(parallel_${PROJECT_NAME} PROPERTIES COMPILE_FLAGS "-fopenmp  -DPARALLEL" LINKER_FLAG "-fopenmp -lgomp -DPARALLEL")
+
+
+
 #add_library(mmd_link SHARED shared.c)
 #target_compile_definitions(mmd_link PUBLIC mmd_link)
 
index b10bb33..d8ce61f 100644 (file)
@@ -195,7 +195,14 @@ int main (void)
   
   freeMemory(&myPrimalMesh, &myDualMesh);
   todo();
-  
+
+#ifdef MMD_PARALLEL
+#pragma omp parallel
+{
+  printf("Hello from process: %d\n", omp_get_thread_num());  
+
+}
+#endif
 }
 
 
index fc8a7bc..5cef64e 100644 (file)
--- a/src/mmd.h
+++ b/src/mmd.h
@@ -35,6 +35,7 @@ Author
 
 #include <stdio.h>
 #include <stdlib.h>
+#include <stdbool.h>
 #include <string.h>
 
 #include <math.h>
@@ -47,6 +48,10 @@ Author
 
 #include <vtkCellType.h>
 
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 #define NBAFFICHE 10
 
 #define DIM3D 3
index 59ef227..f340625 100644 (file)
@@ -4,13 +4,59 @@ include_directories (${MMD_SOURCE_DIR}/src)
 add_definitions("-DDEBUG")
 set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O0 -ggdb -DDEBUG")
 add_executable(test01.out test01.c ${DATA_FILE})
-#SET_TARGET_PROPERTIES(test01.out PROPERTIES COMPILE_FLAGS "-DDEBUG")
+set_target_properties( test01.out
+    PROPERTIES
+    ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
+    )
 target_link_libraries (test01.out
-                       memory
-                       mmd
-                       primal
-                       listop
-                       m
-                       )
+    memory
+    mmd
+    primal
+    listop
+    m
+    )
+
+
+find_package(MPI)
+SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")                   
+include_directories(SYSTEM ${MPI_INCLUDE_PATH})
+add_executable(test_mpi.out test_mpi.c ${DATA_FILE})
+set_target_properties(test_mpi.out
+    PROPERTIES
+    ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
+    )
+target_link_libraries (test_mpi.out
+    memory
+    mmd
+    primal
+    listop
+    m
+    ${MPI_C_LIBRARIES}
+    )
+
+
+
+find_package(MPI)
+SET(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fopenmp")                   
+include_directories(SYSTEM ${MPI_INCLUDE_PATH})
+add_executable(test_openmp.out test_openmp.c ${DATA_FILE})
+set_target_properties(test_openmp.out
+    PROPERTIES
+    ARCHIVE_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    LIBRARY_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/lib"
+    RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/bin"
+    )
+target_link_libraries (test_openmp.out
+    memory
+    mmd
+    primal
+    listop
+    m
+    ${MPI_C_LIBRARIES}
+    )
+
+
diff --git a/test/test_hybrid.c b/test/test_hybrid.c
new file mode 100644 (file)
index 0000000..51db61d
--- /dev/null
@@ -0,0 +1,27 @@
+#include <stdio.h>
+#include <mpi.h>
+#include <omp.h>
+
+int main(int argc, char *argv[])
+{
+    int numprocs, rank, namelen;
+    char processor_name[MPI_MAX_PROCESSOR_NAME];
+    int iam = 0, np = 1;
+
+    MPI_Init(&argc, &argv);
+    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
+    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Get_processor_name(processor_name, &namelen);
+
+    #pragma omp parallel default(shared) private(iam, np)
+    {
+        np = omp_get_num_threads();
+        iam = omp_get_thread_num();
+        printf("Hybrid: Hello from thread %d out of %d from process %d out of %d on %s \n",
+                iam, np, rank, numprocs, processor_name );
+    }
+
+    MPI_Finalize();
+
+    return 0;
+}
diff --git a/test/test_mpi.c b/test/test_mpi.c
new file mode 100644 (file)
index 0000000..a0a003d
--- /dev/null
@@ -0,0 +1,451 @@
+# include <math.h>
+# include <mpi.h>
+# include <stdio.h>
+# include <stdlib.h>
+# include <string.h>
+# include <time.h>
+
+double L = 1.0;                        /* linear size of square region */
+int N = 32;                    /* number of interior points per dim */
+
+double *u, *u_new;             /* linear arrays to hold solution */
+
+/* macro to index into a 2-D (N+2)x(N+2) array */
+#define INDEX(i,j) ((N+2)*(i)+(j))
+
+int my_rank;                   /* rank of this process */
+
+int *proc;                     /* process indexed by vertex */
+int *i_min, *i_max;            /* min, max vertex indices of processes */
+int *left_proc, *right_proc;   /* processes to left and right */
+
+/*
+  Functions:
+*/
+int main ( int argc, char *argv[] );
+void allocate_arrays ( );
+void jacobi ( int num_procs, double f[] );
+void make_domains ( int num_procs );
+double *make_source ( );
+void timestamp ( );
+
+/******************************************************************************/
+
+int main ( int argc, char *argv[] ) 
+
+{
+  double change;
+  double epsilon = 1.0E-03;
+  double *f;
+  char file_name[100];
+  int i;
+  int j;
+  double my_change;
+  int my_n;
+  int n;
+  int num_procs;
+  int step;
+  double *swap;
+  double wall_time;
+/*
+  MPI initialization.
+*/
+  MPI_Init ( &argc, &argv );
+
+  MPI_Comm_size ( MPI_COMM_WORLD, &num_procs );
+
+  MPI_Comm_rank ( MPI_COMM_WORLD, &my_rank );
+/*
+  Read commandline arguments, if present.
+*/
+  if ( 1 < argc )
+  {
+    sscanf ( argv[1], "%d", &N );
+  }
+  else
+  {
+    N = 32;
+  }
+
+  if ( 2 < argc )
+  {
+    sscanf ( argv[2], "%lf", &epsilon );
+  }
+  else
+  {
+    epsilon = 1.0E-03;
+  }
+  if ( 3 < argc )
+  {
+    strcpy ( file_name, argv[3] );
+  }
+  else
+  {
+    strcpy ( file_name, "poisson_mpi.out" );
+  }
+/*
+  Print out initial information.
+*/
+  if ( my_rank == 0 ) 
+  {
+    timestamp ( );
+    printf ( "\n" );
+    printf ( "POISSON_MPI:\n" );
+    printf ( "  C version\n" );
+    printf ( "  2-D Poisson equation using Jacobi algorithm\n" );
+    printf ( "  ===========================================\n" );
+    printf ( "  MPI version: 1-D domains, non-blocking send/receive\n" );
+    printf ( "  Number of processes         = %d\n", num_procs );
+    printf ( "  Number of interior vertices = %d\n", N );
+    printf ( "  Desired fractional accuracy = %f\n", epsilon );
+    printf ( "\n" );
+  }
+
+  allocate_arrays ( );
+  f = make_source ( );
+  make_domains ( num_procs );
+
+  step = 0;
+/*
+  Begin timing.
+*/
+  wall_time = MPI_Wtime ( );
+/*
+  Begin iteration.
+*/
+  do 
+  {
+    jacobi ( num_procs, f );
+    ++step;
+/* 
+  Estimate the error 
+*/
+    change = 0.0;
+    n = 0;
+
+    my_change = 0.0;
+    my_n = 0;
+
+    for ( i = i_min[my_rank]; i <= i_max[my_rank]; i++ )
+    {
+      for ( j = 1; j <= N; j++ )
+      {
+        if ( u_new[INDEX(i,j)] != 0.0 ) 
+        {
+          my_change = my_change 
+            + fabs ( 1.0 - u[INDEX(i,j)] / u_new[INDEX(i,j)] );
+
+          my_n = my_n + 1;
+        }
+      }
+    }
+    MPI_Allreduce ( &my_change, &change, 1, MPI_DOUBLE, MPI_SUM,
+      MPI_COMM_WORLD );
+
+    MPI_Allreduce ( &my_n, &n, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
+
+    if ( n != 0 )
+    {
+      change = change / n;
+    }
+    if ( my_rank == 0 && ( step % 10 ) == 0 ) 
+    {
+      printf ( "  N = %d, n = %d, my_n = %d, Step %4d  Error = %g\n", 
+        N, n, my_n, step, change );
+    }
+/* 
+  Interchange U and U_NEW.
+*/
+    swap = u;
+    u = u_new;
+    u_new = swap;
+  } while ( epsilon < change );
+
+/* 
+  Here is where you can copy the solution to process 0 
+  and print to a file.
+*/
+
+/*
+  Report on wallclock time.
+*/
+  wall_time = MPI_Wtime() - wall_time;
+  if ( my_rank == 0 )
+  {
+    printf ( "\n" );
+    printf ( "  Wall clock time = %f secs\n", wall_time );
+  }
+/*
+  Terminate MPI.
+*/
+  MPI_Finalize ( );
+/*
+  Free memory.
+*/
+  free ( f );
+/*
+  Terminate.
+*/
+  if ( my_rank == 0 )
+  {
+    printf ( "\n" );
+    printf ( "POISSON_MPI:\n" );
+    printf ( "  Normal end of execution.\n" );
+    printf ( "\n" );
+    timestamp ( );
+  }
+  return 0;
+}
+/******************************************************************************/
+
+void allocate_arrays ( ) 
+
+{
+  int i;
+  int ndof;
+
+  ndof = ( N + 2 ) * ( N + 2 );
+
+  u = ( double * ) malloc ( ndof * sizeof ( double ) );
+  for ( i = 0; i < ndof; i++)
+  {
+    u[i] = 0.0;
+  }
+
+  u_new = ( double * ) malloc ( ndof * sizeof ( double ) );
+  for ( i = 0; i < ndof; i++ )
+  {
+    u_new[i] = 0.0;
+  }
+
+  return;
+}
+/******************************************************************************/
+
+void jacobi ( int num_procs, double f[] ) 
+
+{
+  double h;
+  int i;
+  int j;
+  MPI_Request request[4];
+  int requests;
+  MPI_Status status[4];
+/*
+  H is the lattice spacing.
+*/
+  h = L / ( double ) ( N + 1 );
+/* 
+  Update ghost layers using non-blocking send/receive 
+*/
+  requests = 0;
+
+  if ( left_proc[my_rank] >= 0 && left_proc[my_rank] < num_procs ) 
+  {
+    MPI_Irecv ( u + INDEX(i_min[my_rank] - 1, 1), N, MPI_DOUBLE,
+      left_proc[my_rank], 0, MPI_COMM_WORLD,
+      request + requests++ );
+
+    MPI_Isend ( u + INDEX(i_min[my_rank], 1), N, MPI_DOUBLE,
+      left_proc[my_rank], 1, MPI_COMM_WORLD,
+      request + requests++ );
+  }
+
+  if ( right_proc[my_rank] >= 0 && right_proc[my_rank] < num_procs ) 
+  {
+    MPI_Irecv ( u + INDEX(i_max[my_rank] + 1, 1), N, MPI_DOUBLE,
+      right_proc[my_rank], 1, MPI_COMM_WORLD,
+      request + requests++ );
+
+    MPI_Isend ( u + INDEX(i_max[my_rank], 1), N, MPI_DOUBLE,
+      right_proc[my_rank], 0, MPI_COMM_WORLD,
+      request + requests++ );
+  }
+/* 
+  Jacobi update for internal vertices in my domain.
+*/
+  for ( i = i_min[my_rank] + 1; i <= i_max[my_rank] - 1; i++ )
+  {
+    for ( j = 1; j <= N; j++ )
+    {
+      u_new[INDEX(i,j)] =
+        0.25 * ( u[INDEX(i-1,j)] + u[INDEX(i+1,j)] +
+                 u[INDEX(i,j-1)] + u[INDEX(i,j+1)] +
+                 h * h * f[INDEX(i,j)] );
+    }
+  }
+/* 
+  Wait for all non-blocking communications to complete.
+*/
+  MPI_Waitall ( requests, request, status );
+/* 
+  Jacobi update for boundary vertices in my domain.
+*/
+  i = i_min[my_rank];
+  for ( j = 1; j <= N; j++ )
+  {
+    u_new[INDEX(i,j)] =
+      0.25 * ( u[INDEX(i-1,j)] + u[INDEX(i+1,j)] +
+               u[INDEX(i,j-1)] + u[INDEX(i,j+1)] +
+               h * h * f[INDEX(i,j)] );
+  }
+
+  i = i_max[my_rank];
+  if (i != i_min[my_rank])
+  {
+    for (j = 1; j <= N; j++)
+    {
+      u_new[INDEX(i,j)] =
+        0.25 * ( u[INDEX(i-1,j)] + u[INDEX(i+1,j)] +
+                 u[INDEX(i,j-1)] + u[INDEX(i,j+1)] +
+                 h * h * f[INDEX(i,j)] );
+    }
+  }
+
+  return;
+}
+/******************************************************************************/
+
+void make_domains ( int num_procs ) 
+
+
+{
+  double d;
+  double eps;
+  int i;
+  int p;
+  double x_max;
+  double x_min;
+/* 
+  Allocate arrays for process information.
+*/
+  proc = ( int * ) malloc ( ( N + 2 ) * sizeof ( int ) );
+  i_min = ( int * ) malloc ( num_procs * sizeof ( int ) );
+  i_max = ( int * ) malloc ( num_procs * sizeof ( int ) );
+  left_proc = ( int * ) malloc ( num_procs * sizeof ( int ) );
+  right_proc = ( int * ) malloc ( num_procs * sizeof ( int ) );
+/* 
+  Divide the range [(1-eps)..(N+eps)] evenly among the processes.
+*/
+  eps = 0.0001;
+  d = ( N - 1.0 + 2.0 * eps ) / ( double ) num_procs;
+
+  for ( p = 0; p < num_procs; p++ )
+  {
+/* 
+  The I indices assigned to domain P will satisfy X_MIN <= I <= X_MAX.
+*/
+    x_min = - eps + 1.0 + ( double ) ( p * d );
+    x_max = x_min + d;
+/* 
+  For the node with index I, store in PROC[I] the process P it belongs to.
+*/
+    for ( i = 1; i <= N; i++ )
+    {
+      if ( x_min <= i && i < x_max )
+      {
+        proc[i] = p;
+      }
+    }
+  }
+/* 
+  Now find the lowest index I associated with each process P.
+*/
+  for ( p = 0; p < num_procs; p++ )
+  {
+    for ( i = 1; i <= N; i++ )
+    {
+      if ( proc[i] == p )
+      {
+        break;
+      }
+    }
+    i_min[p] = i;
+/* 
+  Find the largest index associated with each process P.
+*/
+    for ( i = N; 1 <= i; i-- )
+    {
+      if ( proc[i] == p )
+      {
+        break;
+      }
+    }
+    i_max[p] = i;
+/* 
+  Find the processes to left and right. 
+*/
+    left_proc[p] = -1;
+    right_proc[p] = -1;
+
+    if ( proc[p] != -1 ) 
+    {
+      if ( 1 < i_min[p] && i_min[p] <= N )
+      {
+        left_proc[p] = proc[i_min[p] - 1];
+      }
+      if ( 0 < i_max[p] && i_max[p] < N )
+      {
+        right_proc[p] = proc[i_max[p] + 1];
+      }
+    }
+  }
+
+  return;
+}
+/******************************************************************************/
+
+double *make_source ( ) 
+
+{
+  double *f;
+  int i;
+  int j;
+  int k;
+  double q;
+
+  f = ( double * ) malloc ( ( N + 2 ) * ( N + 2 ) * sizeof ( double ) );
+
+  for ( i = 0; i < ( N + 2 ) * ( N + 2 ); i++ )
+  {
+    f[i] = 0.0;
+  }
+/* 
+  Make a dipole.
+*/
+  q = 10.0;
+
+  i = 1 + N / 4;
+  j = i;
+  k = INDEX ( i, j );
+  f[k] = q;
+
+  i = 1 + 3 * N / 4;
+  j = i;
+  k = INDEX ( i, j );
+  f[k] = -q;
+
+  return f;
+}
+/******************************************************************************/
+
+void timestamp ( )
+
+{
+# define TIME_SIZE 40
+
+  static char time_buffer[TIME_SIZE];
+  const struct tm *tm;
+  time_t now;
+
+  now = time ( NULL );
+  tm = localtime ( &now );
+
+  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
+
+  printf ( "%s\n", time_buffer );
+
+  return;
+# undef TIME_SIZE
+}
diff --git a/test/test_openmp.c b/test/test_openmp.c
new file mode 100644 (file)
index 0000000..97e2844
--- /dev/null
@@ -0,0 +1,615 @@
+# include <stdlib.h>
+# include <stdio.h>
+# include <math.h>
+# include <time.h>
+# include <omp.h>
+
+# define NX 361
+# define NY 361
+
+int main ( int argc, char *argv[] );
+double r8mat_rms ( int nx, int ny, double a[NX][NY] );
+void rhs ( int nx, int ny, double f[NX][NY] );
+void sweep ( int nx, int ny, double dx, double dy, double f[NX][NY], 
+             int itold, int itnew, double u[NX][NY], double unew[NX][NY] );
+void timestamp ( );
+double u_exact ( double x, double y );
+double uxxyy_exact ( double x, double y );
+
+/******************************************************************************/
+
+int main ( int argc, char *argv[] )
+
+/******************************************************************************/
+/*
+  Purpose:
+  
+    MAIN is the main program for POISSON_OPENMP.
+    
+  Discussion:
+  
+    POISSON_OPENMP is a program for solving the Poisson problem.
+    
+    This program uses OpenMP for parallel execution.
+    
+    The Poisson equation
+    
+      - DEL^2 U(X,Y) = F(X,Y)
+      
+    is solved on the unit square [0,1] x [0,1] using a grid of NX by
+    NX evenly spaced points.  The first and last points in each direction
+    are boundary points.
+    
+    The boundary conditions and F are set so that the exact solution is
+    
+      U(x,y) = sin ( pi * x * y )
+      
+    so that
+    
+      - DEL^2 U(x,y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
+      
+    The Jacobi iteration is repeatedly applied until convergence is detected.
+    
+    For convenience in writing the discretized equations, we assume that NX = NY.
+    
+  Licensing:
+  
+    This code is distributed under the GNU LGPL license. 
+    
+  Modified:
+  
+    14 December 2011
+    
+  Author:
+  
+    John Burkardt
+*/
+{
+  int converged;
+  double diff;
+  double dx;
+  double dy;
+  double error;
+  double f[NX][NY];
+  int i;
+  int id;
+  int itnew;
+  int itold;
+  int j;
+  int nx = NX;
+  int ny = NY;
+  double tolerance = 0.000001;
+  double u[NX][NY];
+  double u_norm;
+  double udiff[NX][NY];
+  double uexact[NX][NY];
+  double unew[NX][NY];
+  double unew_norm;
+  double wtime;
+  double x;
+  double y;
+  
+  dx = 1.0 / ( double ) ( nx - 1 );
+  dy = 1.0 / ( double ) ( ny - 1 );
+  /*
+  Print a message.
+*/
+  timestamp ( );
+  printf ( "\n" );
+  printf ( "POISSON_OPENMP:\n" );
+  printf ( "  C version\n" );
+  printf ( "  A program for solving the Poisson equation.\n" );
+  printf ( "\n" );
+  printf ( "  Use OpenMP for parallel execution.\n" );
+  printf ( "  The number of processors is %d\n", omp_get_num_procs ( ) );
+# pragma omp parallel
+  {
+    id = omp_get_thread_num ( );
+    if ( id == 0 )
+      {
+        printf ( "  The maximum number of threads is %d\n", omp_get_num_threads ( ) ); 
+      }
+  }
+  printf ( "\n" );
+  printf ( "  -DEL^2 U = F(X,Y)\n" );
+  printf ( "\n" );
+  printf ( "  on the rectangle 0 <= X <= 1, 0 <= Y <= 1.\n" );
+  printf ( "\n" );
+  printf ( "  F(X,Y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )\n" );
+  printf ( "\n" );
+  printf ( "  The number of interior X grid points is %d\n", nx );
+  printf ( "  The number of interior Y grid points is %d\n", ny );
+  printf ( "  The X grid spacing is %f\n", dx );
+  printf ( "  The Y grid spacing is %f\n", dy );
+  /*
+  Set the right hand side array F.
+*/
+  rhs ( nx, ny, f );
+  /*
+  Set the initial solution estimate UNEW.
+  We are "allowed" to pick up the boundary conditions exactly.
+*/
+# pragma omp for
+  for ( j = 0; j < ny; j++ )
+    {
+      for ( i = 0; i < nx; i++ )
+        {
+          if ( i == 0 || i == nx - 1 || j == 0 || j == ny - 1 )
+            {
+              unew[i][j] = f[i][j];
+            }
+          else
+            {
+              unew[i][j] = 0.0;
+            }
+        }
+    }
+  unew_norm = r8mat_rms ( nx, ny, unew );
+  /*
+  Set up the exact solution UEXACT.
+*/
+# pragma omp for
+  for ( j = 0; j < ny; j++ )
+    {
+      y = ( double ) ( j ) / ( double ) ( ny - 1 );
+      for ( i = 0; i < nx; i++ )
+        {
+          x = ( double ) ( i ) / ( double ) ( nx - 1 );
+          uexact[i][j] = u_exact ( x, y );
+        }
+    }
+  u_norm = r8mat_rms ( nx, ny, uexact );
+  printf ( "  RMS of exact solution = %g\n", u_norm );
+  /*
+  Do the iteration.
+*/
+  converged = 0;
+  
+  printf ( "\n" );
+  printf ( "  Step    ||Unew||     ||Unew-U||     ||Unew-Exact||\n" );
+  printf ( "\n" );
+  
+  for ( j = 0; j < ny; j++ )
+    {
+      for ( i = 0; i < nx; i++ )
+        {
+          udiff[i][j] = unew[i][j] - uexact[i][j];
+        }
+    }
+  error = r8mat_rms ( nx, ny, udiff );
+  printf ( "  %4d  %14g                  %14g\n", 0, unew_norm, error );
+  
+  wtime = omp_get_wtime ( );
+  
+  itnew = 0;
+  
+  for ( ; ; )
+    {
+      itold = itnew;
+      itnew = itold + 500;
+      /*
+  SWEEP carries out 500 Jacobi steps in parallel before we come
+  back to check for convergence.
+*/
+      sweep ( nx, ny, dx, dy, f, itold, itnew, u, unew );
+      /*
+  Check for convergence.
+*/
+      u_norm = unew_norm;
+      unew_norm = r8mat_rms ( nx, ny, unew );
+      
+# pragma omp for
+      for ( j = 0; j < ny; j++ )
+        {
+          for ( i = 0; i < nx; i++ )
+            {
+              udiff[i][j] = unew[i][j] - u[i][j];
+            }
+        }
+      diff = r8mat_rms ( nx, ny, udiff );
+      
+# pragma omp for
+      for ( j = 0; j < ny; j++ )
+        {
+          for ( i = 0; i < nx; i++ )
+            {
+              udiff[i][j] = unew[i][j] - uexact[i][j];
+            }
+        }
+      error = r8mat_rms ( nx, ny, udiff );
+      
+      printf ( "  %4d  %14g  %14g  %14g\n", itnew, unew_norm, diff, error );
+      
+      if ( diff <= tolerance )
+        {
+          converged = 1;
+          break;
+        }
+      
+    }
+  
+  if ( converged )
+    {
+      printf ( "  The iteration has converged.\n" );
+    }
+  else
+    {
+      printf ( "  The iteration has NOT converged.\n" );
+    }
+  
+  wtime = omp_get_wtime ( ) - wtime;
+  printf ( "\n" );
+  printf ( "  Elapsed seconds = %g\n", wtime );
+  /*
+  Terminate.
+*/
+  printf ( "\n" );
+  printf ( "POISSON_OPENMP:\n" );
+  printf ( "  Normal end of execution.\n" );
+  printf ( "\n" );
+  timestamp ( );
+  
+  return 0;
+}
+/******************************************************************************/
+
+double r8mat_rms ( int nx, int ny, double a[NX][NY] )
+
+/******************************************************************************/
+/*
+  Purpose:
+  
+    R8MAT_RMS returns the RMS norm of a vector stored as a matrix.
+    
+  Licensing:
+  
+    This code is distributed under the GNU LGPL license.
+    
+  Modified:
+  
+    01 March 2003
+    
+  Author:
+  
+    John Burkardt
+    
+  Parameters:
+  
+    Input, int NX, NY, the number of rows and columns in A.
+    
+    Input, double A[NX][NY], the vector.
+    
+    Output, double R8MAT_RMS, the root mean square of the entries of A.
+*/
+{
+  int i;
+  int j;
+  double v;
+  
+  v = 0.0;
+  
+# pragma omp for
+  for ( j = 0; j < ny; j++ )
+    {
+      for ( i = 0; i < nx; i++ )
+        {
+          v = v + a[i][j] * a[i][j];
+        }
+    }
+  v = sqrt ( v / ( double ) ( nx * ny )  );
+  
+  return v;
+}
+/******************************************************************************/
+
+void rhs ( int nx, int ny, double f[NX][NY] )
+
+/******************************************************************************/
+/*
+  Purpose:
+  
+    RHS initializes the right hand side "vector".
+    
+  Discussion:
+  
+    It is convenient for us to set up RHS as a 2D array.  However, each
+    entry of RHS is really the right hand side of a linear system of the
+    form
+    
+      A * U = F
+      
+    In cases where U(I,J) is a boundary value, then the equation is simply
+    
+      U(I,J) = F(i,j)
+      
+    and F(I,J) holds the boundary data.
+    
+    Otherwise, the equation has the form
+    
+      (1/DX^2) * ( U(I+1,J)+U(I-1,J)+U(I,J-1)+U(I,J+1)-4*U(I,J) ) = F(I,J)
+      
+    where DX is the spacing and F(I,J) is the value at X(I), Y(J) of
+    
+      pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
+      
+  Licensing:
+  
+    This code is distributed under the GNU LGPL license. 
+    
+  Modified:
+  
+    28 October 2011
+    
+  Author:
+  
+    John Burkardt
+    
+  Parameters:
+  
+    Input, int NX, NY, the X and Y grid dimensions.
+    
+    Output, double F[NX][NY], the initialized right hand side data.
+*/
+{
+  double fnorm;
+  int i;
+  int j;
+  double x;
+  double y;
+  /*
+  The "boundary" entries of F store the boundary values of the solution.
+  The "interior" entries of F store the right hand sides of the Poisson equation.
+*/
+# pragma omp for
+  for ( j = 0; j < ny; j++ )
+    {
+      y = ( double ) ( j ) / ( double ) ( ny - 1 );
+      for ( i = 0; i < nx; i++ )
+        {
+          x = ( double ) ( i ) / ( double ) ( nx - 1 );
+          if ( i == 0 || i == nx - 1 || j == 0 || j == ny - 1 )
+            {
+              f[i][j] = u_exact ( x, y );
+            }
+          else
+            {
+              f[i][j] = - uxxyy_exact ( x, y );
+            }
+        }
+    }
+  
+  fnorm = r8mat_rms ( nx, ny, f );
+  
+  printf ( "  RMS of F = %g\n", fnorm );
+  
+  return;
+}
+/******************************************************************************/
+
+void sweep ( int nx, int ny, double dx, double dy, double f[NX][NY], 
+             int itold, int itnew, double u[NX][NY], double unew[NX][NY] )
+
+/******************************************************************************/
+/*
+  Purpose:
+  
+   SWEEP carries out several steps of the Jacobi iteration.
+   
+  Discussion:
+  
+    Assuming DX = DY, we can approximate
+    
+      - ( d/dx d/dx + d/dy d/dy ) U(X,Y) 
+      
+    by
+    
+      ( U(i-1,j) + U(i+1,j) + U(i,j-1) + U(i,j+1) - 4*U(i,j) ) / dx / dy
+      
+    The discretization employed below will not be correct in the general
+    case where DX and DY are not equal.  It's only a little more complicated
+    to allow DX and DY to be different, but we're not going to worry about 
+    that right now.
+    
+  Licensing:
+  
+    This code is distributed under the GNU LGPL license. 
+    
+  Modified:
+  
+    14 December 2011
+    
+  Author:
+  
+    John Burkardt
+    
+  Parameters:
+  
+    Input, int NX, NY, the X and Y grid dimensions.
+    
+    Input, double DX, DY, the spacing between grid points.
+    
+    Input, double F[NX][NY], the right hand side data.
+    
+    Input, int ITOLD, the iteration index on input.
+    
+    Input, int ITNEW, the desired iteration index
+    on output.
+    
+    Input, double U[NX][NY], the solution estimate on 
+    iteration ITNEW-1.
+    
+    Input/output, double UNEW[NX][NY], on input, the solution 
+    estimate on iteration ITOLD.  On output, the solution estimate on 
+    iteration ITNEW.
+*/
+{
+  int i;
+  int it;
+  int j;
+  
+# pragma omp parallel \
+  shared ( dx, dy, f, itnew, itold, nx, ny, u, unew ) \
+  private ( i, it, j )
+  
+  for ( it = itold + 1; it <= itnew; it++ )
+    {
+      /*
+  Save the current estimate.
+*/
+# pragma omp for
+      for ( j = 0; j < ny; j++ )
+        {
+          for ( i = 0; i < nx; i++ )
+            {
+              u[i][j] = unew[i][j];
+            }
+        }
+      /*
+  Compute a new estimate.
+*/
+#pragma push
+#pragma O3
+# pragma omp for
+      for ( j = 0; j < ny; j++ )
+        {
+          for ( i = 0; i < nx; i++ )
+            {
+              if ( i == 0 || j == 0 || i == nx - 1 || j == ny - 1 )
+                {
+                  unew[i][j] = f[i][j];
+                }
+              else
+                { 
+                  //            #pragma optimize("", on) 
+                  unew[i][j] = 0.25 * ( 
+                        u[i-1][j] + u[i][j+1] + u[i][j-1] + u[i+1][j] + f[i][j] * dx * dy );
+                }
+            }
+        }
+#pragma pop
+    }
+  return;
+}
+/******************************************************************************/
+
+void timestamp ( )
+
+/******************************************************************************/
+/*
+  Purpose:
+  
+    TIMESTAMP prints the current YMDHMS date as a time stamp.
+    
+  Example:
+  
+    31 May 2001 09:45:54 AM
+    
+  Licensing:
+  
+    This code is distributed under the GNU LGPL license. 
+    
+  Modified:
+  
+    24 September 2003
+    
+  Author:
+  
+    John Burkardt
+    
+  Parameters:
+  
+    None
+*/
+{
+# define TIME_SIZE 40
+  
+  static char time_buffer[TIME_SIZE];
+  const struct tm *tm;
+  time_t now;
+  
+  now = time ( NULL );
+  tm = localtime ( &now );
+  
+  strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );
+  
+  printf ( "%s\n", time_buffer );
+  
+  return;
+# undef TIME_SIZE
+}
+/******************************************************************************/
+
+double u_exact ( double x, double y )
+
+/******************************************************************************/
+/*
+  Purpose:
+  
+    U_EXACT evaluates the exact solution.
+    
+  Licensing:
+  
+    This code is distributed under the GNU LGPL license.
+    
+  Modified:
+  
+    25 October 2011
+    
+  Author:
+  
+    John Burkardt
+    
+  Parameters:
+  
+    Input, double X, Y, the coordinates of a point.
+    
+    Output, double U_EXACT, the value of the exact solution 
+    at (X,Y).
+*/
+{
+  double r8_pi = 3.141592653589793;
+  double value;
+  
+  value = sin ( r8_pi * x * y );
+  
+  return value;
+}
+/******************************************************************************/
+
+double uxxyy_exact ( double x, double y )
+
+/******************************************************************************/
+/*
+  Purpose:
+  
+    UXXYY_EXACT evaluates ( d/dx d/dx + d/dy d/dy ) of the exact solution.
+    
+  Licensing:
+  
+    This code is distributed under the GNU LGPL license.
+    
+  Modified:
+  
+    25 October 2011
+    
+  Author:
+  
+    John Burkardt
+    
+  Parameters:
+  
+    Input, double X, Y, the coordinates of a point.
+    
+    Output, double UXXYY_EXACT, the value of 
+    ( d/dx d/dx + d/dy d/dy ) of the exact solution at (X,Y).
+*/
+{
+  double r8_pi = 3.141592653589793;
+  double value;
+  
+  value = - r8_pi * r8_pi * ( x * x + y * y ) * sin ( r8_pi * x * y );
+  
+  return value;
+}
+# undef NX
+# undef NY