Merged all the code + minor bugfixes for jacobi

mpi-assign4/merged_code/.#Makefile

+bogischv@vmhalle0695.vdi.informatik.tu-muenchen.de.9776:1339164629
\ No newline at end of file

mpi-assign4/merged_code/Makefile

+# Intel compiler
+CC = mpicc
+CFLAGS =-O3
+
+MPICC = mpicc
+
+all: heat
+
+heat : heat.o input.o misc.o timing.o relax_gauss.o relax_jacobi.o
+	$(CC) $(CFLAGS) -o $@ $+ -lm 
+
+%.o : %.c %.h
+	$(CC) $(CFLAGS) -c -o $@ $<
+
+clean:
+	rm -f *.o heat *~ *.ppm
+
+remake : clean all

mpi-assign4/merged_code/heat.c

+/*
+ * heat.h
+ *
+ * Iterative solver for heat distribution
+ */
+
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#define FALSE 0
+#define TRUE 1
+
+#include "input.h"
+#include "heat.h"
+#include "timing.h"
+
+
+void usage( char *s )
+{
+	fprintf(stderr, 
+			"Usage: %s <input file> [result file]\n\n", s);
+}
+
+
+int main( int argc, char *argv[] )
+{
+	unsigned iter;
+	FILE *infile, *resfile;
+	char *resfilename;
+	int periods[2]={FALSE, FALSE};
+	MPI_Comm comm2d;
+
+	int nthreads, rank;
+	int thdx, thdy;
+
+	double *tmp;
+
+	// algorithmic parameters
+	algoparam_t param;
+	int np;
+
+	double runtime, flop;
+	double residual;
+	
+	int coords[2];
+
+	MPI_Init(&argc, &argv);
+	MPI_Comm_size(MPI_COMM_WORLD, &nthreads);
+	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+	
+	// check arguments
+	if( argc < 2 )
+	{
+		usage( argv[0] );
+		MPI_Finalize();
+		return 1;
+	}
+
+	// check input file
+	if( !(infile=fopen(argv[1], "r"))  ) 
+	{
+		fprintf(stderr, 
+				"\nError: Cannot open \"%s\" for reading.\n\n", argv[1]);
+
+		usage(argv[0]);
+		MPI_Finalize();
+		return 1;
+	}
+
+    // check input
+ 
+    if( !read_input(infile, &param) )
+    {
+	fprintf(stderr, "\nError: Error parsing input file.\n\n");
+      
+	usage(argv[0]);
+	MPI_Finalize();
+	return 1;
+    }
+
+
+    MPI_Cart_create(MPI_COMM_WORLD, 2, param.thread_dims, periods, FALSE, &comm2d);
+    MPI_Cart_coords(comm2d, rank, 2, coords);
+
+    if(rank==0) 
+	print_params(&param);
+
+    param.u     = 0;
+    param.uhelp = 0;
+    param.sendbuf_left =0;
+    param.sendbuf_right = 0;
+    param.recbuf_left = 0;
+    param.recbuf_right =0;
+    param.sendbuf_top =0;
+    param.sendbuf_bottom = 0;
+    param.recbuf_top = 0;
+    param.recbuf_bottom =0;
+
+    param.act_res = param.initial_res;
+ 
+    // loop over different resolutions
+    while(1) {
+
+	// free allocated memory of previous experiment
+	if (param.u != 0)
+	    finalize(&param);
+
+	if( !initialize(&param, coords) )
+	{
+	    fprintf(stderr, "Error in Jacobi initialization.\n\n");
+
+	    usage(argv[0]);
+	}
+
+	fprintf(stderr, "Resolution: %5u\r", param.act_res);
+
+	// full size (param.act_res are only the inner points)
+	np = param.act_res + 2;
+    
+	// starting time
+	runtime = MPI_Wtime();
+	residual = 999999999;
+ 	
+	iter = 0;
+	while(1) {
+
+	    switch( param.algorithm ) {
+
+		case 0: // JACOBI
+      
+		    residual = relax_jacobi_return_residual(&param, coords , comm2d);
+			tmp = param.u;
+			param.u = param.uhelp;
+			param.uhelp = tmp;
+		    break;
+
+		case 1: // GAUSS
+		    if (param.thread_dims[0] == 1)
+		    	residual = relax_gauss_return_residual_1dy(&param, param.maxiter, coords , comm2d);
+		    else if (param.thread_dims[1] == 1)
+                        residual = relax_gauss_return_residual_1dx(&param, param.maxiter, coords , comm2d);
+		    else 
+			residual = relax_gauss_return_residual_2d(&param, param.maxiter, coords , comm2d);
+		    break;
+	    }
+	    if (param.algorithm == 0) // JACOBI
+		iter += 1;
+	    else // GAUSS  
+	    	iter+= param.maxiter;
+
+	    // solution good enough ?
+	    if (residual < 0.000005) break;
+
+	    // max. iteration reached ? (no limit with maxiter=0)
+	    if (param.maxiter>0 && iter>=param.maxiter) break;
+	    
+	    if (iter % 100 == 0)
+		fprintf(stderr, "residual %f, %d iterations\n", residual, iter);
+	}
+	
+	// Flop count after <i> iterations
+	flop = iter * 7.0 * param.act_res * param.act_res;
+	// stopping time
+	runtime = MPI_Wtime() - runtime;
+
+	if (rank == 0)
+	{
+		fprintf(stderr, "Resolution: %5u, ", param.act_res);
+		fprintf(stderr, "Time: %04.3f ", runtime);
+		fprintf(stderr, "(%3.3f GFlop => %6.2f MFlop/s, ", 
+			flop/1000000000.0,
+			flop/runtime/1000000);
+		fprintf(stderr, "residual %f, %d iterations)\n", residual, iter);
+
+		// for plot...
+		printf("%5d %f\n", param.act_res, flop/runtime/1000000);
+	}
+
+	if (param.act_res + param.res_step_size > param.max_res) break;
+	param.act_res += param.res_step_size;
+    }
+
+    finalize( &param );
+    MPI_Finalize();
+    return 0;
+}

