Skip to content

Commit c7bf1f5

Browse files
IthanilIthanil
committed
First kinda working (result correct, but errors) version of MPI wrapper for MCI + ex1 version with MPI
1 parent 37d047b commit c7bf1f5

File tree

4 files changed

+290
-0
lines changed

4 files changed

+290
-0
lines changed

‎config_template.sh

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,9 @@ LIBNAME="mci"
66
# C++ compiler
77
CC="g++"
88

9+
# MPI compiler wrapper
10+
MPICC="mpic++"
11+
912
# C++ flags (std=c++11 is necessary)
1013
FLAGS="-std=c++11 -Wall -Werror"
1114

‎examples/testmpi/main.cpp

Lines changed: 196 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,196 @@
1+
#include "mpi.h"
2+
#include <iostream>
3+
#include <cmath>
4+
#include <math.h>
5+
#include <fstream>
6+
7+
#include "MCIntegrator.hpp"
8+
#include "MCIWrapperMPI.hpp"
9+
10+
11+
// Observable functions
12+
class Parabola: public MCIObservableFunctionInterface{
13+
public:
14+
Parabola(const int &ndim): MCIObservableFunctionInterface(ndim, 1) {}
15+
16+
void observableFunction(const double * in, double * out){
17+
out[0] = 4.*in[0] - in[0]*in[0];
18+
}
19+
};
20+
21+
class NormalizedParabola: public MCIObservableFunctionInterface{
22+
public:
23+
NormalizedParabola(const int &ndim): MCIObservableFunctionInterface(ndim, 1) {}
24+
25+
void observableFunction(const double * in, double * out){
26+
out[0] = (4. - in[0]) * 5.;
27+
if (std::signbit(in[0])) out[0] = -out[0];
28+
}
29+
};
30+
31+
32+
33+
// Sampling function
34+
// the 48 is for normalization (even if not strictly necessary)
35+
class NormalizedLine: public MCISamplingFunctionInterface{
36+
public:
37+
NormalizedLine(const int &ndim): MCISamplingFunctionInterface(ndim, 1) {}
38+
39+
void samplingFunction(const double * in, double * protovalue){
40+
protovalue[0] = 0.2 * abs(in[0]);
41+
}
42+
43+
double getAcceptance(const double * protoold, const double * protonew){
44+
return protonew[0] / protoold[0];
45+
}
46+
};
47+
48+
49+
50+
int main() {
51+
using namespace std;
52+
53+
// intro
54+
cout << "We want to compute the integral" << endl;
55+
cout << " Integral[-1:3] dx (4x - x^2)" << endl << endl;
56+
cout << "Notice that we can re-write the integral as" << endl;
57+
cout << " Integral[-1:3] dx (5*sign(x)*(4-x) * |x|/5)" << endl;
58+
cout << "where g(x)=|x|/5 is a candidate sampling function since it's positive on the domain [-1:3] and its integral on the domain is equal to 1." << endl << endl;
59+
60+
cout << "We start by initializing MPI and setting the MCI:" << endl;
61+
62+
MPIMCI::init();
63+
64+
// --- From now on we run in parallel ---
65+
66+
int myrank = MPI::COMM_WORLD.Get_rank();
67+
68+
const int ndim = 1;
69+
MCI * mci = new MCI(ndim);
70+
71+
if (myrank == 0) cout << "ndim = " << mci->getNDim() << endl;
72+
73+
// set the integration range to [-1:3]
74+
double ** irange = new double*[ndim];
75+
irange[0] = new double[2];
76+
irange[0][0] = -1.;
77+
irange[0][1] = 3.;
78+
mci->setIRange(irange);
79+
80+
if (myrank == 0) cout << "irange = [ " << mci->getIRange(0, 0) << " ; " << mci->getIRange(0, 1) << " ]" << endl;
81+
82+
// initial walker position
83+
double * initpos = new double[ndim];
84+
initpos[0] = -0.5;
85+
mci->setX(initpos);
86+
87+
if (myrank == 0) cout << "initial walker position = " << mci->getX(0) << endl;
88+
89+
90+
// initial MRT2 step
91+
double * step = new double[ndim];
92+
step[0] = 0.5;
93+
mci->setMRT2Step(step);
94+
95+
if (myrank == 0) cout << "MRT2 step = " << mci->getMRT2Step(0) << endl;
96+
97+
// target acceptance rate
98+
double * targetacceptrate = new double[1];
99+
targetacceptrate[0] = 0.7;
100+
mci->setTargetAcceptanceRate(targetacceptrate);
101+
102+
if (myrank == 0) {
103+
cout << "Acceptance rate = " << mci->getTargetAcceptanceRate() << endl;
104+
cout << endl << endl;
105+
106+
// first way of integrating
107+
cout << "We first compute the integral without setting a sampling function." << endl;
108+
cout << " f(x) = (4x - x^2) " << endl;
109+
cout << " g(x) = - " << endl << endl;
110+
}
111+
112+
// observable
113+
MCIObservableFunctionInterface * obs = new Parabola(ndim);
114+
mci->addObservable(obs);
115+
116+
if (myrank == 0) {
117+
cout << "Number of observables set = " << mci->getNObs() << endl;
118+
// sampling function
119+
cout << "Number of sampling function set = " << mci->getNSampF() << endl;
120+
}
121+
122+
// integrate
123+
const long Nmc = 1000000;
124+
double * average;
125+
double * error;
126+
if (myrank == 0) { // allocate only for root
127+
average = new double[mci->getNObsDim()];
128+
error = new double[mci->getNObsDim()];
129+
}
130+
MPIMCI::integrate(mci, Nmc, average, error);
131+
132+
if (myrank == 0) {
133+
cout << "The integral gives as result = " << average[0] << " +- " << error[0] << endl;
134+
cout << "--------------------------------------------------------" << endl << endl;
135+
136+
// second way of integrating
137+
cout << "Now we compute the integral using a sampling function." << endl;
138+
cout << " f(x) = 5*sign(x)*(4-x) " << endl;
139+
cout << " g(x) = |x|/5 " << endl << endl;
140+
}
141+
142+
// observable
143+
delete obs;
144+
obs = new NormalizedParabola(ndim);
145+
mci->clearObservables(); // we first remove the old observable
146+
mci->addObservable(obs);
147+
148+
if (myrank == 0) cout << "Number of observables set = " << mci->getNObs() << endl;
149+
150+
151+
// sampling function
152+
MCISamplingFunctionInterface * sf = new NormalizedLine(ndim);
153+
mci->addSamplingFunction(sf);
154+
155+
if (myrank == 0) cout << "Number of sampling function set = " << mci->getNSampF() << endl;
156+
157+
158+
// integrate
159+
MPIMCI::integrate(mci, Nmc, average, error);
160+
161+
// deallocate per-thread allocations
162+
delete sf;
163+
164+
delete obs;
165+
166+
delete[] targetacceptrate;
167+
168+
delete[] step;
169+
170+
delete[] initpos;
171+
172+
delete[] irange[0];
173+
delete[] irange;
174+
175+
MPIMCI::finalize(mci); // also deletes all mci
176+
177+
178+
179+
// --- From now on we are single threaded again ---
180+
181+
cout << "The integral gives as result = " << average[0] << " +- " << error[0] << endl;
182+
cout << "--------------------------------------------------------" << endl << endl;
183+
184+
185+
// final comments
186+
cout << "Using a sampling function in this case gives worse performance. In fact, the error bar is larger." << endl;
187+
cout << "This implies that the variance of the re-factored f(x) written for introducing a sampling function, is larger than the original f(x)." << endl;
188+
189+
190+
// deallocate
191+
delete[] average;
192+
delete[] error;
193+
194+
// end
195+
return 0;
196+
}

