From 3e31ed1254362dd89c7760638296a07c623d0e90 Mon Sep 17 00:00:00 2001 From: Anatoly Date: Mon, 24 Jun 2019 19:36:24 +0300 Subject: [PATCH] trying to find the parralel answer --- src/CMakeLists.txt | 1 + src/spacetimecorr.cpp | 157 ++++++++++++++++++++++++++++++++---------- 2 files changed, 122 insertions(+), 36 deletions(-) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 295035e..eaa54df 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,4 +1,5 @@ set(GROMACS_CXX_FLAGS "${GROMACS_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") +set( CMAKE_CXX_FLAGS "-fsanitize=thread -O2 -g") add_executable(spacetimecorr spacetimecorr.cpp) include_directories( ${GROMACS_INCLUDE_DIRS} diff --git a/src/spacetimecorr.cpp b/src/spacetimecorr.cpp index d3b13ef..2466e9d 100644 --- a/src/spacetimecorr.cpp +++ b/src/spacetimecorr.cpp @@ -134,13 +134,56 @@ void make_rout_file(double crl_border, std::vector< int > indx, std::vector< std std::fprintf(file, "correlations >= %0.2f\n\n", crl_border); for (unsigned int i = 0; i < rout.size(); i++) { for (unsigned int j = 0; j < rout[i].size(); j++) { - std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.15\n", indx[rout[i][j].first]/* + 1*/, indx[rout[i][j].second]/* + 1*/); + std::fprintf(file, "cgo_arrow (id %3d), (id %3d), radius=0.15\n", indx[rout[i][j].first] + 1, indx[rout[i][j].second] + 1); } std::fprintf(file, "\n\n"); } std::fclose(file); } +void make_table_routs(double crl_border, std::vector< int > indx, std::vector< std::vector< std::pair< unsigned int, unsigned int > > > rout, const char* file_name) +{ + FILE *file; + file = std::fopen(file_name, "w+"); + std::fprintf(file, "correlations >= %0.2f\n\n", crl_border); + std::vector< std::tuple< int, int, std::vector< int > > > table; + + for (unsigned int i1 = 0; i1 < rout.size(); i1++) { + for (unsigned int i2 = 0; i2 < rout[i1].size(); i2++) { + bool flag = true; + for (unsigned int i3 = 0; i3 < table.size(); i3++) { + if (std::get<0>(table[i3]) == indx[rout[i1][i2].first] + 1) { + std::get<1>(table[i3])++; + std::get<2>(table[i3]).push_back(indx[rout[i1][i2].second] + 1); + flag = false; + } else if (std::get<0>(table[i3]) == indx[rout[i1][i2].second] + 1) { + std::get<1>(table[i3])++; + std::get<2>(table[i3]).push_back(indx[rout[i1][i2].first] + 1); + flag = false; + } + } + if (flag) { + std::vector< int > a; + a.resize(0); + a.push_back(indx[rout[i1][i2].second] + 1); + table.push_back(std::make_tuple(indx[rout[i1][i2].first] + 1, 1, a)); + a.resize(0); + a.push_back(indx[rout[i1][i2].first] + 1); + table.push_back(std::make_tuple(indx[rout[i1][i2].second] + 1, 1, a)); + } + } + } + + for (unsigned int i = 0; i < table.size(); i++) { + std::fprintf(file, "atomN %d connections %d | ", std::get<0>(table[i]), std::get<1>(table[i])); + for (unsigned int j = 0; j < std::get<2>(table[i]).size(); j++) { + std::fprintf(file, "%d ", std::get<2>(table[i])[j]); + } + std::fprintf(file, "\n"); + } + std::fclose(file); +} + void make_best_corrs_graphics(std::vector< std::vector< std::vector< double > > > correlations, std::vector< std::vector< std::pair< unsigned int, unsigned int > > > rout_pairs, std::vector< int > indx, @@ -198,49 +241,89 @@ const std::vector< std::vector< std::vector< double > > > correlation_evaluation crl.resize(tauE + 1); for (unsigned int i = 0; i < crl.size(); i++) { crl[i].resize(traj.front().size()); - for (unsigned int j = 0; j < crl[i].size(); j++) { - crl[i][j].resize(traj.front().size(), 0); + for (unsigned int j = 0; j < traj.front().size(); j++) { + //crl[i][j].resize(traj.front().size(), 0); + for (unsigned int k = 0; k < traj.front().size(); k++) { + crl[i][j].push_back(0.); + } } } - #pragma omp parallel for schedule(dynamic) shared(crl, traj, ref) - for (unsigned int i = 0; i <= tauE; i++) { - RVec temp1, temp2; - std::vector< std::vector< double > > a, b, c; - std::vector< double > d; - d.resize(traj.front().size(), 0); - a.resize(traj.front().size(), d); - b.resize(traj.front().size(), d); - c.resize(traj.front().size(), d); - for (unsigned int j = 0; j < traj.size() - i - 1; j++) { - for (unsigned int f1 = 0; f1 < traj.front().size(); f1++) { - for (unsigned int f2 = 0; f2 < traj.front().size(); f2++) { - temp1 = traj[j][f1] - ref[f1]; - temp2 = traj[j + i][f2] - ref[f2]; - a[f1][f2] += (static_cast(temp1[0]) * static_cast(temp2[0]) + - static_cast(temp1[1]) * static_cast(temp2[1]) + - static_cast(temp1[2]) * static_cast(temp2[2])); - b[f1][f2] += (static_cast(temp1[0]) * static_cast(temp1[0]) + - static_cast(temp1[1]) * static_cast(temp1[1]) + - static_cast(temp1[2]) * static_cast(temp1[2])); - c[f1][f2] += (static_cast(temp2[0]) * static_cast(temp2[0]) + - static_cast(temp2[1]) * static_cast(temp2[1]) + - static_cast(temp2[2]) * static_cast(temp2[2])); + int trjsize = traj.size(), trjsize2 = traj.front().size(); + + #pragma omp parallel for ordered schedule(dynamic) shared(crl, traj, ref, trjsize, trjsize2) + for (unsigned int i = 0; i <= tauE; i++) { + + #pragma omp ordered + { + + + RVec temp1, temp2; + std::vector< std::vector< double > > a, b, c; + std::vector< double > d; + d.resize(0); + for (unsigned int j = 0; j < traj.front().size(); j++) { + a.push_back(d); + b.push_back(d); + c.push_back(d); + for (unsigned int k = 0; k < traj.front().size(); k++) { + a[j].push_back(0.); + b[j].push_back(0.); + c[j].push_back(0.); + } + } + /*d.resize(traj.front().size(), 0); + a.resize(traj.front().size(), d); + b.resize(traj.front().size(), d); + c.resize(traj.front().size(), d);*/ + + for (unsigned int j = 0; j < trjsize - i - 1; j++) { + for (unsigned int f1 = 0; f1 < trjsize2; f1++) { + for (unsigned int f2 = 0; f2 < trjsize2; f2++) { +#pragma omp critical +{ + RVec temp3 = traj[j][f1]; + RVec temp4 = ref[f1]; + temp1 = temp3 - temp4; + + temp3 = traj[j + i][f2]; + temp4 = ref[f2]; + + temp2 = temp3 - temp4; + + //temp1 = traj[j][f1] - ref[f1]; + //temp2 = traj[j + i][f2] - ref[f2]; + + a[f1][f2] += (static_cast(temp1[0]) * static_cast(temp2[0]) + + static_cast(temp1[1]) * static_cast(temp2[1]) + + static_cast(temp1[2]) * static_cast(temp2[2])); + b[f1][f2] += (static_cast(temp1[0]) * static_cast(temp1[0]) + + static_cast(temp1[1]) * static_cast(temp1[1]) + + static_cast(temp1[2]) * static_cast(temp1[2])); + c[f1][f2] += (static_cast(temp2[0]) * static_cast(temp2[0]) + + static_cast(temp2[1]) * static_cast(temp2[1]) + + static_cast(temp2[2]) * static_cast(temp2[2])); +} + } } } - } - #pragma omp critical - for (unsigned int j = 0; j < traj.front().size(); j++) { - for (unsigned int f = 0; f < traj.front().size(); f++) { - crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f])); + #pragma omp critical + { + for (unsigned int j = 0; j < traj.front().size(); j++) { + for (unsigned int f = 0; f < traj.front().size(); f++) { + crl[i][j][f] = a[j][f] / (std::sqrt(b[j][f] * c[j][f])); + } + } } + std::cout << i << " | "; + if (i % 100 == 0) { + std::cout << "\n"; + } + + } + } - std::cout << i << " | "; - if (i % 100 == 0) { - std::cout << "\n"; - } - } #pragma omp barrier std::cout << "\n"; return crl; @@ -647,6 +730,8 @@ SpaceTimeCorr::finishAnalysis(int nframes) make_best_corrs_graphics(crltns, rout_new, index, (OutPutName + "_routs_graphics.txt").c_str()); std::cout << "corelation routs' pairs' graphics printed\n"; + make_table_routs(static_cast< double >(crl_border), index, rout_new, (OutPutName + "_routs_table.txt").c_str()); + /*std::cout << "extra params\n"; std::vector< double > diffusion; evaluate_diffusion(trajectory, diffusion); -- 2.22.0