==== Idea 1 - Jacobi Method for 2D Poisson Problem ==== This topic is based on the following pdf: https://math.berkeley.edu/~wilken/228A.F07/chr_lecture.pdf Background: Poisson's equation according to Wikipedia: "In mathematics, Poisson's equation is a partial differential equation of elliptic type with broad utility in mechanical engineering and theoretical physics. It arises, for instance, to describe the potential field caused by a given charge or mass density distribution; with the potential field known, one can then calculate gravitational or electrostatic field. It is a generalization of Laplace's equation, which is also frequently seen in physics." According to the intro in the pdf above, finite-difference and finite-element methods are the solution techniques of choice for solving elliptic PDE problems. Regardless of which type of technique you choose, you will end up with sets of linear relationships between various variables, and will need to solve for the solution u_k which satisfies all the linear relationships prescribed by the PDE.  This can be written, according to the pdf, as a matrix "Au = b, where we wish to find a solution u, given that A is a matrix capturing the differentiation operator, and b corresponds to any source or boundary terms". The Jacobi method can be used to solve this matrix, and is used in the code sample later. Jacobi method according to Wikipedia: "In numerical linear algebra, the Jacobi method (or Jacobi iterative method[1]) is an algorithm for determining the solutions of a diagonally dominant system of linear equations. Each diagonal element is solved for, and an approximate value is plugged in. The process is then iterated until it converges." <nowiki>*******************</nowiki> CODE - Based on PDF <nowiki>*******************</nowiki> The following code is based on the jacobi.cc and common.cc code at the bottom of the source PDF, but there are some changes for these tests. The files use .cpp extensions instead of .cc, the code in common.cc was added to jacobi.cpp so it's a single file, and code for carrying out the Jacobi iterations in the main were placed into a seperate function called dojacobi() for gprof profiling. The std::cout line inside the output_and_error function was commented out to avoid excessive printing to the terminal. Compilation on Linux:g++ -std=c++0x -O2 -g -pg -o jacobi jacobi.cpp  <source>// Load standard libraries#include <cstdio>#include <cstdlib>#include <iostream>#include <fstream>#include <cmath>using namespace std; // Set grid size and number of iterationsconst int save_iters=20;const int total_iters=5000;const int error_every=2;const int m=330,n=330;const double xmin=-1,xmax=1;const double ymin=-1,ymax=1; // Compute useful constantsconst double pi=3.1415926535897932384626433832795;const double omega=2/(1+sin(2*pi/n));const double dx=(xmax-xmin)/(m-1);const double dy=(ymax-ymin)/(n-1);const double dxxinv=1/(dx*dx);const double dyyinv=1/(dy*dy);const double dcent=1/(2*(dxxinv+dyyinv)); // Input functioninline double f(int i,int j) { double x=xmin+i*dx,y=ymin+j*dy; return abs(x)>0.5||abs(y)>0.5?0:1;} // Common output and error routinevoid output_and_error(char* filename,double *a,const int sn) { // Computes the error if sn%error every==0 if(sn%error_every==0) { double z,error=0;int ij; for(int j=1;j<n-1;j++) { for(int i=1;i<m-1;i++) { ij=i+m*j; z=f(i,j)-a[ij]*(2*dxxinv+2*dyyinv) +dxxinv*(a[ij-1]+a[ij+1]) +dyyinv*(a[ij-m]+a[ij+m]); error+=z*z; } } //cout << sn << " " << error*dx*dy << endl; }  // Saves the matrix if sn<=save iters if(sn<=save_iters) { int i,j,ij=0,ds=sizeof(float); float x,y,data_float;const char *pfloat; pfloat=(const char*)&data_float;  ofstream outfile; static char fname[256]; sprintf(fname,"%s.%d",filename,sn); outfile.open(fname,fstream::out | fstream::trunc | fstream::binary);  data_float=m;outfile.write(pfloat,ds);  for(i=0;i<m;i++) { x=xmin+i*dx; data_float=x;outfile.write(pfloat,ds); }  for(j=0;j<n;j++) { y=ymin+j*dy; data_float=y; outfile.write(pfloat,ds);  for(i=0;i<m;i++) { data_float=a[ij++]; outfile.write(pfloat,ds); } }  outfile.close(); }} void dojacobi(int i, int j, int ij, int k, double error, double u[], double v[], double z, double* a, double* b) {  // Set initial guess to be identically zero for(insert topicij=0;ij<m*n;ij++) u[ij]=v[ij]=0; output_and_error("jacobi out",u,0);  // Carry out Jacobi iterations for(k=1;k<=total_iters;k++) { // Alternately flip input and output matrices if (k%2==0) {a=u;b=v;} else {a=v;b=u;}  // Compute Jacobi iteration for(j=1;j<n-1;j++) { for(i=1;i<m-1;i++) { ij=i+m*j; a[ij]=(f(i,j)+dxxinv*(b[ij-1]+b[ij+1]) +dyyinv*(b[ij-m]+b[ij+m]))*dcent; } }  // Save and compute error if necessary output_and_error("jacobi out",a,k); } } int main() { int i,j,ij,k; double error,u[m*n],v[m*n],z; double *a,*b;  dojacobi(i, j, ij, k, error, u, v, z, a, b);  return 0;}</source> <nowiki>********************</nowiki> Inputs / Performance <nowiki>********************</nowiki> In a 2D PDE such as the 2D Poisson problem, the matrix will have m * n gridpoints. At m = 33, n = 33, iterations = 5000: <source>Flat profile: Each sample counts as 0.01 seconds. % cumulative self self total time seconds seconds calls Ts/call Ts/call name 100.00 0.05 0.05 dojacobi(int, int, int, int, double, double*, double*, double, double*, double*) 0.00 0.05 0.00 5001 0.00 0.00 output_and_error(char*, double*, int) 0.00 0.05 0.00 1 0.00 0.00 _GLOBAL__sub_I__Z16output_and_errorPcPdi </source> At m = 165, n = 165, iterations =5000: <source>Flat profile: Each sample counts as 0.01 seconds. % cumulative self self total time seconds seconds calls us/call us/call name 99.15 1.16 1.16 dojacobi(int, int, int, int, double, double*, double*, double, double*, double*) 0.85 1.17 0.01 5001 2.00 2.00 output_and_error(char*, double*, int) 0.00 1.17 0.00 1 0.00 0.00 _GLOBAL__sub_I__Z16output_and_errorPcPdi </source> At m =330, n =330, iterations =5000: <source>Flat profile: Each sample counts as 0.01 seconds. % cumulative self self total time seconds seconds calls us/call us/call name 99.43 5.26 5.26 dojacobi(int, int, int, int, double, double*, double*, double, double*, double*) 0.57 5.29 0.03 5001 6.00 6.00 output_and_error(char*, double*, int) 0.00 5.29 0.00 1 0.00 0.00 _GLOBAL__sub_I__Z16output_and_errorPcPdi </source> The hotspot seems to clearly be the triple for-loop in the Jacobi iterations code of the dojacobi() function. I believe these matrix calculations could be parallelized for improved performance.