Init

README.md

-# poisson
+# A Program for the Poisson Equation in a Rectangle
 
-Poisson computes an approximate solution to the Poisson equation in a rectangle.
+This program computes an approximate
+solution to the Poisson equation in a rectangular region.
 
+The version of Poisson's equation being solved here is
+
+            - ( d/dx d/dx + d/dy d/dy ) U(x,y) = F(x,y)
+          
+
+over the rectangle 0 <= X <= 1, 0 <= Y <= 1, with exact
+solution
+
+            U(x,y) = sin ( pi * x * y )
+          
+
+so that
+
+            F(x,y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )
+          
+
+and with Dirichlet boundary conditions along the lines x = 0, x = 1, y =
+0 and y = 1. (The boundary conditions will actually be zero in this
+case, but we write up the problem as though we didn't know that, which
+makes it easy to change the problem later.)
+
+We compute an approximate solution by discretizing the geometry,
+assuming that DX = DY, and approximating the Poisson operator by
+
+            ( U(i-1,j) + U(i+1,j) + U(i,j-1) + U(i,j+1) - 4*U(i,j) ) / dx /dy
+          
+
+Along with the boundary conditions at the boundary nodes, we have a
+linear system for U. We can apply the Jacobi iteration to estimate the
+solution to the linear system.
+
+
+## STEP 0: Create git repository (10%)
+
+Your code should be forked from this repository and hosted on code.it4i.cz as a private project with access for all teachers.
+
+## STEP 1:  Building the library (10%)
+
+Provide compilation script for your application (the script should run independently on a current path). Script should load all necessary modules and call `cmake`.
+
+## STEP 2:  Analysis of the application (10%)
+
+Use `Arm map` (module Forge) to analyze a sequential run of your application with given use case (`tests/large`). Identify the most time consuming regions that can be parallelized by OpenMP.
+
+## STEP 3:  Use OpenMP to run the application in parallel (10%)
+
+Put OpenMP pragmas to a correct positions with appropriate variables visibility in order to utilize more threads effectively.
+
+## STEP 4:  Test the correctness of the code (10%)
+
+Create script that automatically check correctness of your application for at least 3 different test cases. Comparison can be implemented as comparison of outputs of sequential and parallel runs.
+
+## STEP 5:  Test the behavior of the code on the Karolina cluster (40%)
+
+1. Implement time measurement for all parallel regions using omp_get_wtime(). 
+2. Create script for run strong scalability measurement (PBS script).
+3. Evaluate strong scalability of measured regions up to 128 threads and different thread affinity (compact, scatter, balanced, none).
+4. Evaluate behaviour for domain specific parameters of your applications (project dependent, maybe none). (Blocking implementation for this project)
+5. Prepare charts for all measurements.
+
+## STEP 6:  Presentation of your project (10%)
+
+Prepare presentation in form of slides (pptx, pdf). The slides should address all topics requested above.

src/main.cpp

+# include "poisson.h"
+
+//****************************************************************************80
+
+int main(int argc, char* argv[])
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//    MAIN is the main program for POISSON_SERIAL.
+//
+//  Discussion:
+//
+//    POISSON_SERIAL is a program for solving the Poisson problem.
+//
+//    This program runs serially.  Its output is used as a benchmark for
+//    comparison with similar programs run in a parallel environment.
+//
+//    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:
+//
+//    26 October 2011
+//
+//  Author:
+//
+//    John Burkardt
+//
+{
+	bool converged;
+	double diff;
+	double dx;
+	double dy;
+	double error;
+	double** f;
+	int i;
+	int it;
+	int it_max = 100000;
+	int j;
+	int nx = 256;
+	int ny = 128;
+	double tolerance = 0.000001;
+	double** u;
+	double u_norm;
+	double** udiff;
+	double** uexact;
+	double** unew;
+	double unew_norm;
+	double x;
+	double y;
+
+	//Check args
+	if (argc < 2) {
+		printf("USAGE:\tpoisson 0 (tiny)\n\tpoisson 1 (large)\n");
+		return 0;
+	}
+
+	// Read Input
+	int type = atoi(argv[1]);
+	//tiny
+	if (type == 0) {
+		nx = 11;
+		ny = 11;
+		it_max = 1000;
+	}
+	//large
+	if (type == 1) {
+		nx = 256;
+		ny = 128;
+		it_max = 100000;
+	}
+
+	//Alloc memory
+	f = new double* [nx];
+	u = new double* [nx];
+	udiff = new double* [nx];
+	uexact = new double* [nx];
+	unew = new double* [nx];
+	for (i = 0; i < nx; i++)
+	{
+		f[i] = new double[ny];
+		u[i] = new double[ny];
+		udiff[i] = new double[ny];
+		uexact[i] = new double[ny];
+		unew[i] = new double[ny];
+	}
+
+	dx = 1.0 / (double)(nx - 1);
+	dy = 1.0 / (double)(ny - 1);
+	//
+	//  Print a message.
+	//
+	cout << "\n";
+	cout << "POISSON_SERIAL:\n";
+	cout << "  C++ version\n";
+	cout << "  A program for solving the Poisson equation.\n";
+	cout << "\n";
+	cout << "  -DEL^2 U = F(X,Y)\n";
+	cout << "\n";
+	cout << "  on the rectangle 0 <= X <= 1, 0 <= Y <= 1.\n";
+	cout << "\n";
+	cout << "  F(X,Y) = pi^2 * ( x^2 + y^2 ) * sin ( pi * x * y )\n";
+	cout << "\n";
+	cout << "  The number of interior X grid points is " << nx << "\n";
+	cout << "  The number of interior Y grid points is " << ny << "\n";
+	cout << "  The X grid spacing is " << dx << "\n";
+	cout << "  The Y grid spacing is " << dy << "\n";
+	//
+	//  Initialize the data.
+	//
+	rhs(nx, ny, f);
+	//
+	//  Set the initial solution estimate.
+	//  We are "allowed" to pick up the boundary conditions exactly.
+	//
+	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.
+	//
+	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);
+	cout << "  RMS of exact solution = " << u_norm << "\n";
+	//
+	//  Do the iteration.
+	//
+	converged = false;
+
+	cout << "\n";
+	cout << "  Step    ||Unew||     ||Unew-U||     ||Unew-Exact||\n";
+	cout << "\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);
+	cout << "  " << setw(4) << 0
+		<< "  " << setw(14) << unew_norm
+		<< "  " << "              "
+		<< "  " << setw(14) << error << "\n";
+
+	for (it = 1; it <= it_max; it++)
+	{
+		for (j = 0; j < ny; j++)
+		{
+			for (i = 0; i < nx; i++)
+			{
+				u[i][j] = unew[i][j];
+			}
+		}
+		//
+		//  UNEW is derived from U by one Jacobi step.
+		//
+		sweep(nx, ny, dx, dy, f, u, unew);
+		//
+		//  Check for convergence.
+		//
+		u_norm = unew_norm;
+		unew_norm = r8mat_rms(nx, ny, unew);
+
+		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);
+
+		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);
+
+		if (diff <= tolerance)
+		{
+			cout << "  " << setw(4) << it
+				<< "  " << setw(14) << unew_norm
+				<< "  " << setw(14) << diff
+				<< "  " << setw(14) << error << "\n";
+
+			converged = true;
+			break;
+		}
+	}
+
+	if (converged)
+	{
+		cout << "  The iteration has converged.\n";
+	}
+	else
+	{
+		cout << "  The iteration has NOT converged.\n";
+	}
+	//
+	//  Terminate.
+	//
+	cout << "\n";
+	cout << "POISSON_SERIAL:\n";
+	cout << "  Normal end of execution.\n";
+	cout << "\n";
+
+	//Free memory
+	for (i = 0; i < nx; i++)
+	{
+		delete[] f[i];
+		delete[] u[i];
+		delete[] udiff[i];
+		delete[] uexact[i];
+		delete[] unew[i];
+	}
+
+	delete[] f;
+	delete[] u;
+	delete[] udiff;
+	delete[] uexact;
+	delete[] unew;
+
+	return 0;
+}
+//****************************************************************************80

src/poisson.cpp

+# include "poisson.h"
+
+double r8mat_rms ( int nx, int ny, double **a )
+
+//****************************************************************************80
+//
+//  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;
+
+  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;
+}
+//****************************************************************************80
+
+void rhs ( int nx, int ny, double **f )
+
+//****************************************************************************80
+//
+//  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.
+//
+  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 );
+
+  cout << "  RMS of F = " << fnorm << "\n";
+
+  return;
+}
+//****************************************************************************80
+
+void sweep ( int nx, int ny, double dx, double dy, double **f, 
+  double **u, double **unew )
+
+//****************************************************************************80
+//
+//  Purpose:
+//
+//   SWEEP carries out one step of the Jacobi iteration.
+//
+//  Licensing:
+//
+//    This code is distributed under the GNU LGPL license. 
+//
+//  Modified:
+//
+//    26 October 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, double U[NX][NY], the previous solution estimate.
+//
+//    Output, double UNEW[NX][NY], the updated solution estimate.
+//
+{
+  int i;
+  int j;
+
+  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
+      { 
+        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 );
+      }
+    }
+  }
+  return;
+}
+
+//****************************************************************************80
+
+double u_exact ( double x, double y )
+
+//****************************************************************************80
+//
+//  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 pi = 3.141592653589793;
+  double value;
+
+  value = sin ( pi * x * y );
+
+  return value;
+}
+//****************************************************************************80
+
+double uxxyy_exact ( double x, double y )
+
+//****************************************************************************80
+//
+//  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 pi = 3.141592653589793;
+  double value;
+
+  value = - pi * pi * ( x * x + y * y ) * sin ( pi * x * y );
+
+  return value;
+}

src/poisson.h

+#pragma once
+
+# include <cstdlib>
+# include <iostream>
+# include <iomanip>
+# include <ctime>
+# include <cmath>
+
+using namespace std;
+
+double r8mat_rms ( int m, int n, double **a );
+void rhs ( int nx, int ny, double **f );
+void sweep ( int nx, int ny, double dx, double dy, double **f, 
+  double **u, double **unew);
+double u_exact ( double x, double y );
+double uxxyy_exact ( double x, double y );
\ No newline at end of file