mpi-assign4/merged_code/heat.h

+/*
+ * heat.h
+ *
+ * Global definitions for the iterative solver
+ */
+
+#ifndef JACOBI_H_INCLUDED
+#define JACOBI_H_INCLUDED
+
+#include <stdio.h>
+#include <mpi.h>
+
+// configuration
+
+typedef struct
+{
+    float posx;
+    float posy;
+    float range;
+    float temp;
+}
+heatsrc_t;
+
+typedef struct
+{
+    unsigned maxiter;       // maximum number of iterations
+    unsigned act_res;
+    unsigned max_res;       // spatial resolution
+    unsigned initial_res;
+    unsigned res_step_size;
+    unsigned blocksize_x;
+    unsigned blocksize_y;
+    int algorithm;          // 0=>Jacobi, 1=>Gauss
+
+  
+    double *u, *uhelp;
+    double *sendbuf_left, *recbuf_left, *sendbuf_right, *recbuf_right;
+    double *sendbuf_top, *recbuf_top, *sendbuf_bottom, *recbuf_bottom;
+
+    unsigned   numsrcs;     // number of heat sources
+    heatsrc_t *heatsrcs;
+    int thread_dims[2]; 	//x*y dimensions of thread
+
+}
+algoparam_t;
+
+
+// function declarations
+
+//Get neighbour ranks
+int get_rank_from_2dcoords(int coord_x, int coord_y, int *thread_dims, MPI_Comm comm2d);
+
+// misc.c
+int initialize( algoparam_t *param, int *coords );
+int finalize( algoparam_t *param );
+
+// Gauss-Seidel: relax_gauss.c
+double relax_gauss_return_residual_1dy( algoparam_t *param, int interleaving_count, int *coords, MPI_Comm comm2d);
+double relax_gauss_return_residual_1dx( algoparam_t *param, int interleaving_count, int *coords, MPI_Comm comm2d);
+double relax_gauss_return_residual_2d( algoparam_t *param, int interleaving_count, int *coords, MPI_Comm comm2d);
+// Jacobi: relax_jacobi.c
+
+double relax_jacobi_return_residual( algoparam_t *param, int *coords, MPI_Comm comm2d ); 
+
+
+#endif // JACOBI_H_INCLUDED

