Changes

Jump to: navigation, search

GPUSquad

1,409 bytes added, 12:12, 7 April 2018
→‎Idea 1 - Jacobi Method for 2D Poisson Problem
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. The error_and_output function was rewritten as outputImage to get rid of terminal and file I/O during calculations, and only output the final image at the end of the calculations (this would be necessary for the parallel versions, so the changes were applied to the serial code as well to make things fair).
<source>
#include <fstream>
#include <cmath>
#include <chrono>
using namespace std;
// Set grid size and number of iterations
const int save_iters=20;const int total_iters=5000;const int error_every=2;const int m=33032,n=3301024;const double float xmin=-1,xmax=1;const double float ymin=-1,ymax=1;
// Compute useful constants
const double float pi=3.1415926535897932384626433832795;const double float omega=2/(1+sin(2*pi/n));const double float dx=(xmax-xmin)/(m-1);const double float dy=(ymax-ymin)/(n-1);const double float dxxinv=1/(dx*dx);const double float dyyinv=1/(dy*dy);const double float dcent=1/(2*(dxxinv+dyyinv));
// Input function
inline double float f(int i,int j) { double float x=xmin+i*dx,y=ymin+j*dy; return abs(x)>0.5||abs(y)>0.5?0:1;
}
// Common output and error routine
void output_and_error(char* filename,double float *a,const int sn) { // Computes the error if sn%error every==0 if(sn%error_every==0) { double float 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 dojacobioutputImage(int i, int j, int ij, int k, double error, double u[]char* filename, double v[], double z, doublefloat * a, double* b) { // Computes the error if sn%error every==0
// Set initial guess to be identically zero
for(ij=0;ij<m*n;ij++) u[ij]=v[ij]=0;
output_and_error("jacobi out",u,0);
// Carry out Jacobi iterations for(k=1;kSaves the matrix if sn<=total_iters;k++) {save iters // Alternately flip input and output matrices if (k%2 int i, j, ij =0, ds =0sizeof(float) {a=u;b=v float x, y, data_float;} else {a=vconst char *pfloat;b pfloat =u(const char*)&data_float;}
// Compute Jacobi iteration ofstream outfile; for(j=1;j<n-1 static char fname[256];j++) { for sprintf(i=1;i<m-1;i++fname, "%s.%d", filename, 101) { ij=i+m*j; a[ij]=(f outfile.open(ifname,j)+dxxinv*(b[ij-1]+b[ij+1])fstream::out +dyyinv*(b[ij-m]+b[ij+m] | fstream::trunc | fstream::binary))*dcent; } }
data_float = m; outfile.write(pfloat, ds);  for (i = 0; i < m; i++) { x = xmin + i*dx; data_float = x; outfile.write(pfloat, ds);//will need to be deferred }  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, float error, float u[], float v[], float z, float* a, float* b) {  // Set initial guess to be identically zero for (ij = 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); } outputImage("jacobi out", a);
}
int main() {
int i,j,ij,k; double //float error,u[m*n],v[m*n],z; double float error, z; float *a,*b; float *u; float *v; u = new float[m*n]; v = new float[m*n];
std::chrono::steady_clock::time_point ts, te; ts = std::chrono::steady_clock::now(); dojacobi(i, j, ij, k, error, u, v, z, a, b); delete[] u; delete[] v;
te = std::chrono::steady_clock::now(); std::chrono::steady_clock::duration duration = te - ts; auto ms = std::chrono::duration_cast<std::chrono::milliseconds>(duration); std::cout << "Serial Code Time: " << ms.count() << " ms" << std::endl;  return 0;
}
 
</source>
Tsarkarcd
93
edits
Retrieved from "https://wiki.cdot.senecapolytechnic.ca/wiki/Special:MobileDiff/133031"

Navigation menu