‎examples/testmpi/run.sh

Lines changed: 40 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,40 @@
1+
#!/bin/bash
2+
source ../../config.sh
3+
OS_NAME=$(uname)
4+
5+
# Delete old compiled files
6+
\rm -f exe
7+
\rm -f *.o
8+
9+
#runtime dynamic library path
10+
RPATH="$(pwd)/../.."
11+
12+
FLAGS_TO_USE=$OPTFLAGS
13+
14+
# Build the main executable
15+
echo "$MPICC $FLAGS $FLAGS_TO_USE -I$(pwd)/../../src/ -c *.cpp"
16+
$MPICC $FLAGS $FLAGS_TO_USE -Wall -I$(pwd)/../../src/ -c *.cpp
17+
18+
case ${OS_NAME} in
19+
"Darwin")
20+
echo "$MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -o exe *.o -l${LIBNAME}"
21+
$MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -o exe *.o -l${LIBNAME}
22+
;;
23+
"Linux")
24+
echo "$MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -Wl,-rpath=${RPATH} -o exe *.o -l${LIBNAME}"
25+
$MPICC $FLAGS $FLAGS_TO_USE -L$(pwd)/../.. -Wl,-rpath=${RPATH} -o exe *.o -l${LIBNAME}
26+
;;
27+
esac
28+
29+
echo "Rebuilt the executable file"
30+
echo ""
31+
echo ""
32+
33+
# Run the debugging executable
34+
echo "Ready to run!"
35+
echo ""
36+
echo "--------------------------------------------------------------------------"
37+
echo ""
38+
echo ""
39+
echo ""
40+
mpirun -np 4 exe