mpi-assign4/merged_code/input.c

+//
+// input.c
+// 
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "input.h"
+
+#define BUFSIZE 100
+
+int read_input( FILE *infile, algoparam_t *param )
+{
+  int i, n;
+  char buf[BUFSIZE];
+
+  fgets(buf, BUFSIZE, infile);
+  n = sscanf( buf, "%d %d", &(param->thread_dims[0]), &(param->thread_dims[1]) );
+  if( n!=2 )
+    return 0;
+
+  fgets(buf, BUFSIZE, infile);
+  n = sscanf( buf, "%u", &(param->maxiter) );
+  if( n!=1 )
+    return 0;
+
+  fgets(buf, BUFSIZE, infile);
+  n = sscanf( buf, "%u", &(param->initial_res) );
+  if( n!=1 )
+    return 0;
+
+  fgets(buf, BUFSIZE, infile);
+  n = sscanf( buf, "%u", &(param->max_res) );
+  if( n!=1 )
+    return 0;
+
+  fgets(buf, BUFSIZE, infile);
+  n = sscanf( buf, "%u", &(param->res_step_size) );
+  if( n!=1 )
+    return 0;
+
+
+  fgets(buf, BUFSIZE, infile);
+  n = sscanf(buf, "%d", &(param->algorithm) );
+  if( n!=1 )
+    return 0;
+
+  fgets(buf, BUFSIZE, infile);
+  n = sscanf(buf, "%u", &(param->numsrcs) );
+  if( n!=1 )
+    return 0;
+
+  (param->heatsrcs) = 
+    (heatsrc_t*) malloc( sizeof(heatsrc_t) * (param->numsrcs) );
+  
+  for( i=0; i<param->numsrcs; i++ )
+    {
+      fgets(buf, BUFSIZE, infile);
+      n = sscanf( buf, "%f %f %f %f",
+		  &(param->heatsrcs[i].posx),
+		  &(param->heatsrcs[i].posy),
+		  &(param->heatsrcs[i].range),
+		  &(param->heatsrcs[i].temp) );
+
+      if( n!=4 )
+	return 0;
+    }
+
+  return 1;
+}
+
+
+void print_params( algoparam_t *param )
+{
+  int i;
+
+  fprintf(stderr, "Number of Threads: (%d x %d)\n", param->thread_dims[0], param->thread_dims[1]);
+  fprintf(stderr, "Resolutions       : (%u, %u, ... %u)\n",
+	  param->initial_res,
+	  param->initial_res + param->res_step_size,
+	  param->max_res);
+  fprintf(stderr, "Iterations        : %u\n", param->maxiter);
+  fprintf(stderr, "Algorithm         : %d (%s)\n",
+	  param->algorithm,
+	  (param->algorithm == 0) ? "Jacobi":"Gauss-Jacobi" );
+  fprintf(stderr, "Num. Heat sources : %u\n", param->numsrcs);
+
+  for( i=0; i<param->numsrcs; i++ )
+    {
+      fprintf(stderr, "  %2d: (%2.2f, %2.2f) %2.2f %2.2f \n",
+	     i+1,
+	     param->heatsrcs[i].posx,
+	     param->heatsrcs[i].posy,
+	     param->heatsrcs[i].range,
+	     param->heatsrcs[i].temp );
+    }
+}

mpi-assign4/merged_code/input.h

+//
+// input.h
+//
+
+#ifndef INPUT_H_INCLUDED
+#define INPUT_H_INCLUDED
+
+#include <stdio.h>
+
+#include "heat.h"
+
+int read_input( FILE *infile, algoparam_t *param );
+void print_params( algoparam_t *param );
+
+
+#endif // INPUT_H_INCLUDED

