jacobi iteration with openmp

jacobi_openmp.c

/*****************************************************
* Exercise 4
* File: jacobi.c
*
* description:
*
*  Jacobi iteration on a NxN matrix with MAX_ITER iterations
*
* Language: c
*
* Compile:
*	cc -o bin/jacobi -O3 -xopenmp src/jacobi.c -lm
*
* Usage:
*	export OMP_NUM_THREADS=<num_threads>
*	./bin/jacobi
*
*	<num_threads> number of threads
*
*
* Martin Perdacher
* a1142622
*
*****************************************************/
#include <stdio.h>
#include <sys/time.h>
#include <omp.h>
#include <math.h>
#define N 1000
#define MAX_ITER 5000

typedef float MATRIX_V[N+1][N+1];

int timeval_subtract(struct timeval *result, struct timeval *t2, struct timeval *t1);

int main(int argc, char* argv) {
	int i,j,k;
	double diffnorm=0.0;
	int ncore= omp_get_num_procs(); // get max processors
	int nthreads,iter=0;
	MATRIX_V V, V_tmp;
    struct timeval start, end, diff;
	float *p_local, *p_local_tmp, *swap;
	
	
	// initialize all points with reasonable starting values for the potential
	for (i = 1; i < N; i++)
		for (j = 1; j < N; j++)
			V[i][j] = 51.0;
                
	// boundary conditions
	
	for (j = 0; j <= N; j++) V[0][j] =  30;   /// left/
	
	for (j = 0; j <= N; j++) V[N][j] =  70;   /// right/
	
	for (i = 0; i <= N; i++) V[i][0] =   0;   /// bottom/
	
	for (i = 0; i <= N; i++) V[i][N] = 100;   /// top/


	// copy whole array once
	for (i = 0 ; i <= N; i++)
   		for (j = 0 ; j <= N ; j++)
			V[i][j] = V_tmp[i][j];

	p_local=&V[0][0];
	p_local_tmp=&V_tmp[0][0];

	// get starting time
	gettimeofday(&start, NULL);
	

	#pragma omp parallel 
	{	
	
		nthreads = omp_get_num_threads();
		
		// numerically solve Laplace's equation
		do{
			diffnorm = 0.0;
			// calculate new values
			#pragma omp parallel for private(i,j) shared(V, V_tmp,p_local_tmp,p_local)
			for (i = 1; i < N; i++){
	      		for (j = 1; j < N; j++){
//					V_tmp[i][j] = 0.25* (V[i-1][j] + V[i+1][j] + V[i][j-1] + V[i][j+1]);
					*(p_local_tmp+i*N+j) = 0.25 * ( *(p_local+((i-1)*N)+j) + *(p_local+((i+1)*N)+j) + 
												    *(p_local+(i*N)+j-1) + *(p_local+(i*N)+j+1));
					
					diffnorm += ( *(p_local_tmp+i*N+j) - *(p_local+i*N+j) ) * 
								( *(p_local_tmp+i*N+j) - *(p_local+i*N+j) );				
				}
			}
		    
			// copy array
//			for (i = 1; i < N; i++)
//	      		for (j = 1; j < N; j++)
//							V[i][j] = V_tmp[i][j];
							
							
			swap = p_local_tmp; 
			p_local_tmp = p_local;
			p_local = swap;

			diffnorm = sqrt(diffnorm);
			iter++;
		}while( diffnorm > 1.0e-2 && iter < MAX_ITER );
	
	}

	// get end time
	gettimeofday(&end, NULL);
	timeval_subtract(&diff, &end, &start);

	printf("%d;%d;%d;",N,ncore,nthreads);
	printf("%ld.%06ld;%f;%d\n", diff.tv_sec, diff.tv_usec,diffnorm,iter);

}

int timeval_subtract(struct timeval *result, struct timeval *t2, struct timeval *t1)
{
    long int diff = (t2->tv_usec + 1000000 * t2->tv_sec) - (t1->tv_usec + 1000000 * t1->tv_sec);
    result->tv_sec = diff / 1000000;
    result->tv_usec = diff % 1000000;

    return (diff<0);
}