Я моделирую стохастическое дифференциальное уравнение методом Монте-Карло, который в принципе идеально подходит для openMP, поскольку разные реализации не зависят друг от друга. К сожалению, я столкнулся с некоторыми проблемами с моим кодом, который дает неверный результат, как только я включаю openMP. Без этого все работает отлично. Мой «критический» цикл выглядит так:
double price = 0.0
#pragma omp parallel for private(VOld, VNew)
for (long i = 0; i < NSim; ++i){
VOld = S_0;
for (long index = 0; index < Nt; ++index){
VNew = VOld + (dt * r * VOld) + (sqrdt * sig * VOld * dW());
VOld = VNew;
}
double tmp = myOption.PayOff(VNew);
price += (tmp)/double(NSim);
}
Буду очень признателен за любую помощь. Заранее спасибо 🙂
Распространенная ошибка — забывать, что у каждого потока должен быть свой генератор случайных чисел. Если это не так, то каждый вызов dW будет портить внутреннее состояние (общего, а не частного) генератора случайных чисел.
Надеюсь, это поможет.
Ну, одна проблема, которую я вижу, состоит в том, что у вас есть условие гонки по переменной price
, Вы должны делать сокращение
#pragma omp parallel for private(VOld, VNew) reduction(+:price)
То же самое касается вашей переменной OptionPrice
Кроме того, мне кажется, что rng все еще общедоступный, а не приватный. Вы должны определить его в параллельном блоке, если вы хотите, чтобы он был закрытым или объявить его закрытым (для закрытых переменных я предпочитаю объявлять их в параллельном блоке, который автоматически делает их закрытыми, а не объявлять их закрытыми).
Итак, основываясь на ответах @jmbr и @raxman, я переместил внутренний цикл в отдельную функцию и убедился, что rng
сейчас действительно личное. Также обратите внимание на трюк с посевом, который становится жизненно важным. Кроме того, я представил reduction
на OptionPrice
, Код ниже работает нормально.
double SimulateStockPrice(const double InitialPrize, const double dt, const long Nt, const double r, const double sig, boost::mt19937 *rng){
static unsigned long seed = 0;
boost::mt19937 *rng = new boost::mt19937();
rng -> seed((++seed) + time(NULL));
boost::normal_distribution<> nd(0.0, 1.0);
boost::variate_generator< boost::mt19937, boost::normal_distribution<> > dW(*rng, nd);
double sqrdt = sqrt(dt);
double PriceNew(0.0), PriceOld(InitialPrize);
for (long index = 0; index < Nt; ++index){
PriceNew = PriceOld + (dt * r * PriceOld) + (sqrdt * sig * PriceOld * dW());
PriceOld = PriceNew;
}
delete rng;
return PriceNew;
}
Затем в большой петле я иду с:
#pragma omp parallel for default(none) shared(dt, NSim, Nt, S_0, myOption) reduction(+:OptionPrice)
for (long i = 0; i < NSim; ++i){
double StockPrice = SimulateStockPrice(S_0, dt, Nt, myOption.r, myOption.sig, rng);
double PayOff = myOption.myPayOffFunction(StockPrice);
OptionPrice += PayOff;
}
И пошло-поехало 🙂