mpi-assign4/merged_code/misc.c

+/*
+ * misc.c
+ *
+ * Helper functions for
+ * - initialization
+ * - finalization,
+ * - writing out a picture
+ */
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+#include <float.h>
+
+#include "heat.h"
+
+/*
+ * Initialize the iterative solver
+ * - allocate memory for matrices
+ * - set boundary conditions according to configuration
+ */
+int initialize( algoparam_t *param , int *coords)
+{
+    int i, j;
+    double dist;
+
+    // total number of points 
+    const int np = param->act_res + 2;
+  
+    param->blocksize_x=param->act_res / param->thread_dims[0];
+    param->blocksize_y=param->act_res / param->thread_dims[1];
+
+	//
+    // allocate memory
+    //
+    (param->sendbuf_left) = (double*)calloc( sizeof(double), param->blocksize_y );
+    (param->recbuf_left) = (double*)calloc( sizeof(double), param->blocksize_y );
+    (param->sendbuf_right) = (double*)calloc( sizeof(double), param->blocksize_y );
+    (param->recbuf_right) = (double*)calloc( sizeof(double), param->blocksize_y );
+    (param->sendbuf_top) = (double*)calloc( sizeof(double), param->blocksize_x );
+    (param->recbuf_top) = (double*)calloc( sizeof(double), param->blocksize_x );
+    (param->sendbuf_bottom) = (double*)calloc( sizeof(double), param->blocksize_x );
+    (param->recbuf_bottom) = (double*)calloc( sizeof(double), param->blocksize_x );
+    (param->u)     = (double*)calloc( sizeof(double),(param->blocksize_x + 2)*(param->blocksize_y + 2) );
+    (param->uhelp) = (double*)calloc( sizeof(double),(param->blocksize_x + 2)*(param->blocksize_y + 2) );
+
+    if( !(param->u) || !(param->uhelp) || !(param->sendbuf_left) || !(param->sendbuf_right) || !(param->recbuf_left) || !(param->recbuf_right) )
+    {
+	fprintf(stderr, "Error: Cannot allocate memory\n");
+	return 0;
+    }
+
+    for( i=0; i<param->numsrcs; i++ )
+    {
+	/* top row */
+	if(coords[1]==0)
+	{
+	for( j=0 ; j<param->blocksize_x + 1 ; j++ )
+	{
+	    dist = sqrt( pow((double)(j + (param->blocksize_x * coords[0]) )/ (double)(np-1) - 
+			     param->heatsrcs[i].posx, 2)+
+			 pow(param->heatsrcs[i].posy, 2));
+	  
+	    if( dist <= param->heatsrcs[i].range )
+	    {
+		(param->u)[j] +=
+		    (param->heatsrcs[i].range-dist) /
+		    param->heatsrcs[i].range *
+		    param->heatsrcs[i].temp;
+		(param->uhelp)[j] = (param->u)[j];
+	    }
+	}
+ 	}	
+	/* bottom row */
+	if(coords[1]==param->thread_dims[1]-1)
+	{
+	for( j=0; j<param->blocksize_x + 1 ; j++ )
+	{
+
+            dist = sqrt( pow((double)(j + (param->blocksize_x * coords[0]) )/ (double)(np-1) -
+                             param->heatsrcs[i].posx, 2)+
+                         pow(1-param->heatsrcs[i].posy, 2));
+
+	    if( dist <= param->heatsrcs[i].range )
+	    {
+		(param->u)[(param->blocksize_y+1)*(param->blocksize_x+2) + j] +=
+		    (param->heatsrcs[i].range-dist) / 
+		    param->heatsrcs[i].range * 
+		    param->heatsrcs[i].temp;
+		(param->uhelp)[(param->blocksize_y+1)*(param->blocksize_x+2)+j] = (param->u)[(param->blocksize_y+1)*(param->blocksize_x+2)+j];
+	    }
+	}
+      }
+
+	if(coords[0]==0) 
+	{
+	/* leftmost column */
+	for( j=0; j<param->blocksize_y+1; j++ )
+	{
+	    dist = sqrt( pow(param->heatsrcs[i].posx, 2)+
+			 pow((double)(j + (param->blocksize_y * coords[1]) ) / (double)(np-1) - 
+			     param->heatsrcs[i].posy, 2)); 
+	  
+	    if( dist <= param->heatsrcs[i].range )
+	    {
+		(param->u)[ j*(param->blocksize_x+2) ]+=
+		    (param->heatsrcs[i].range-dist) / 
+		    param->heatsrcs[i].range *
+		    param->heatsrcs[i].temp;
+		(param->uhelp)[ j*(param->blocksize_x+2) ] = (param->u)[ j*(param->blocksize_x+2) ];
+	    }
+	}
+      }
+
+	if(coords[0]==param->thread_dims[0]-1)
+	{	
+	/* rightmost column */
+	for( j=1; j<param->blocksize_y+1; j++ )
+	{
+	 
+            dist = sqrt( pow(1 - param->heatsrcs[i].posx, 2)+
+                         pow((double)(j + (param->blocksize_y * coords[1]) ) / (double)(np-1) -  
+                             param->heatsrcs[i].posy, 2));
+ 
+	    if( dist <= param->heatsrcs[i].range )
+	    {
+		(param->u)[ j*(param->blocksize_x+2)+(param->blocksize_x+1) ]+=
+		    (param->heatsrcs[i].range-dist) /
+		    param->heatsrcs[i].range *
+		    param->heatsrcs[i].temp;
+		(param->uhelp)[ j*(param->blocksize_x+2)+(param->blocksize_x+1) ] = (param->u)[ j*(param->blocksize_x+2)+(param->blocksize_x+1) ];
+	    }
+	}
+	}
+    }
+    return 1;
+}
+
+/*
+ * free used memory
+ */
+int finalize( algoparam_t *param )
+{
+    if( param->u ) 
+	{
+		free(param->u);
+		param->u = 0;
+    }
+
+    if( param->uhelp ) 
+	{
+		free(param->uhelp);
+		param->uhelp = 0;
+    }
+
+	if ( param->sendbuf_left ) 
+	{
+		free(param->sendbuf_left);
+		param->sendbuf_left = 0;
+	}
+
+	if ( param->sendbuf_right ) 
+	{
+		free(param->sendbuf_right);
+		param->sendbuf_right = 0;
+	}
+
+	if ( param->recbuf_left ) 
+	{
+		free(param->recbuf_left);
+		param->recbuf_left = 0;
+	}
+
+	if ( param->recbuf_right ) 
+	{
+		free(param->recbuf_right);
+		param->recbuf_right = 0;
+	}
+
+	if ( param->sendbuf_top ) 
+	{
+		free(param->sendbuf_top);
+		param->sendbuf_top = 0;
+	}
+
+	if ( param->sendbuf_bottom ) 
+	{
+		free(param->sendbuf_bottom);
+		param->sendbuf_bottom = 0;
+	}
+
+	if ( param->recbuf_top ) 
+	{
+		free(param->recbuf_top);
+		param->recbuf_top = 0;
+	}
+
+	if ( param->recbuf_bottom ) 
+	{
+		free(param->recbuf_bottom);
+		param->recbuf_bottom = 0;
+	}
+
+    return 1;
+}
+

mpi-assign4/merged_code/pfmon2cg

+#!/usr/bin/perl
+#
+# pfmon2cg, version 0.02
+# Converter for pfmon output to Cachegrind profile format
+#
+# 4/2007, Josef Weidendorfer
+#
+
+#
+# Read in pfmon data, e.g. written from the command
+#
+#  pfmon --with-header --smpl-outfile=pfmon.out \
+#        -e l2_misses --short-smpl-periods=100000 -- <command>
+#  (sampling points every 100000 L2 Misses)
+#
+# This gives a list of sampled instruction addresses ($ipps)
+ 
+while(<>) {
+	if (/^#/) {
+		if (/command line\s*:\s*(.*)$/) { $command = $1; }
+		if (/PMD\d+: (.*?)\s*(,|=)/) { $events .= "$1 "; }
+		next; 
+	}
+	if (/^entry/) {
+		($iip) = ($entry =~ /IIP:(\S+)/);
+		if ($iip ne "") { $iips .= "$iip\n"; }
+		$entry = "";
+	}
+	chomp;
+	$entry .= $_;
+}
+($iip) = ($entry =~ /IIP:(\S+)/);
+if ($iip ne "") { $iips .= "$iip\n"; }
+
+($pid) = ($entry =~ /PID:(\S+)/);
+
+open OUT, ">.p2cg";
+print OUT $iips;
+close OUT;
+
+
+# Let addr2line convert the instruction addresses into
+# source code references (needs debug info)
+#
+# This creates hashes with event sums
+
+# Strip command line arguments
+$onlycommand = $command;
+$onlycommand =~ s/\s.*//;
+
+$cmd = "addr2line -f -e $onlycommand < .p2cg";
+#print "#Run Cmd $cmd\n";
+open ATOL, "$cmd < .p2cg|";
+while(<ATOL>) {
+	chomp;
+	$func = $_;
+	$fileline = <ATOL>;
+	chomp $fileline;
+	($file,$line) = ($fileline =~ /(.*?):(.*)/);
+	#print "Func $func, file $file, line $line\n";
+
+	if ($file{$func} eq "") {
+		push @funcs, $func;
+		$file{$func} = $file;
+	}
+	if ($count{$func,$line} eq "") {
+		$tmp = $lines->{$func};
+		$tmp = [] unless defined $tmp;
+		push @$tmp, $line;
+		$lines->{$func} = $tmp;
+	}
+	$count{$func,$line}++;
+	$totals++;
+}
+close ATOL;
+
+#
+# Dump out Cachegrind format
+#
+
+# First header info...
+print "version: 1\n";
+print "creator: pf2mon 0.02\n";
+print "pid: $pid\n";
+print "cmd: $command\n\n";
+
+print "positions: line\n";
+print "events: $events\n";
+print "summary: $totals\n\n";
+
+
+# Now profiling data.
+# 
+# We only support static linked binaries for now, so all code
+# is in the binary
+print "ob=$onlycommand\n";
+
+foreach $f (@funcs) {
+	print "\nfl=$file{$f}\n";
+	print "fn=$f\n";
+	$tmp = $lines->{$f};
+	foreach $l (@$tmp) {
+		print "$l $count{$f,$l}\n";
+	}
+}
+

