OpenMP Directives
Assignment
Report
Matrix Multiplication (mm)
Directives
I used a parallel do for the outer loop of the matrix initialization. There is no dependency since data is only written. I did not make the inner loop parallel because each thread should have some work to do, not only setting 2 array elements.
I also used a parallel do for the outer loop of the matrix multiplication. There is no data dependency, since data is only read from arrays a and b and only written to array c. Again, I only used the outer loop to put enough work into one thread and avoid excessive overhead through thread creation.
Speed-up Data
Speed-up Data Plots
Comments
Since calculation is very independent, each chunk of work is equally hard and the OpenMP overhead is quite small we are able to achieve a almost optimal speedup (see "550"-"ideal for 550" in first figure). So, I guess the point is that we can really simply get very good speedup results with OpenMP.
Source Code
mm.f: See end of report
Mandelbrot Calculation (mandel)
Directives
I used a parallel for construct for the outermost loop of the calculation, which is decomposition at row level. The reason was as before: Enough work for each thread. I think making the inner loop parallel would introduce too much overhead (thread creation).
Speed-up Data
Speed-up Data Plots
Comments
As can be seen from the figures the speedup is really bad. This is because the workload of each pixel is different and we are using a static decomposition. See section 5 (Extra Credit A) for a better solution.
Source Code
mandel.c: See end of report
2D Heat Equation (heat2d)
Directives
Again, I used a parallel for to parallize only one loop, for the same reasons as stated in the previous examples. But this time it was not the outermost loop, because it can't be parallized, as explained in addendum of Homework 3. Therefore I took the second but outermost loop.
Speed-up Data
Speed-up Data Plots
Comments
The problem in this example is that the outermost loop can not be parallized. We therefore pay the thread management overhead in every iteration of the outermost loop. Since the program takes quite a while to execute we can see what happens at small problem sizes.
Source Code
heat2d: See end of report
Poor Parallelization
Below are some test results to show how bad it is to have a large serial outer loop (iterations, fixed) and a small, but parallized, inner loop (problem size). As you can see, the speedup stays below 1 and even drops with the number of threads. Therefore don't parallize this kind of problem.
Comments and Conclusions
Parallezing code with OpenMP is really simple. One can concentrate on what to make parallel and leave the how to the compiler. PVM seems like an "assembler language," while OpenMP lifts parallelism to compiler level. But, as usual, you pay some overhead penalty for the convenience and portability.
When moving from a cluster to a shared memory architecture one must be aware that smaller decomposition sizes are possible and useful since the communication overhead is much smaller.
The ideal times given are calculated by serial_time/number_of_threads.
Extra Credit A (dynamic mandel)
Directives
I used the same directives as for the former test (see section 2, Mandelbrot Calculation (mandel)). In addition I specified a dynamic schedule with a chunk size of 5.
Speed-up Data
Speed-up Data Plots
Comments
As expected, the dynamic schedule balanced the workload for the threads. The main problem was to find a reasonable chunk size. I used a simple trial-and-error approach to find a good chunk size. The minimum problem size is 200 and we have at most 8 threads. So, starting with a chunk size of 25 (=200/8) seemed reasonable. I changed the chunk size from 25 to 20, 10 and finally 5, which gave a good result, ie almost optimal speedup (see first figures). So, the point is to balance the workload through dynamic decomposition.
Source Code
mandel_dyn.c: See end of report
Source Code
mm.f
* A program to multiply two matrices
*
* the program requests the user enter the dimensions of the matrices.
* it then allocates space for the arrays and initializes them.
* a timing function computes the elapsed time for the matrix multiply
*
program main
************************************************************
implicit none
integer n,m
double precision , dimension (:, :), ALLOCATABLE :: a,b,c
double precision sum
integer i,j,k
integer status
real gettime
integer, dimension(8):: tdat
double precision :: t1, t2, et
*
* Read dimension information
*
write(*,*) "Input array dimension n"
read(5,*) n
*
* allocate arrays dynamically
*
allocate (a(n, n), STAT = status)
if (status > 0) then
print*, "allocation error"
stop
endif
allocate (b(n, n), STAT = status)
if (status > 0) then
print*, "allocation error"
stop
endif
allocate (c(n, n), STAT = status)
if (status > 0) then
print*, "allocation error"
stop
endif
call initialize (n,a,b)
t1=gettime();
*
* multiply the matrices
*
!$omp parallel do private(i,j,sum,K) shared(a,b,c)
do i=1,n
do j = 1,n
sum = 0
do K=1,n
sum = sum + a(i,k)*b(k,j)
enddo
c(i,j) = sum
enddo
enddo
t2=gettime();
et = t2-t1
print*, "Dimension: ", n, " time: ",et
* do i=1,n
* do j=1,n
* print*, c(i,j)
* enddo
* enddo
stop
end
subroutine initialize (n,u,v)
******************************************************
* Initializes data
*
******************************************************
implicit none
integer n
double precision u(n,n),v(n,n)
integer omp_get_thread_num
integer i,j
* Initilize initial condition and RHS
!$omp parallel do private(j,i) shared(u,v)
do j = 1,n
do i = 1,n
u(i,j) = i+j
v(i,j) = i-j
enddo
enddo
return
end
real function gettime()
integer, dimension(8):: tdat
CALL DATE_AND_TIME(VALUES=tdat)
gettime=tdat(5)*3600 + tdat(6)*60+ tdat(7)+tdat(8)/1000.0
return
end
mandel.c
#include <stdio.h>
#include <string.h>
#include <math.h>
#define X_RESN 800 /* x resolution */
#define Y_RESN 800 /* y resolution */
#define X_MIN -2.0
#define X_MAX 2.0
#define Y_MIN -2.0
#define Y_MAX 2.0
typedef struct complextype
{
float real, imag;
} Compl;
int main ( int argc, char* argv[])
{
/* Mandlebrot variables */
int i, j, k;
Compl z, c;
float lengthsq, temp;
int maxIterations;
int res[X_RESN][Y_RESN];
/* Calculate and draw points */
maxIterations = atoi( argv[1] );
#pragma omp parallel for shared(res,maxIterations) private(i,j,z,c,k,temp,lengthsq)
for(i=0; i < Y_RESN; i++)
for(j=0; j < X_RESN; j++) {
z.real = z.imag = 0.0;
c.real = X_MIN + j * (X_MAX - X_MIN)/X_RESN;
c.imag = Y_MAX - i * (Y_MAX - Y_MIN)/Y_RESN;
k = 0;
do { /* iterate for pixel color */
temp = z.real*z.real - z.imag*z.imag + c.real;
z.imag = 2.0*z.real*z.imag + c.imag;
z.real = temp;
lengthsq = z.real*z.real+z.imag*z.imag;
k++;
} while (lengthsq < 4.0 && k < maxIterations);
if (k >= maxIterations) res[i][j] = 0;
else res[i][j] = 1;
}
printf("%d\n", res[0][0]);
/* Program Finished */
return 0;
}
heat2d.c
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#define CRESN 302 /* x, y resolution */
#define RESN 300
#define MAX_COLORS 8
#define MAX_HEAT 20.0
double solution[2][CRESN][CRESN], diff_constant;
int cur_gen, next_gen;
double sim_time, final = 10000.0, time_step = 0.1, diff = 10000.0;
void compute_one_iteration (double);
void setup ();
void
main (int argc, char *argv[])
{
int temp;
/* set window position */
/* Calculate and draw points */
setup ();
for (sim_time = 0; sim_time < final; sim_time += time_step)
{
compute_one_iteration (sim_time);
temp = cur_gen;
cur_gen = next_gen;
next_gen = temp;
}
}
void
setup ()
{
int i, j;
#pragma omp parallel for shared(solution) private(i,j)
for (i = 0; i < CRESN; i++)
for (j = 0; j < CRESN; j++)
solution[0][i][j] = solution[1][i][j] = 0.0;
cur_gen = 0;
next_gen = 1;
diff_constant = diff * time_step / ((float) RESN * (float) RESN);
}
void
compute_one_iteration (double t)
{
int i, j;
/* set boundary values */
for (i = 0; i < CRESN; i++)
{
if (i < 256 || i > 768)
solution[cur_gen][i][0] = solution[cur_gen][i][1];
else
solution[cur_gen][i][0] = MAX_HEAT;
}
for (i = 0; i < CRESN; i++)
{
solution[cur_gen][i][CRESN - 1] = solution[cur_gen][i][CRESN - 2];
}
for (i = 0; i < CRESN; i++)
{
solution[cur_gen][0][i] = solution[cur_gen][1][i];
solution[cur_gen][CRESN - 1][i] = solution[cur_gen][CRESN - 2][i];
}
/* corners ? */
#pragma omp parallel for shared(solution,cur_gen,next_gen,diff_constant) private(i,j)
for (i = 1; i <= RESN; i++)
for (j = 1; j <= RESN; j++)
solution[next_gen][i][j] = solution[cur_gen][i][j] +
(solution[cur_gen][i + 1][j] +
solution[cur_gen][i - 1][j] +
solution[cur_gen][i][j + 1] +
solution[cur_gen][i][j - 1] -
4.0 * solution[cur_gen][i][j]) *
diff_constant;
}
mandel_dyn.c
#include <stdio.h>
#include <string.h>
#include <math.h>
#define X_RESN 800 /* x resolution */
#define Y_RESN 800 /* y resolution */
#define X_MIN -2.0
#define X_MAX 2.0
#define Y_MIN -2.0
#define Y_MAX 2.0
typedef struct complextype
{
float real, imag;
} Compl;
int main ( int argc, char* argv[])
{
/* Mandlebrot variables */
int i, j, k;
Compl z, c;
float lengthsq, temp;
int maxIterations;
int res[X_RESN][Y_RESN];
/* Calculate and draw points */
maxIterations = atoi( argv[1] );
#pragma omp parallel for shared(res,maxIterations) private(i,j,z,c,k,temp,lengthsq) schedule (dynamic,5)
for(i=0; i < Y_RESN; i++)
for(j=0; j < X_RESN; j++) {
z.real = z.imag = 0.0;
c.real = X_MIN + j * (X_MAX - X_MIN)/X_RESN;
c.imag = Y_MAX - i * (Y_MAX - Y_MIN)/Y_RESN;
k = 0;
do { /* iterate for pixel color */
temp = z.real*z.real - z.imag*z.imag + c.real;
z.imag = 2.0*z.real*z.imag + c.imag;
z.real = temp;
lengthsq = z.real*z.real+z.imag*z.imag;
k++;
} while (lengthsq < 4.0 && k < maxIterations);
if (k >= maxIterations) res[i][j] = 0;
else res[i][j] = 1;
}
printf("%d\n", res[0][0]);
/* Program Finished */
return 0;
}
