Skip to content
Snippets Groups Projects
Commit 643698f9 authored by Mactavish's avatar Mactavish
Browse files

remove unused files

parent 4661ceac
Branches
No related tags found
No related merge requests found
template < typename T >
T **Allocate2DArray( int nRows, int nCols)
{
//(step 1) allocate memory for array of elements of column
T **ppi = (T **) malloc(sizeof(T *)*nRows);
//(step 2) allocate memory for array of elements of each row
T *curPtr = (T *) malloc(sizeof(T) * nRows * nCols);
// Now point the pointers in the right place
for( int i = 0; i < nRows; ++i)
{
*(ppi + i) = curPtr;
curPtr += nCols;
}
return ppi;
}
template < typename T >
void Free2DArray(T** Array)
{
free(*Array);
free(Array);
}
\ No newline at end of file
# copy to make.def
UNAME := $(shell uname)
CC = gcc
CFLAGS = -fopenmp
LIBS =
BINS = hello pi matmul pi_mc prod_cons \
matmul_recur mandel linked
BINS = hello pi #pi_spmd_simple_soln pi_spmd_final_soln pi_loop_soln
#
ifeq ($(UNAME), Darwin)
CFLAGS = -Xpreprocessor -fopenmp
LIBS += -lomp
......@@ -21,18 +18,9 @@ all: $(BINS)
$(BINS): %: %.o
$(CC) $(CFLAGS) -o $@ $@.o $(LIBS)
pi_mc: pi_mc.o random.o
$(CC) $(CFLAGS) -o $@ $^ $(LIBS)
test: $(BINS)
./hello
./pi
./matmul
./pi_mc
./prod_cons
./matmul_recur
./mandel
./linked
clean:
$(RM) $(BINS) *.o
# OpenMP
Here are some small exercises that help you understand the OpenMP basics.
More details in the slides and in the tutorial.
## Build
```bash
# build all
make
# test all executables
make test
# clean all compilation files
make clean
```
## Credit
The exercises are based on [this tutorial](https://github.com/tgmattso/OpenMP_intro_tutorial)
#include <stdlib.h>
#include <stdio.h>
#include <omp.h>
#ifndef N
#define N 5
#endif
#ifndef FS
#define FS 38
#endif
struct node {
int data;
int fibdata;
struct node* next;
};
int fib(int n) {
int x, y;
if (n < 2) {
return (n);
} else {
x = fib(n - 1);
y = fib(n - 2);
return (x + y);
}
}
void processwork(struct node* p)
{
int n;
n = p->data;
p->fibdata = fib(n);
}
struct node* init_list(struct node* p) {
int i;
struct node* head = NULL;
struct node* temp = NULL;
head = malloc(sizeof(struct node));
p = head;
p->data = FS;
p->fibdata = 0;
for (i=0; i< N; i++) {
temp = malloc(sizeof(struct node));
p->next = temp;
p = temp;
p->data = FS + i + 1;
p->fibdata = i+1;
}
p->next = NULL;
return head;
}
int main(int argc, char *argv[]) {
double start, end;
struct node *p=NULL;
struct node *temp=NULL;
struct node *head=NULL;
printf("Process linked list\n");
printf(" Each linked list node will be processed by function 'processwork()'\n");
printf(" Each ll node will compute %d fibonacci numbers beginning with %d\n",N,FS);
p = init_list(p);
head = p;
start = omp_get_wtime();
{
while (p != NULL) {
processwork(p);
p = p->next;
}
}
end = omp_get_wtime();
p = head;
while (p != NULL) {
printf("%d : %d\n",p->data, p->fibdata);
temp = p->next;
free (p);
p = temp;
}
free (p);
printf("Compute Time: %f seconds\n", end - start);
return 0;
}
/*
** PROGRAM: Mandelbrot area
**
** PURPOSE: Program to compute the area of a Mandelbrot set.
** Correct answer should be around 1.510659.
** WARNING: this program may contain errors
**
** USAGE: Program runs without input ... just run the executable
**
** HISTORY: Written: (Mark Bull, August 2011).
** Changed "comples" to "d_comples" to avoid collsion with
** math.h complex type (Tim Mattson, September 2011)
*/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <omp.h>
# define NPOINTS 1000
# define MAXITER 1000
void testpoint(void);
struct d_complex{
double r;
double i;
};
struct d_complex c;
int numoutside = 0;
int main(){
int i, j;
double area, error, eps = 1.0e-5;
// Loop over grid of points in the complex plane which contains the Mandelbrot set,
// testing each point to see whether it is inside or outside the set.
#pragma omp parallel for default(shared) private(c,eps)
for (i=0; i<NPOINTS; i++) {
for (j=0; j<NPOINTS; j++) {
c.r = -2.0+2.5*(double)(i)/(double)(NPOINTS)+eps;
c.i = 1.125*(double)(j)/(double)(NPOINTS)+eps;
testpoint();
}
}
// Calculate area of set and error estimate and output the results
area=2.0*2.5*1.125*(double)(NPOINTS*NPOINTS-numoutside)/(double)(NPOINTS*NPOINTS);
error=area/(double)NPOINTS;
printf("Area of Mandlebrot set = %12.8f +/- %12.8f\n",area,error);
printf("Correct answer should be around 1.510659\n");
}
void testpoint(void){
// Does the iteration z=z*z+c, until |z| > 2 when point is known to be outside set
// If loop count reaches MAXITER, point is considered to be inside the set
struct d_complex z;
int iter;
double temp;
z=c;
for (iter=0; iter<MAXITER; iter++){
temp = (z.r*z.r)-(z.i*z.i)+c.r;
z.i = z.r*z.i*2+c.i;
z.r = temp;
if ((z.r*z.r+z.i*z.i)>4.0) {
numoutside++;
break;
}
}
}
/*
** PROGRAM: Matrix Multiply
**
** PURPOSE: This is a simple matrix multiply program.
** It will compute the product
**
** C = A * B
**
** A and B are set to constant matrices so we
** can make a quick test of the multiplication.
**
** USAGE: Right now, I hardwire the martix dimensions.
** later, I'll take them from the command line.
**
** HISTORY: Written by Tim Mattson, Nov 1999.
*/
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#define ORDER 1000
#define AVAL 3.0
#define BVAL 5.0
#define TOL 0.001
int main(int argc, char **argv)
{
int Ndim, Pdim, Mdim; /* A[N][P], B[P][M], C[N][M] */
int i,j,k;
double *A, *B, *C, cval, tmp, err, errsq;
double dN, mflops;
double start_time, run_time;
Ndim = ORDER;
Pdim = ORDER;
Mdim = ORDER;
A = (double *)malloc(Ndim*Pdim*sizeof(double));
B = (double *)malloc(Pdim*Mdim*sizeof(double));
C = (double *)malloc(Ndim*Mdim*sizeof(double));
/* Initialize matrices */
for (i=0; i<Ndim; i++)
for (j=0; j<Pdim; j++)
*(A+(i*Ndim+j)) = AVAL;
for (i=0; i<Pdim; i++)
for (j=0; j<Mdim; j++)
*(B+(i*Pdim+j)) = BVAL;
for (i=0; i<Ndim; i++)
for (j=0; j<Mdim; j++)
*(C+(i*Ndim+j)) = 0.0;
/* Do the matrix product */
start_time = omp_get_wtime();
for (i=0; i<Ndim; i++){
for (j=0; j<Mdim; j++){
tmp = 0.0;
for(k=0;k<Pdim;k++){
/* C(i,j) = sum(over k) A(i,k) * B(k,j) */
tmp += *(A+(i*Ndim+k)) * *(B+(k*Pdim+j));
}
*(C+(i*Ndim+j)) = tmp;
}
}
/* Check the answer */
run_time = omp_get_wtime() - start_time;
printf(" Order %d multiplication in %f seconds \n", ORDER, run_time);
dN = (double)ORDER;
mflops = 2.0 * dN * dN * dN/(1000000.0* run_time);
printf(" Order %d multiplication at %f mflops\n", ORDER, mflops);
cval = Pdim * AVAL * BVAL;
errsq = 0.0;
for (i=0; i<Ndim; i++){
for (j=0; j<Mdim; j++){
err = *(C+i*Ndim+j) - cval;
errsq += err * err;
}
}
if (errsq > TOL)
printf("\n Errors in multiplication: %f",errsq);
else
printf("\n Hey, it worked");
printf("\n all done \n");
}
// Several versions of serial codes for matrix-matrix multiplication
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include "2DArray.h"
// define sizes of matrices to be used
#define MM 1000
#define NN 1000
#define PP 1000
double dabs(double d){return (d<0.0?d:(-d));}
// Default triple-nested loop for matrix-matrix multiplication
void matmult1(int m, int n, int p, double **A, double **B, double **C)
{
int i, j, k;
for (i = 0; i < m; i++)
for (j = 0; j < n; j++){
C[i][j]=0;
for (k = 0; k < p; k++)
C[i][j] += A[i][k]*B[k][j];
}
}
/*
Recursive code for matrix multiplication.
The recursion uses the formula
C00 = A00*B00 + A01*B10
C01 = A00*B01 + B01*B11
C10 = A10*B00 + A11*B10
C11 = A10*B01 + A11*B11
*/
void matmultleaf(int mf, int ml, int nf, int nl, int pf, int pl, double **A, double **B, double **C)
/*
subroutine that uses the simple triple loop to multiply
a submatrix from A with a submatrix from B and store the
result in a submatrix of C.
(We could use a tiled version,for better performance)
*/
// mf, ml; /* first and last+1 i index */
// nf, nl; /* first and last+1 j index */
// pf, pl; /* first and last+1 k index */
{
int i,j,k;
for (i = mf; i < ml; i++)
for (j = nf; j < nl; j++)
for (k = pf; k < pl; k++)
C[i][j] += A[i][k]*B[k][j];
}
#define GRAIN 32768 /* product size below which matmultleaf is used */
void matmultrec(int mf, int ml, int nf, int nl, int pf, int pl, double **A, double **B, double **C)
/*
recursive subroutine to compute the product of two
submatrices of A and B and store the result in C
*/
// mf, ml; /* first and last+1 i index */
// nf, nl; /* first and last+1 j index */
// pf, pl; /* first and last+1 k index */
{
//
// Check sizes of matrices;
// if below threshold then compute product w/o recursion
//
if ((ml-mf)*(nl-nf)*(pl-pf) < GRAIN)
matmultleaf(mf, ml, nf, nl, pf, pl, A, B, C);
else {
//
// Apply OpenMP tasks to the eight recursive calls below
// be sure to not create data races between tasks
//
// C00 += A00 * B00
matmultrec(mf, mf+(ml-mf)/2, nf, nf+(nl-nf)/2, pf, pf+(pl-pf)/2, A, B, C);
// C01 += A00 * B01
matmultrec(mf, mf+(ml-mf)/2, nf+(nl-nf)/2, nl, pf, pf+(pl-pf)/2, A, B, C);
// C00 += A01 * B10
matmultrec(mf, mf+(ml-mf)/2, nf, nf+(nl-nf)/2, pf+(pl-pf)/2, pl, A, B, C);
// C01 += A01 * B11
matmultrec(mf, mf+(ml-mf)/2, nf+(nl-nf)/2, nl, pf+(pl-pf)/2, pl, A, B, C);
// C10 += A10 * B00
matmultrec(mf+(ml-mf)/2, ml, nf, nf+(nl-nf)/2, pf, pf+(pl-pf)/2, A, B, C);
// C11 += A10 * B01
matmultrec(mf+(ml-mf)/2, ml, nf+(nl-nf)/2, nl, pf, pf+(pl-pf)/2, A, B, C);
// C10 += A11 * B10
matmultrec(mf+(ml-mf)/2, ml, nf, nf+(nl-nf)/2, pf+(pl-pf)/2, pl, A, B, C);
// C11 += A11 * B11
matmultrec(mf+(ml-mf)/2, ml, nf+(nl-nf)/2, nl, pf+(pl-pf)/2, pl, A, B, C);
}
}
//
// "Helper" function to intialize C and start recursive routine
//
void matmultr(int m, int n, int p, double **A, double **B, double **C)
{
int i,j;
for (i = 0; i < m; i++)
for (j=0; j < n; j++)
C[i][j] = 0;
matmultrec(0, m, 0, n, 0, p, A, B, C);
}
int CheckResults(int m, int n, double **C, double **C1)
{
#define ERR_THRESHOLD 0.001
int code = 0;
//
// May need to take into consideration the floating point roundoff error
// due to parallel execution
//
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
if (dabs(C[i][j] - C1[i][j]) > ERR_THRESHOLD ) {
printf("%f %f at [%d][%d]\n", C[i][j], C1[i][j], i, j);
code = 1;
}
}
}
return code;
}
int main(int argc, char* argv[])
{
int i, j;
double start, time1, time2;
int M = MM;
int N = NN;
int P = PP;
//
// If 3 values on command line, use those for matrix sizes
//
if (argc != 4) {
printf("Suggested Usage: %s <M> <N> <P> \n", argv[0]);
printf("Using default values\n");
}
else {
M = atoi(argv[1]);
N = atoi(argv[2]);
P = atoi(argv[3]);
}
double **A = Allocate2DArray< double >(M, P);
double **B = Allocate2DArray< double >(P, N);
double **C1 = Allocate2DArray< double >(M, N);
double **C4 = Allocate2DArray< double >(M, N);
//
// Initialize with random values
//
for (i = 0; i < M; i++) {
for (j = 0; j < P; j++) {
A[i][j] = (double)(rand()%100) / 10.0;
}
}
for (i = 0; i < P; i++) {
for (j = 0; j < N; j++) {
B[i][j] = (double)(rand()%100) / 10.0;
}
}
printf("Matrix Dimensions: M = %d P = %d N = %d\n\n", M, P, N);
printf("Execute matmult1\n");
start = omp_get_wtime();
matmult1(M, N, P, A, B, C1);
time1 = omp_get_wtime() - start;
printf("Time = %f seconds\n\n",time1);
printf("Execute matmultr\n");
start = omp_get_wtime();
matmultr(M, N, P, A, B, C4);
time2 = omp_get_wtime() - start;
printf("Time = %f seconds\n\n",time2);
printf("Checking...");
if (CheckResults(M, N, C1, C4))
printf("Error in Recursive Matrix Multiplication\n\n");
else {
printf("OKAY\n\n");
printf("Speedup = %5.1fX\n", time1/time2);
}
Free2DArray< double >(A);
Free2DArray< double >(B);
Free2DArray< double >(C1);
Free2DArray< double >(C4);
return 0;
}
......@@ -20,6 +20,7 @@ static long num_steps = 100000000;
double step;
int main() {
// TODO: parallelize this code using OpenMP
int i;
double x, pi, sum = 0.0;
double start_time, run_time;
......
/*
NAME:
Pi_mc: PI Monte Carlo
Purpose:
This program uses a Monte Carlo algorithm to compute PI as an
example of how random number generators are used to solve problems.
Note that if your goal is to find digits of pi, there are much
better algorithms you could use.
Usage:
To keep the program as simple as possible, you must edit the file
and change the value of num_trials to change the number of samples
used. Then compile and run the program.
Algorithm:
The basic idea behind the algorithm is easy to visualize. Draw a
square on a wall. Inside the square, draw a circle. Now randomly throw
darts at the wall. some darts will land inside the square. Of those,
some will fall inside the circle. The probability of landing inside
the circle or the square is proportional to their areas.
We can use a random number generator to "throw the darts" and count
how many "darts" fall inside the square and how many inside the
cicle. Dividing these two numbers gives us the ratio of their areas
and from that we can compute pi.
Algorithm details:
To turn this into code, I need a bit more detail. Assume the circle
is centered inside the square. the circle will have a radius of r and
each side of the square will be of area 2*r (i.e. the diameter of the
circle).
A(circle) = pi * r^2
A(square) = (2*r)*(2*r) = 4*r^2
ratio = A(circle)/A(square) = pi/4
Since the probability (P) of a dart falling inside a figure (i.e. the square
or the circle) is proportional to the area, we have
ratio = P(circle)/P(square) = pi/4
If I throw N darts as computed by random numbers evenly distributed
over the area of the square
P(sqaure) = N/N .... i.e. every dart lands in the square
P(circle) = N(circle)/N
ratio = (N(circle)/N)/(N/N) = N(circle)/N
Hence, to find the area, I compute N random "darts" and count how many fall
inside the circle. The equation for a circle is
x^2 + y^2 = r^2
So I randomly compute "x" and "y" evenly distributed from -r to r and
count the "dart" as falling inside the cicle if
x^2 + y^2 < or = r
Results:
Remember, our goal is to demonstrate a simple monte carlo algorithm,
not compute pi. But just for the record, here are some results (Intel compiler
version 10.0, Windows XP, core duo laptop)
100 3.160000
1000 3.148000
10000 3.154000
100000 3.139920
1000000 3.141456
10000000 3.141590
100000000 3.141581
As a point of reference, the first 7 digits of the true value of pi
is 3.141592
History:
Written by Tim Mattson, 9/2007.
*/
#include <stdio.h>
#include <omp.h>
#include "random.h"
//
// The monte carlo pi program
//
static long num_trials = 10000;
int main ()
{
long i; long Ncirc = 0;
double pi, x, y, test;
double r = 1.0; // radius of circle. Side of squrare is 2*r
seed(-r, r); // The circle and square are centered at the origin
#pragma omp parallel for private(x,y,test) reduction(+:Ncirc)
for(i=0;i<num_trials; i++)
{
x = drandom();
y = drandom();
test = x*x + y*y;
if (test <= r*r) Ncirc++;
}
pi = 4.0 * ((double)Ncirc/(double)num_trials);
printf("\n %ld trials, pi is %lf \n",num_trials, pi);
return 0;
}
/*
** PROGRAM: A simple serial producer/consumer program
**
** One function generates (i.e. produces) an array of random values.
** A second functions consumes that array and sums it.
**
** HISTORY: Written by Tim Mattson, April 2007.
*/
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#define N 10000
/* Some random number constants from numerical recipies */
#define SEED 2531
#define RAND_MULT 1366
#define RAND_ADD 150889
#define RAND_MOD 714025
int randy = SEED;
/* function to fill an array with random numbers */
void fill_rand(int length, double *a)
{
int i;
for (i=0;i<length;i++) {
randy = (RAND_MULT * randy + RAND_ADD) % RAND_MOD;
*(a+i) = ((double) randy)/((double) RAND_MOD);
}
}
/* function to sum the elements of an array */
double Sum_array(int length, double *a)
{
int i; double sum = 0.0;
for (i=0;i<length;i++) sum += *(a+i);
return sum;
}
int main()
{
double *A, sum, runtime;
int flag = 0;
A = (double *)malloc(N*sizeof(double));
runtime = omp_get_wtime();
fill_rand(N, A); // Producer: fill an array of data
sum = Sum_array(N, A); // Consumer: sum the array
runtime = omp_get_wtime() - runtime;
printf(" In %f seconds, The sum is %f \n",runtime,sum);
}
//**********************************************************
// Pseudo random number generator:
// double random
// void seed (lower_limit, higher_limit)
//**********************************************************
//
// A simple linear congruential random number generator
// (Numerical Recipies chapter 7, 1st ed.) with parameters
// from the table on page 198j.
//
// Uses a linear congruential generator to return a value between
// 0 and 1, then scales and shifts it to fill the desired range. This
// range is set when the random number generator seed is called.
//
// USAGE:
//
// pseudo random sequence is seeded with a range
//
// void seed(lower_limit, higher_limit)
//
// and then subsequent calls to the random number generator generates values
// in the sequence:
//
// double drandom()
//
// History:
// Written by Tim Mattson, 9/2007.
// changed to drandom() to avoid collision with standard libraries, 11/2011
static long MULTIPLIER = 1366;
static long ADDEND = 150889;
static long PMOD = 714025;
long random_last = 0;
double random_low, random_hi;
double drandom()
{
long random_next;
double ret_val;
//
// compute an integer random number from zero to mod
//
random_next = (MULTIPLIER * random_last + ADDEND)% PMOD;
random_last = random_next;
//
// shift into preset range
//
ret_val = ((double)random_next/(double)PMOD)*(random_hi-random_low)+random_low;
return ret_val;
}
//
// set the seed and the range
//
void seed(double low_in, double hi_in)
{
if(low_in < hi_in)
{
random_low = low_in;
random_hi = hi_in;
}
else
{
random_low = hi_in;
random_hi = low_in;
}
random_last = PMOD/ADDEND; // just pick something
}
//**********************************************************
// end of pseudo random generator code.
//**********************************************************
double drandom();
void seed(double low_in, double hi_in);
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment