Commit efeb86fd authored by Thomas Haener's avatar Thomas Haener
Browse files

added solution03

parent 5558e25b
all: simple_sampling_serial simple_sampling_omp
simple_sampling_serial: simple_sampling_serial.cpp
g++ -Wall -O3 --std=c++11 -o simple_sampling_serial simple_sampling_serial.cpp
simple_sampling_omp: simple_sampling_omp.cpp
g++ -Wall -O3 --std=c++11 -fopenmp -o simple_sampling_omp simple_sampling_omp.cpp
clean:
rm -f simple_sampling_serial simple_sampling_omp
2 0.630109 0.370499
4 0.528984 0.213223
8 0.342852 0.131797
16 0.24691 0.0708502
32 0.290587 0.0507888
64 0.327917 0.0438269
128 0.27963 0.0278127
256 0.314442 0.0219146
512 0.289087 0.0151105
1024 0.303712 0.010204
2048 0.301339 0.007341
4096 0.302506 0.00520078
8192 0.302555 0.00371337
16384 0.295586 0.00259773
32768 0.300143 0.00184166
from math import pi
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
def isfloat(x):
try:
float(x)
return True
except:
return False
data = np.loadtxt('integral.dat', unpack=True)
x = data[0]
y = data[1]
y_errors = data[2]
# === approximations
plt.xscale('log')
plt.yscale('linear')
current_plot ,= plt.plot(x,y,'-o',markersize=4.0)
plt.errorbar(x,y,yerr=y_errors,fmt='o',markersize=4.0,color=current_plot.get_color())
plt.axhline(y = 0.3, linestyle='dashed', color='k')
plt.ylabel(r'$I$',fontsize=16)
plt.xlabel(r'$N$',fontsize=16)
plt.savefig('integral.pdf',dpi=600)
plt.show()
# === error convergence
plt.xscale('log')
plt.yscale('log')
current_plot ,= plt.plot(x,y_errors,'-o',markersize=4.0)
plt.ylabel(r'$\Delta I$',fontsize=16)
plt.xlabel(r'$N$',fontsize=16)
plt.savefig('integral_errors.pdf',dpi=600)
plt.show()
# === scaling
#directly parse file
x = []
y = []
error= []
with open("scaling.dat") as f:
for line in f:
line=line.strip().split()
if(len(line)!=0):
if(line[0]=="N:"):
x.append([])
y.append([])
error.append([])
elif(isfloat(line[0])):
x[-1].append(int(line[0]))
y[-1].append(float(line[1]))
error[-1].append(float(line[2]))
plt.xscale('linear')
plt.yscale('linear')
##change color
cmap = mpl.cm.rainbow
#optimal scaling
optimal = np.arange(1,25,1)
plt.plot(optimal,optimal,label='optimal', color="black")
#now plot each line
for i in range(0,len(y),2):
y[i] = y[i][0] /np.array(y[i])
print(y[i])
plt.plot(x[i],y[i],label='2^%i'%(i+1), color=cmap(i/float(len(y))))
plt.legend (loc = 'best')
xmin, xmax = plt.xlim()
plt.ylim( (xmin, xmax) )
plt.ylabel(r'$T_1/T_N$',fontsize=16)
plt.xlabel(r'#threads',fontsize=16)
plt.savefig('scaling.pdf',dpi=600)
plt.show()
# runtimes
********************
N: 1
------------
1 0.000226294 2.73905e-06
2 0.000228364 3.78358e-06
4 0.000228328 3.15366e-06
8 0.000258876 7.37445e-06
16 0.000272379 1.1455e-05
24 0.00302437 7.44037e-05
********************
N: 2
------------
1 0.00037433 9.89554e-08
2 0.000202858 6.78554e-06
4 0.000202452 5.75721e-06
8 0.000213709 3.48326e-06
16 0.000223005 4.2986e-06
24 0.00101741 0.000286535
********************
N: 4
------------
1 0.000614641 2.58053e-06
2 0.000305689 6.6777e-06
4 0.000202504 8.22366e-06
8 0.000210961 3.67558e-06
16 0.000216846 3.55935e-06
24 0.000871304 0.000267495
********************
N: 8
------------
1 0.00117046 1.25358e-05
2 0.000736727 7.57006e-06
4 0.000501195 7.21369e-06
8 0.000351494 4.03051e-06
16 0.000364358 4.38023e-06
24 0.00036692 5.32095e-06
********************
N: 16
------------
1 0.00210755 1.45065e-05
2 0.00125737 1.73688e-05
4 0.000677491 9.04715e-06
8 0.000497466 5.16312e-06
16 0.000447241 8.31082e-06
24 0.00162532 0.000334028
********************
N: 32
------------
1 0.00423347 1.56254e-05
2 0.00226469 2.52928e-05
4 0.00131854 9.01054e-06
8 0.000857339 5.42914e-06
16 0.000571035 8.35197e-06
24 0.000951861 0.000250897
********************
N: 64
------------
1 0.0102499 1.66829e-05
2 0.00575999 2.11096e-05
4 0.00295561 6.81385e-06
8 0.00217856 6.23578e-06
16 0.00190438 1.75784e-05
24 0.00171414 8.39708e-05
********************
N: 128
------------
1 0.022585 1.92096e-05
2 0.0119533 2.84173e-05
4 0.00738901 1.05856e-05
8 0.00405907 3.65872e-06
16 0.00278451 1.00695e-05
24 0.00491851 0.00024345
********************
N: 256
------------
1 0.0432157 2.32481e-05
2 0.0238573 5.76296e-05
4 0.0125464 1.23535e-05
8 0.00732915 8.13726e-06
16 0.00439049 8.04745e-06
24 0.00377393 0.000367446
********************
N: 512
------------
1 0.0845487 1.93963e-05
2 0.043503 8.13873e-05
4 0.0223268 1.78677e-05
8 0.0143931 1.68831e-05
16 0.0081564 9.10133e-06
24 0.00623385 7.56613e-05
********************
N: 1024
------------
1 0.181464 3.91246e-05
2 0.0926455 7.38742e-05
4 0.0472356 2.43914e-05
8 0.0278302 1.87887e-05
16 0.0145261 1.14959e-05
24 0.0110035 0.000544944
********************
N: 2048
------------
1 0.355321 5.4037e-05
2 0.181339 0.000150155
4 0.0940211 4.0211e-05
8 0.0520544 1.73453e-05
16 0.0272996 1.29613e-05
24 0.0200894 0.000521937
********************
N: 4096
------------
1 0.704936 4.44657e-05
2 0.353712 0.000260601
4 0.182444 4.75366e-05
8 0.100196 3.41458e-05
16 0.0546285 1.59652e-05
24 0.037717 0.00065204
********************
N: 8192
------------
1 1.38353 7.95007e-05
2 0.690551 0.000172427
4 0.351261 0.000260311
8 0.19538 2.76196e-05
16 0.103456 1.75067e-05
24 0.0695684 1.44762e-05
********************
N: 16384
------------
1 2.76779 0.000135056
2 1.38361 0.000312325
4 0.701693 0.000272662
8 0.385656 9.31814e-05
16 0.202718 1.46731e-05
24 0.152669 0.00624627
********************
N: 32768
------------
1 5.57993 0.000238185
2 2.79908 0.000102401
4 1.40716 0.000354378
8 0.792049 4.8954e-05
16 0.408063 1.51241e-05
24 0.276663 1.63749e-05
********************
N: 65536
------------
1 11.1809 0.000381863
2 5.58191 0.000297387
4 2.83357 0.000804982
8 1.55447 0.000399105
16 0.803619 2.31566e-05
24 0.597334 0.0216178
********************
N: 131072
------------
1 22.347 0.000685921
2 11.1601 0.000384305
4 5.63894 0.00186011
8 3.09329 0.000601403
16 1.61003 2.12438e-05
24 1.30639 0.0541003
#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>
#include <random>
#include <cassert>
#include <iostream>
#include <iterator>
#include <chrono>
#include <omp.h>
template <class BC>
class RandomWalk {
typedef std::mt19937 generator_type;
public:
RandomWalk ( const std::vector<double> ipos_, const double d_ )
: ipos (ipos_)
, d (d_)
{ }
double operator () (const std::size_t seed) {
std::vector<double> pos (ipos);
generator_type generator (seed);
std::function <double()> rng = std::bind ( std::uniform_real_distribution <> (0, 1), std::ref (generator) );
double direction;
while ( bc.inside (pos) ) {
direction = rng ();
pos [0] += d * std::cos (2 * M_PI * direction);
pos [1] += d * std::sin (2 * M_PI * direction);
}
return bc.evaluate ( pos );
}
private:
BC bc;
const std::vector<double> ipos;
double d;
};
template <class T>
class MC {
public:
MC (const T& func_, std::size_t seed_)
: func (func_)
, seed (seed_)
{}
void sample (const std::size_t N, const unsigned num_threads) {
double sum_ = 0.0, sum2_ = 0.0;
#pragma omp parallel num_threads(num_threads)
#pragma omp for reduction(+:sum_) reduction(+:sum2_) schedule(static)
//#pragma omp for reduction(+:sum_) reduction(+:sum2_) schedule(guided,4)
for (std::size_t i = 0; i < N; ++i) {
double value = func (seed * N + i);
sum_ += value;
sum2_ += value * value;
}
sum = sum_;
sum2 = sum2_;
}
double mean (const std::size_t N) const
{
return sum / double (N);
}
double error (const std::size_t N) const
{
return std::sqrt ( ( sum2 - sum * sum / double(N) ) / double ( N * (N-1) ) );
}
private:
const std::function <double(std::size_t)> func;
const std::size_t seed;
double sum, sum2;
};
struct BC {
const double a = 1;
const double b = 1;
bool inside (const std::vector<double> &pos) const {
assert (pos.size() == 2);
return pos[0] > 0 && pos[0] < a && pos [1] > 0 && pos [1] < b;
}
virtual double evaluate (const std::vector<double>& pos) const { return 0; }
double project (const std::vector<double>& pos_) {
std::vector<double> pos (pos_);
if ( pos [0] < 0 ) pos [0] = 0;
if ( pos [1] < 0 ) pos [1] = 0;
if ( pos [0] > a ) pos [0] = a;
if ( pos [1] > b ) pos [1] = b;
return evaluate (pos);
}
};
struct BC_Dirichlet1 : BC {
double evaluate (const std::vector<double>& pos) const {
assert (pos.size() == 2);
return pos[0];
}
};
int main()
{
const double d = 0.01;
const std::vector<double> ipos = {0.3, 0.4};
const std::size_t reps = 18;
const std::size_t seed = 17;
RandomWalk < BC_Dirichlet1 > rw (ipos, d);
MC < RandomWalk <BC_Dirichlet1> > mc (rw, seed);
const std::vector<unsigned> threads_vec = {1,2,4,8,16,24};
std::chrono::time_point<std::chrono::high_resolution_clock> start,end;
std::cout << "# runtimes" << std::endl;
for (std::size_t i = 0; i < reps; i++) {
const std::size_t N = std::size_t(1) << i;
std::cout << "********************" << std::endl;
std::cout << "N: " << N << std::endl;
std::cout << "------------" << std::endl;
for (const auto num_threads : threads_vec) {
std::vector<double> runtimes;
for (unsigned j = 0; j < 20; ++j) {
start = std::chrono::high_resolution_clock::now();
mc.sample (N, num_threads);
end = std::chrono::high_resolution_clock::now();
const double elapsed = std::chrono::duration<double>(end-start).count();
runtimes.push_back(elapsed);
}
//std::sort (runtimes.begin(), runtimes.end());
//const double time_min = *std::min_element (runtimes.begin(), runtimes.end());
//const double time_median = runtimes[5];
//mean and error
const double time_mean = std::accumulate(runtimes.begin(),runtimes.end(),(double) 0.0)/runtimes.size();
double time_err = 0.0;
std::for_each(runtimes.begin(),runtimes.end(), [&time_err,&time_mean](const double d){
time_err += (d-time_mean)*(d-time_mean);});
time_err = std::sqrt(time_err/(runtimes.size()*(runtimes.size()-1)));
std::cout << num_threads << ' ' << time_mean << ' ' << time_err << std::endl;
}
}
}
#include <vector>
#include <algorithm>
#include <functional>
#include <cmath>
#include <random>
#include <cassert>
#include <iostream>
#include <iterator>
template <class BC>
class RandomWalk {
typedef std::mt19937 generator_type;
public:
RandomWalk ( const std::vector<double> ipos_, const double d_ )
: ipos (ipos_)
, d (d_)
{ }
double operator () (const std::size_t seed) {
std::vector<double> pos (ipos);
generator_type generator (seed);
std::function <double()> rng = std::bind ( std::uniform_real_distribution <> (0, 1), std::ref (generator) );
double direction;
while ( bc.inside (pos) ) {
direction = rng ();
pos [0] += d * std::cos (2 * M_PI * direction);
pos [1] += d * std::sin (2 * M_PI * direction);
}
return bc.evaluate ( pos );
}
private:
BC bc;
const std::vector<double> ipos;
double d;
};
template <class T>
class MC {
public:
MC (const T& func_, std::size_t seed_)
: func (func_)
, seed (seed_)
{}
void sample (const std::size_t N) {
sum = sum2 = 0.0;
for (std::size_t i = 0; i < N; ++i) {
double value = func (seed * N + i);
sum += value;
sum2 += value * value;
}
}
double mean (const std::size_t N) const
{
return sum / double (N);
}
double error (const std::size_t N) const
{
return std::sqrt ( ( sum2 - sum * sum / double(N) ) / double ( N * (N-1) ) );
}
private:
const std::function <double(std::size_t)> func;
const std::size_t seed;
double sum, sum2;
};
struct BC {
const double a = 1;
const double b = 1;
bool inside (const std::vector<double> &pos) const {
assert (pos.size() == 2);
return pos[0] > 0 && pos[0] < a && pos [1] > 0 && pos [1] < b;
}
virtual double evaluate (const std::vector<double>& pos) const { return 0; }
double project (const std::vector<double>& pos_) {
std::vector<double> pos (pos_);
if ( pos [0] < 0 ) pos [0] = 0;
if ( pos [1] < 0 ) pos [1] = 0;
if ( pos [0] > a ) pos [0] = a;
if ( pos [1] > b ) pos [1] = b;
return evaluate (pos);
}
};
struct BC_Dirichlet1 : BC {
double evaluate (const std::vector<double>& pos) const {
assert (pos.size() == 2);
return pos[0];
}
};
int main()
{
const double d = 0.01;
const std::vector<double> ipos = {0.3, 0.3};
const std::size_t reps = 16;
const std::size_t seed = 1;
RandomWalk < BC_Dirichlet1 > rw (ipos, d);
MC < RandomWalk <BC_Dirichlet1> > mc (rw, seed);
std::cout << "# samples mean error" << std::endl;
for (std::size_t i = 1; i < reps; i++) {
const std::size_t N = std::size_t(1) << i;
mc.sample(N);
std::cout << N << ' ' << mc.mean(N) <<