mpi-assign4/merged_code/relax_jacobi.c

+/*
+ * relax_jacobi.c
+ *
+ * Jacobi Relaxation
+ *
+ */
+
+#include "heat.h"
+
+
+/*
+ * One Jacobi iteration step
+ */
+double relax_jacobi_return_residual( algoparam_t *param, int *coords, MPI_Comm comm2d )
+{
+	int i, j;
+	double tmp, diff, local_residual, global_residual;
+	int size_x, size_y;
+	int rank_left, rank_right, rank_top, rank_bottom;
+	MPI_Request sendrq_right, recrq_right, sendrq_left, recrq_left, sendrq_top, recrq_top, sendrq_bottom, recrq_bottom;
+	MPI_Status status;
+	
+	rank_left = get_rank_from_2dcoords(coords[0]-1, coords[1], param->thread_dims, comm2d);
+	rank_right = get_rank_from_2dcoords(coords[0]+1, coords[1], param->thread_dims, comm2d);
+	rank_top = get_rank_from_2dcoords(coords[0], coords[1]-1, param->thread_dims, comm2d);
+	rank_bottom = get_rank_from_2dcoords(coords[0], coords[1]+1, param->thread_dims, comm2d);
+	local_residual=0.0;
+	size_x = param->blocksize_x + 2;
+	size_y = param->blocksize_y + 2;
+
+	// do the calculation for all 4 borders
+		i= 1;
+		for( j=1; j<size_x-1; j++ )
+		{
+			param->uhelp[i*size_x + j]= 0.25 * (param->u[ i*size_x+j -1 ]+  // left
+						param->u[ i*size_x+j +1 ]+  // right
+						param->u[ (i-1)*size_x + j ]+  // top
+						param->u[ (i+1)*size_x + j ]); // bottom
+
+				diff = param->uhelp[i*size_x + j] - param->u[i*size_x +j];
+				local_residual += diff * diff;
+		}
+		j = 1;
+		for (i = 2; i<size_y-2; i++ )
+		{
+			param->uhelp[i*size_x + j]= 0.25 * (param->u[ i*size_x+j -1 ]+  // left
+						param->u[ i*size_x+j +1 ]+  // right
+						param->u[ (i-1)*size_x + j ]+  // top
+						param->u[ (i+1)*size_x + j ]); // bottom
+
+				diff = param->uhelp[i*size_x + j] - param->u[i*size_x +j];
+				local_residual += diff * diff;
+		}	
+		
+		i = size_y-2;
+		for( j=1; j<size_x-1; j++ )
+		{
+			param->uhelp[i*size_x + j]= 0.25 * (param->u[ i*size_x+j -1 ]+  // left
+						param->u[ i*size_x+j +1 ]+  // right
+						param->u[ (i-1)*size_x + j ]+  // top
+						param->u[ (i+1)*size_x + j ]); // bottom
+
+				diff = param->uhelp[i*size_x + j] - param->u[i*size_x +j];
+				local_residual += diff * diff;
+		}
+		j = size_x-2;
+		for (i = 2; i<size_y-2; i++ )
+		{
+			param->uhelp[i*size_x + j]= 0.25 * (param->u[ i*size_x+j -1 ]+  // left
+						param->u[ i*size_x+j +1 ]+  // right
+						param->u[ (i-1)*size_x + j ]+  // top
+						param->u[ (i+1)*size_x + j ]); // bottom
+
+				diff = param->uhelp[i*size_x + j] - param->u[i*size_x +j];
+				local_residual += diff * diff;
+		}
+
+
+	// Send border values to right neighbor
+		if (rank_right != -1)
+		{		
+			for (i=1; i < size_y-1; i++)
+				param->sendbuf_right[i-1] = param->uhelp[i*size_x + size_x-2];
+			MPI_Isend(param->sendbuf_right, size_y - 2, MPI_DOUBLE, rank_right, 10, comm2d, &sendrq_right);
+			
+		}
+	// Send border values to bottom neighbor
+		if (rank_bottom != -1)
+		{
+			for (j=1; j < size_x -1; j++)
+				param->sendbuf_bottom[j-1] = param->uhelp[(size_y-2)*size_x + j];
+			MPI_Isend(param->sendbuf_bottom, size_x - 2, MPI_DOUBLE, rank_bottom, 10, comm2d, &sendrq_bottom);
+		}
+	// Send border values to left neighbor
+		if (rank_left != -1)
+		{
+			for (i=1; i < size_y-1; i++)
+				param->sendbuf_left[i-1] = param->uhelp[i*size_x + 1];
+			MPI_Isend(param->sendbuf_left, size_y - 2, MPI_DOUBLE, rank_left, 10, comm2d, &sendrq_left);
+		}
+	// Send border values to top neighbor
+		if (rank_top != -1)
+		{
+			for (j=1; j < size_x-1; j++)
+				param->sendbuf_top[j-1] = param->uhelp[1*size_x + j];
+			MPI_Isend(param->sendbuf_top, size_x - 2, MPI_DOUBLE, rank_top, 10, comm2d, &sendrq_top);	
+		}
+
+	// Request receive from right neighbor
+		if (rank_right != -1)
+		{
+			MPI_Irecv(param->recbuf_right, size_y - 2, MPI_DOUBLE, rank_right, 10, comm2d, &recrq_right);
+			
+		}
+	// Request receive from bottom neighbor
+		if (rank_bottom != -1)
+		{						
+			MPI_Irecv(param->recbuf_bottom, size_x - 2, MPI_DOUBLE, rank_bottom, 10, comm2d, &recrq_bottom);
+		}
+	// Request receive from left neighbor
+		if (rank_left != -1)
+		{
+			MPI_Irecv(param->recbuf_left, size_y - 2, MPI_DOUBLE, rank_left, 10, comm2d, &recrq_left);
+			
+		}
+	// Request receive from top neighbor
+		if (rank_top != -1)
+		{						
+			MPI_Irecv(param->recbuf_top, size_x - 2, MPI_DOUBLE, rank_top, 10, comm2d, &recrq_top);
+		}
+	
+	// Calculate the inner part
+		for( i=2; i<size_y-2; i++ )
+		{
+			for( j=2; j<size_x-2; j++ )
+		
+			{
+				param->uhelp[i*size_x + j]= 0.25 * (param->u[ i*size_x+j -1 ]+  // left
+						param->u[ i*size_x+j +1 ]+  // right
+						param->u[ (i-1)*size_x + j ]+  // top
+						param->u[ (i+1)*size_x + j ]); // bottom
+
+				diff = param->uhelp[i*size_x + j] - param->u[i*size_x +j];
+				local_residual += diff * diff;
+			}
+		
+		}
+
+
+	// Wait for data from right neighbor
+		if (rank_right != -1)
+		{
+			MPI_Wait(&recrq_right, &status);
+			for (i=1; i < size_y-1; i++)
+				param->uhelp[i*size_x+ size_x-1] = param->recbuf_right[i-1];
+		}
+	// Wait for data from bottom neighbor
+		if (rank_bottom != -1)
+		{
+			MPI_Wait(&recrq_bottom, &status);
+			for (j=1; j < size_x-1; j++)
+				param->uhelp[(size_y-1)*size_x+ j] = param->recbuf_bottom[j-1];
+		}
+	// Wait for data from left neighbor
+		if (rank_left != -1)
+		{
+			MPI_Wait(&recrq_left, &status);
+			for (i=1; i < size_y-1; i++)
+				param->uhelp[i*size_x] = param->recbuf_left[i-1];
+		}
+	// Wait for data from top neighbor
+		if (rank_top != -1)
+		{
+			MPI_Wait(&recrq_top, &status);
+			for (j=1; j < size_x -1; j++)
+				param->uhelp[j] = param->recbuf_top[j-1];
+		}
+
+
+	
+	MPI_Reduce(&local_residual, &global_residual, 1, MPI_DOUBLE, MPI_SUM, 0, comm2d);
+	MPI_Bcast(&global_residual, 1, MPI_DOUBLE, 0, comm2d);
+	return global_residual;
+}

mpi-assign4/merged_code/test.dat

+2 2    # x and y dimensions for threads
+50     # iterations
+256    # initial resolution
+4096   # max resolution (spatial resolution)
+3840    # resolution step size
+1      # Algorithm 0=Jacobi 1=Gauss 
+2                     # number of heat sources
+0.0  0.0  1.0  1.0    # (x,y), size temperature
+1.0  1.0  1.0  0.5 

mpi-assign4/merged_code/timing.c

+//
+// timing.c 
+// 
+
+
+#include <sys/time.h>
+#include "timing.h"
+
+
+double wtime()
+{
+  struct timeval tv;
+  gettimeofday(&tv, 0);
+
+  return tv.tv_sec+1e-6*tv.tv_usec;
+}
+

mpi-assign4/merged_code/timing.h

+//
+// timing.h
+//
+
+#ifndef TIMING_H_INCLUDED
+#define TIMING_H_INCLUDED
+
+double wtime();
+
+#endif // TIMING_H_IINCLUDED