‎src/MCIWrapperMPI.hpp

Lines changed: 51 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,51 @@
1+
#include <mpi.h>
2+
#include <iostream>
3+
#include "MCIntegrator.hpp"
4+
5+
namespace MPIMCI
6+
{
7+
void init()
8+
{
9+
if (MPI::Is_initialized()) throw std::runtime_error("MPI already initialized!");
10+
MPI::Init();
11+
MPI::COMM_WORLD.Set_errhandler(MPI::ERRORS_THROW_EXCEPTIONS);
12+
}
13+
14+
void integrate(MCI * const mci, const long &Nmc, double * average, double * error)
15+
{
16+
// make sure the user has MPI in the correct state
17+
if (!MPI::Is_initialized()) throw std::runtime_error("MPI not initialized!");
18+
if (MPI::Is_finalized()) throw std::runtime_error("MPI already finalized!");
19+
20+
const int myrank = MPI::COMM_WORLD.Get_rank();
21+
const int nranks = MPI::COMM_WORLD.Get_size();
22+
23+
// the results are stored in myAverage/Error and then reduced into average/error for root process
24+
double myAverage[mci->getNObsDim()];
25+
double myError[mci->getNObsDim()];
26+
27+
mci->integrate(Nmc, myAverage, myError);
28+
29+
for (int i=0; i<mci->getNObsDim(); ++i) {
30+
MPI::COMM_WORLD.Reduce(&myAverage[i], &average[i], 1, MPI::DOUBLE, MPI::SUM, 0);
31+
32+
myError[i] *= myError[i];
33+
MPI::COMM_WORLD.Reduce(&myError[i], &error[i], 1, MPI::DOUBLE, MPI::SUM, 0);
34+
35+
if (myrank == 0) {
36+
average[i] /= nranks;
37+
error[i] = sqrt(error[i]) / nranks;
38+
}
39+
}
40+
}
41+
42+
void finalize(MCI * mci)
43+
{
44+
// make sure the user has MPI in the correct state
45+
if (!MPI::Is_initialized()) throw std::runtime_error("MPI not initialized!");
46+
if (MPI::Is_finalized()) throw std::runtime_error("MPI already finalized!");
47+
48+
delete mci;
49+
MPI::Finalize();
50+
}
51+
};

0 commit comments

Comments
 (0)