我目前正在寫一個物理模擬（t.b.m.p.求解一個隨機微分方程），我需要並行化它。
現在這可以通過MPI來實現，我想我將來需要做一些工作，但是現在我想利用我的本地機器的所有8個內核。對於一個參數集，正常運行需要2-17小時。所以我想利用多線程，特別是下面的函數應該並行執行。對於Nsteps時間步長，此功能基本上解決了相同的SDE Nrep次。將結果取平均值並將每個線程存儲到Nthreads x Nsteps陣列JpmArr的單獨行中。多線程和動態數組/矩陣運算的
double **JpmArr;
void worker(const dtype gamma, const dtype dt, const uint seed, const uint Nsteps, const uint Nrep,\
const ESpMatD& Jplus, const ESpMatD& Jminus, const ESpMatD& Jz, const uint tId){
dtype dW(0), stdDev(sqrt(dt));
std::normal_distribution<> WienerDistr(0, stdDev);
//create the arrays for the values of <t|J+J-|t>
dtype* JpmExpect = JpmArr[tId];
//execute Nrep repetitions of the experiment
for (uint r(0); r < Nrep; ++r) {
//reinitialize the wave function
psiVecArr[tId] = globalIstate;
//<t|J+J-|t>
tmpVecArr[tId] = Jminus* psiVecArr[tId];
JpmExpect[0] += tmpVecArr[tId].squaredNorm();
//iterate over the timesteps
for (uint s(1); s < Nsteps; ++s) {
//get a random number
dW = WienerDistr(RNGarr[tId]);
//execute one step of the RK-s.o. 1 scheme
tmpPsiVecArr[tId] = F2(gamma, std::ref(Jminus), std::ref(psiVecArr[tId]));
tmpVecArr[tId] = psiVecArr[tId] + tmpPsiVecArr[tId] * sqrt(dt);
psiVecArr[tId] = psiVecArr[tId] + F1(gamma, std::ref(Jminus), std::ref(Jplus), std::ref(psiVecArr[tId])) * dt + tmpPsiVecArr[tId] * dW \
+ 0.5 * (F2(gamma, std::ref(Jminus), std::ref(tmpVecArr[tId])) - F2(gamma, std::ref(Jminus), std::ref(psiVecArr[tId]))) *(dW * dW - dt)/sqrt(dt);
//normalise
psiVecArr[tId].normalize();
//compute <t|J+J-|t>
tmpVecArr[tId] = Jminus* psiVecArr[tId];
JpmExpect[s] += tmpVecArr[tId].squaredNorm();
}
}
//average over the repetitions
for (uint j(0); j < Nsteps; ++j) {
JpmExpect[j] /= Nrep;
}
}
我使用Eigen作爲線性代數庫，因而：
typedef Eigen::SparseMatrix<dtype, Eigen::RowMajor> ESpMatD;
typedef Eigen::Matrix<dtype, Eigen::Dynamic, Eigen::RowMajor> VectorXdrm;
被用作類型。上述工函數調用：
VectorXdrm& F1(const dtype a, const ESpMatD& A, const ESpMatD& B, const VectorXdrm& v) {
z.setZero(v.size());
y.setZero(v.size());
// z is for simplification
z = A*v;
//scalar intermediate value c = <v, Av>
dtype c = v.dot(z);
y = a * (2.0 * c * z - B * z - c * c * v);
return y;
}
VectorXdrm& F2(const dtype a, const ESpMatD& A, const VectorXdrm& v) {
//zero the data
z.setZero(v.size());
y.setZero(v.size());
z = A*v;
dtype c = v.dot(z);
y = sqrt(2.0 * a)*(z - c * v);
return y;
}
其中矢量z,yVectorXdrm是類型，並且在相同的文件（模塊全局）被聲明。 所有的數組（RNGarr, JpmArr, tmpPsiVecArr, tmpVecArr, psiVecArr）都在main中初始化（使用main.cpp中的extern聲明）。在完成設置後，我使用std::async運行函數，等待所有完成，然後將JpmArr中的數據收集到main()中的單個數組中，並將其寫入文件。
問題： 如果我使用std::launch::async，結果是無稽之談。 如果我使用std::launch::deferred，計算結果和平均結果匹配（直到數值方法允許）我通過分析手段獲得的結果。
我不知道東西失敗的地方。我曾經使用Armadillo作爲線性代數，但它的normalize例程交付了nan's，所以我切換到Eigen，它提示（在文檔中）可用於多個線程 - 它仍然失敗。 在我花了4天的時間之前沒有和線程一起工作，現在試圖讓這個工作和閱讀的東西。後者導致我使用全局數組RNGarr, JpmArr, tmpPsiVecArr, tmpVecArr, psiVecArr（之前我只是試圖在worker中創建合適的數組，並將結果通過struct workerResult返回到main。）以及使用std::ref()將矩陣Jplus, Jminus, Jz傳遞給該函數。 （爲了簡潔起見，上面的功能被省略了） -
但是我得到的結果仍然是錯誤的，我不再有任何想法了，什麼是錯的，我應該怎麼做才能獲得正確的結果。任何輸入和/或指向此類（線程）問題或參考解決方案示例的指針都將不勝感激。
你能展示實際調用啓動:: async的代碼？能夠運行多個線程的關鍵要求是所有的代碼都應該是線程安全的（這很難確定，但是在C++中改變爲true比在C中更好），並且每個線程應該在其自己的數據結構上工作（或一些程序上的保證，他們沒有彼此踩腳趾） – Soren
'auto ws = std :: async（std :: launch :: deferred，＆worker，gamma，dt，RNGseed [s]，Nsteps，NrepPt，std :: ref（Jplus），std :: ref（Jminus），std :: ref（Jz），tId ++）;'哪裏s'是0,1,2,3現在手動設置 – Nox
Jplus，Jminus - 它們是引用其他矩陣？在工作人員正在運行時他們是否會被另一個線程更新？ – Soren