VM2D  1.12
Vortex methods for 2D flows simulation
knnCPU.cpp
Go to the documentation of this file.
1 
2 #include "knnCPU.h"
3 
4 #include <omp.h>
5 #include <algorithm>
6 
7 
8 const int codeLength = 15;
9 const int twoPowCodeLength = (1 << codeLength);
10 
11 
13 std::pair<Point2D, Point2D> GetGabarit(const std::vector<Vortex2D>& vtx)
14 {
15  auto xx = std::minmax_element(vtx.begin(), vtx.end(), [](const Vortex2D& P1, const Vortex2D& P2) {return P1.r()[0] < P2.r()[0]; });
16  auto yy = std::minmax_element(vtx.begin(), vtx.end(), [](const Vortex2D& P1, const Vortex2D& P2) {return P1.r()[1] < P2.r()[1]; });
17 
18  return { { xx.first->r()[0], yy.first->r()[1] }, { xx.second->r()[0], yy.second->r()[1] } };
19 }
20 
21 
22 
23 //"Разрежение" двоичного представления беззнакового целого, вставляя по одному нулику между всеми битами
24 unsigned int ExpandBits(unsigned int v)
25 {
26  // вставит 1 нуль
27  v = (v | (v << 8)) & 0x00FF00FF; // 00000000`00000000`abcdefgh`ijklmnop
28  // | 00000000`abcdefgh`ijklmnop`00000000
29  // = 00000000`abcdefgh`XXXXXXXX`ijklmnop
30  // & 00000000`11111111`00000000`11111111
31  // = 00000000`abcdefgh`00000000`ijklmnop
32 
33  v = (v | (v << 4)) & 0x0F0F0F0F; // 00000000`abcdefgh`00000000`ijklmnop
34  // | 0000abcd`efgh0000`0000ijkl`mnop0000
35  // = 0000abcd`XXXXefgh`0000ijkl`XXXXmnop
36  // & 00001111`00001111`00001111`00001111
37  // = 0000abcd`0000efgh`0000ijkl`0000mnop
38 
39  v = (v | (v << 2)) & 0x33333333; // 0000abcd`0000efgh`0000ijkl`0000mnop
40  // | 00abcd00`00efgh00`00ijkl00`00mnop00
41  // = 00abXXcd`00efXXgh`00ijXXkl`00mnXXop
42  // & 00110011`00110011`00110011`00110011
43  // = 00ab00cd`00ef00gh`00ij00kl`00mn00op
44 
45  v = (v | (v << 1)) & 0x55555555; // 00ab00cd`00ef00gh`00ij00kl`00mn00op
46  // | 0ab00cd0`0ef00gh0`0ij00kl0`0mn00op0
47  // = 0aXb0cXd`0eXf0gXh`0iXj0kXl`0mXn0oXp
48  // & 01010101`01010101`01010101`01010101
49  // = 0a0b0c0d`0e0f0g0h`0i0j0k0l`0m0n0o0p
50  return v;
51 }
52 
53 
54 //Мортоновский код для пары из чисел типа double
55 //Исходное число - строго в диапазоне [0, 1 / (1.0 + 1/(2^codeLength - 1) ))
56 unsigned int Morton2D(const Point2D& r)
57 {
58  const Point2D& rscale = twoPowCodeLength * r;
59  const unsigned int& xx = ExpandBits((unsigned int)(rscale[0]));
60  const unsigned int& yy = ExpandBits((unsigned int)(rscale[1]));
61  return yy | (xx << 1);
62 }
63 
64 
66 void RSort_Parallel(TParticleCode* m, TParticleCode* m_temp, unsigned int n, unsigned int* s)
67 {
68  if (n == 0)
69  return;
70  // Количество задействованных потоков unsigned char threads = omp_get_max_threads(); //std::cout << "threads = " << (int)threads << std::endl; #pragma omp parallel num_threads(threads) { TParticleCode* source = m; TParticleCode* dest = m_temp; unsigned int l = omp_get_thread_num(); unsigned int div = n / omp_get_num_threads(); unsigned int mod = n % omp_get_num_threads(); unsigned int left_index = l < mod ? (div + (mod == 0 ? 0 : 1)) * l : n - (omp_get_num_threads() - l) * div; unsigned int right_index = left_index + div - (mod > l ? 0 : 1); for (unsigned int digit = 0; digit < sizeof(m->key); ++digit) { unsigned int s_sum[256] = { 0 }; unsigned int s0[256] = { 0 }; unsigned char* b1 = (unsigned char*)&source[right_index].key; unsigned char* b2 = (unsigned char*)&source[left_index].key; while (b1 >= b2) { ++s0[*(b1 + digit)]; b1 -= sizeof(TParticleCode); } for (unsigned int i = 0; i < 256; i++) { s[i + 256 * l] = s0[i]; } #pragma omp barrier for (unsigned int j = 0; j < threads; j++) { for (unsigned int i = 0; i < 256; i++) { s_sum[i] += s[i + 256 * j]; if (j < l) { s0[i] += s[i + 256 * j]; } } } for (unsigned int i = 1; i < 256; ++i) { s_sum[i] += s_sum[i - 1]; s0[i] += s_sum[i - 1]; } unsigned char* b = (unsigned char*)&source[right_index].key + digit; TParticleCode* v1 = &source[right_index]; TParticleCode* v2 = &source[left_index]; while (v1 >= v2) { dest[--s0[*b]] = *v1--; b -= sizeof(TParticleCode); } #pragma omp barrier std::swap(source, dest); } } // Если ключ структуры однобайтовый, просто копируем в исходный массив if (sizeof(m->key) == 1) { memcpy(m, m_temp, n * sizeof(TParticleCode)); } } inline size_t BinSearch(const std::vector<std::pair<double, size_t>>& currentNN, double x, int low, int high) { int mid = -1; if (x > currentNN[high].first) return high; while (low <= high) { mid = (low + high) / 2; if (currentNN[mid].first == x) return mid + 1; //выход из цикла if (currentNN[mid].first < x) low = mid + 1; else high = mid - 1; } return mid; } void newSort(std::vector<std::pair<double, size_t>>& mass, std::vector<std::pair<double, size_t>>& dstKeys) { size_t k = mass.size(); size_t cnt = 0; for (size_t i = 0; i < k; i++) { double elem = mass[i].first; cnt = 0; for (size_t j = 0; j < k; j++) { cnt += (mass[j].first < elem); //std::cout << cnt; } //if (!AssumeDistinct) for (size_t j = 0; j < i; j++) cnt += (mass[j].first == elem); dstKeys[cnt] = mass[i]; } mass.swap(dstKeys); } void newMerge(std::vector<std::pair<double, size_t>>& currentNN, const std::vector<std::pair<double, size_t>>& candidateNN, std::vector<size_t>& loc, std::vector<size_t>& counter, std::vector<size_t>& offset, std::vector<size_t>& counterScan, std::vector<std::pair<double, size_t>>& updateNN ) { const size_t k = candidateNN.size() / 2; //количество ближайших соседей //std::cout << "loc: " << std::endl; for (size_t j = 0; j < 2 * k; ++j) { loc[j] = BinSearch(currentNN, candidateNN[j].first, 0, (int)k - 1); //std::cout << loc[j] << " "; } //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) counter[j] = j & 1; //std::cout << "counter: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counter[j] << " "; //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) { if ((loc[j] > 0) && ((loc[j] == k) || (candidateNN[j].first == currentNN[loc[j] - 1].first))) //можно через second, проверить!! offset[j] = k; else offset[j] = counter[2 * loc[j]]++; } //std::cout << "counter: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counter[j] << " "; //std::cout << std::endl; //std::cout << "offset: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << offset[j] << " "; //std::cout << std::endl; counterScan[0] = 0; for (size_t j = 1; j < 2 * k; ++j) counterScan[j] = counterScan[j - 1] + counter[j - 1]; //std::cout << "counterScan: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counterScan[j] << " "; //std::cout << std::endl; size_t index; for (size_t j = 0; j < k; ++j) { index = counterScan[2 * j + 1]; //std::cout << index << std::endl; if (index < k) updateNN[index] = currentNN[j]; } //std::cout << "updateNN: " << std::endl; //for (size_t j = 0; j < k; ++j) // std::cout << "(" << updateNN[j].first << ", " << updateNN[j].second << "), "; //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) { if (2 * loc[j] < 2 * k) { index = counterScan[2 * loc[j]] + offset[j]; if (index < k) updateNN[index] = candidateNN[j]; } } currentNN.swap(updateNN); } void mainCycle( const int k, const std::vector<Vortex2D>& vtx, const std::pair<Point2D, Point2D>& gab, const double scale, const int sdvig, std::vector<TParticleCode>& mortonCodeBoth, std::vector<unsigned int>& s, std::vector<TParticleCode>& mortonCodeBoth_temp, std::vector<TParticleCode>& mcdata, std::vector<TParticleCode>& mcquery, std::vector<unsigned int>& iq, std::vector<size_t>& initneig, std::vector<std::vector<std::pair<double, size_t>>>& initdist, double* time ) { //масштабирование координат вихрей в [0;0.75)^2, поиск их мортоновских кодов time[0] = omp_get_wtime(); #pragma omp parallel for for (int i = 0; i < vtx.size(); ++i) { //создание массива data mortonCodeBoth[i].key = Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 }); mortonCodeBoth[i].key <<= 1; //mortonCodeBoth[i].isQuery = 0; mortonCodeBoth[i].originNumber = i; //создание массива query mortonCodeBoth[vtx.size() + i].key = mortonCodeBoth[i].key;// Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 }); mortonCodeBoth[vtx.size() + i].key |= 1; mortonCodeBoth[vtx.size() + i].originNumber = i; //std::cout << "code_data[" << i << "] = " << std::bitset<8>(mortonCodeData[i].key) << std::endl; } time[1] = omp_get_wtime(); //std::cout << "----------" << std::endl; //for (size_t i = 0; i < mortonCodeBoth.size(); ++i) //{ // std::cout << "code[" << i << "] = " << std::bitset<2 * codeLength + 1>(mortonCodeBoth[i].key) << std::endl; //} //сортировка единого массива (data и query) по мортоновским кодам RSort_Parallel(mortonCodeBoth.data(), mortonCodeBoth_temp.data(), (int)mortonCodeBoth.size(), s.data()); time[2] = omp_get_wtime(); //для отладки (рандомно меняем местами data и query с одинаковыми мортоновскими кодами) //for (size_t i = 0; i < mortonCodeBoth.size() / 2; ++i) // if (rand() > RAND_MAX / 2) // std::swap(mortonCodeBoth[2 * i], mortonCodeBoth[2 * i + 1]); //std::cout << "----------" << std::endl; //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные) mcdata.resize(0); mcquery.resize(0); mcdata.reserve(vtx.size()); mcquery.reserve(vtx.size()); time[3] = omp_get_wtime(); iq.resize(0); iq.reserve(vtx.size()); for (size_t i = 0, j = -1; i < mortonCodeBoth.size(); ++i) { if (mortonCodeBoth[i].key & 1) { mcquery.push_back(mortonCodeBoth[i]); iq.push_back((unsigned)(j + 1) == 0 ? 0 : (unsigned int)j); } else { ++j; mcdata.push_back(mortonCodeBoth[i]); } } time[4] = omp_get_wtime(); //std::cout << "----------" << std::endl; // //for (size_t i = 0; i < iq.size(); ++i) // std::cout << "iq[" << i << "] = " << iq[i] << std::endl; // //std::cout << "----------" << std::endl; //for (size_t i = 0; i < mcdata.size(); ++i) // std::cout << "mcdata[" << i << "] = " << std::bitset<8>(mcdata[i].key) << ", isQuery = " << (int)mcdata[i].isQuery << std::endl; //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой) //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей #pragma omp parallel for for (int q = 0; q < mcquery.size(); ++q) { size_t pt = iq[q]; size_t left = pt - k + 1;// , right = pt + k; if (pt < k - 1) { left = 0; //right = 2 * k - 1; } if (pt > mcquery.size() - k - 1) { //right = mcquery.size() - 1; left = (mcquery.size() - 1) - (2 * k - 1); } initneig[q] = left; } time[5] = omp_get_wtime(); //std::cout << "----------" << std::endl; //for (size_t i = 0; i < initneig.size(); ++i) // std::cout << "initneig[" << i << "] = from " << initneig[i] << std::endl; //std::cout << "----------" << std::endl; std::vector<std::pair<double, size_t>> dist(2 * k, { -1.0, -1 }); //std::vector<std::pair<double, size_t>> mergeArray(3 * k); std::vector<size_t> loc(2 * k); std::vector<size_t> counter(2 * k, 0); std::vector<size_t> offset(2 * k, 0); std::vector<size_t> counterScan(2 * k, 0); std::vector<std::pair<double, size_t>> updateNN(k, { 0.0, 0 }); std::vector<std::pair<double, size_t>> dstKeys1(2 * k, { 0.0, 0 }); std::vector<std::pair<double, size_t>> dstKeys(k, { 0.0, 0 }); time[6] = omp_get_wtime(); if (sdvig == 0) #pragma omp parallel for firstprivate(dstKeys1) for (int i = 0; i < initdist.size(); ++i) { for (size_t j = 0; j < 2 * k; ++j) initdist[mcquery[i].originNumber][j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber }; //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; }); newSort(initdist[mcquery[i].originNumber], dstKeys1); initdist[mcquery[i].originNumber].resize(k); } else #pragma omp parallel for firstprivate(dist) firstprivate(loc, counter, offset, counterScan, updateNN, dstKeys)//, firstprivate(mergeArray) for (int i = 0; i < initdist.size(); ++i) { for (size_t j = 0; j < 2 * k; ++j) dist[j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber }; //std::sort(dist.begin(), dist.end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; }); //std::cout << "dist[" << i << "] = "; //for (size_t j = 0; j < k; ++j) // std::cout << "(" << dist[j].first << ", " << dist[j].second << "), "; //std::cout << std::endl; //Mymerge(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), dist.begin(), dist.end(), mergeArray.begin()); newMerge(initdist[mcquery[i].originNumber], dist, loc, counter, offset, counterScan, updateNN); newSort(initdist[mcquery[i].originNumber], dstKeys); //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end()); //std::copy(mergeArray.begin(), mergeArray.begin() + k, initdist[mcquery[i].originNumber].begin()); } time[7] = omp_get_wtime(); } void WakekNN(const std::vector<Vortex2D>& vtx, const size_t k, std::vector<std::vector<std::pair<double, size_t>>>& initdist) { const size_t n = vtx.size(); auto gab = GetGabarit(vtx); double scale = std::max(gab.second[0] - gab.first[0], gab.second[1] - gab.first[1]); //создание единого массива data и query std::vector<TParticleCode> mortonCodeBoth(2 * vtx.size()); //сортировка единого массива (data и query) по мортоновским кодам std::vector<unsigned int> s(256 * omp_get_max_threads()); std::vector<TParticleCode> mortonCodeBoth_temp(2 * vtx.size()); //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные) std::vector<TParticleCode> mcdata, mcquery; std::vector<unsigned int> iq; //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой) //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей std::vector<size_t> initneig(vtx.size(), -1); //выбор из 2*k кандидатов k соседей (выбор по расстоянию от данных до запроса) // initdist = {расстояние; до элемента данных} // сортировка по расстоянию /*std::vector < std::vector < std::pair<double, size_t> >> initdist(vtx.size()); for (auto& d : initdist) d.resize(2 * k, { -1.0, -1 });*/ //Сдвигаем все точки double t1 = omp_get_wtime(); double tStart[5], tFinish[5]; double time[5][8]; for (size_t sdvig = 0; sdvig <= 4; ++sdvig) //4 { tStart[sdvig] = omp_get_wtime(); mainCycle((int)k, vtx, gab, scale, (int)sdvig, mortonCodeBoth, s, mortonCodeBoth_temp, mcdata, mcquery, iq, initneig, initdist, time[sdvig]); tFinish[sdvig] = omp_get_wtime(); //std::cout << "time_knn:" // << " " << time[sdvig][1] - time[sdvig][0] // << " " << time[sdvig][2] - time[sdvig][1] // << " " << time[sdvig][3] - time[sdvig][2] // << " " << time[sdvig][4] - time[sdvig][3] // << " " << time[sdvig][5] - time[sdvig][4] // << " " << time[sdvig][6] - time[sdvig][5] // << " " << time[sdvig][7] - time[sdvig][6] << std::endl; }//for sdvig }
71  unsigned char threads = omp_get_max_threads();
72  //std::cout << "threads = " << (int)threads << std::endl;
73 
74 #pragma omp parallel num_threads(threads)
75  {
76  TParticleCode* source = m;
77  TParticleCode* dest = m_temp;
78  unsigned int l = omp_get_thread_num();
79  unsigned int div = n / omp_get_num_threads();
80  unsigned int mod = n % omp_get_num_threads();
81  unsigned int left_index = l < mod ? (div + (mod == 0 ? 0 : 1)) * l : n - (omp_get_num_threads() - l) * div;
82  unsigned int right_index = left_index + div - (mod > l ? 0 : 1);
83 
84  for (unsigned int digit = 0; digit < sizeof(m->key); ++digit)
85  {
86  unsigned int s_sum[256] = { 0 };
87  unsigned int s0[256] = { 0 };
88  unsigned char* b1 = (unsigned char*)&source[right_index].key;
89  unsigned char* b2 = (unsigned char*)&source[left_index].key;
90  while (b1 >= b2)
91  {
92  ++s0[*(b1 + digit)];
93  b1 -= sizeof(TParticleCode);
94  }
95  for (unsigned int i = 0; i < 256; i++)
96  {
97  s[i + 256 * l] = s0[i];
98  }
99 
100 #pragma omp barrier
101  for (unsigned int j = 0; j < threads; j++)
102  {
103  for (unsigned int i = 0; i < 256; i++)
104  {
105  s_sum[i] += s[i + 256 * j];
106  if (j < l)
107  {
108  s0[i] += s[i + 256 * j];
109  }
110  }
111  }
112 
113  for (unsigned int i = 1; i < 256; ++i)
114  {
115  s_sum[i] += s_sum[i - 1];
116  s0[i] += s_sum[i - 1];
117  }
118  unsigned char* b = (unsigned char*)&source[right_index].key + digit;
119  TParticleCode* v1 = &source[right_index];
120  TParticleCode* v2 = &source[left_index];
121  while (v1 >= v2)
122  {
123  dest[--s0[*b]] = *v1--;
124  b -= sizeof(TParticleCode);
125  }
126 #pragma omp barrier
127  std::swap(source, dest);
128  }
129  }
130 
131  // Если ключ структуры однобайтовый, просто копируем в исходный массив
132  if (sizeof(m->key) == 1)
133  {
134  memcpy(m, m_temp, n * sizeof(TParticleCode));
135  }
136 }
137 
138 
139 inline size_t BinSearch(const std::vector<std::pair<double, size_t>>& currentNN, double x, int low, int high)
140 {
141  int mid = -1;
142 
143  if (x > currentNN[high].first)
144  return high;
145 
146  while (low <= high) {
147  mid = (low + high) / 2;
148 
149  if (currentNN[mid].first == x)
150  return mid + 1; //выход из цикла if (currentNN[mid].first < x) low = mid + 1; else high = mid - 1; } return mid; } void newSort(std::vector<std::pair<double, size_t>>& mass, std::vector<std::pair<double, size_t>>& dstKeys) { size_t k = mass.size(); size_t cnt = 0; for (size_t i = 0; i < k; i++) { double elem = mass[i].first; cnt = 0; for (size_t j = 0; j < k; j++) { cnt += (mass[j].first < elem); //std::cout << cnt; } //if (!AssumeDistinct) for (size_t j = 0; j < i; j++) cnt += (mass[j].first == elem); dstKeys[cnt] = mass[i]; } mass.swap(dstKeys); } void newMerge(std::vector<std::pair<double, size_t>>& currentNN, const std::vector<std::pair<double, size_t>>& candidateNN, std::vector<size_t>& loc, std::vector<size_t>& counter, std::vector<size_t>& offset, std::vector<size_t>& counterScan, std::vector<std::pair<double, size_t>>& updateNN ) { const size_t k = candidateNN.size() / 2; //количество ближайших соседей //std::cout << "loc: " << std::endl; for (size_t j = 0; j < 2 * k; ++j) { loc[j] = BinSearch(currentNN, candidateNN[j].first, 0, (int)k - 1); //std::cout << loc[j] << " "; } //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) counter[j] = j & 1; //std::cout << "counter: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counter[j] << " "; //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) { if ((loc[j] > 0) && ((loc[j] == k) || (candidateNN[j].first == currentNN[loc[j] - 1].first))) //можно через second, проверить!! offset[j] = k; else offset[j] = counter[2 * loc[j]]++; } //std::cout << "counter: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counter[j] << " "; //std::cout << std::endl; //std::cout << "offset: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << offset[j] << " "; //std::cout << std::endl; counterScan[0] = 0; for (size_t j = 1; j < 2 * k; ++j) counterScan[j] = counterScan[j - 1] + counter[j - 1]; //std::cout << "counterScan: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counterScan[j] << " "; //std::cout << std::endl; size_t index; for (size_t j = 0; j < k; ++j) { index = counterScan[2 * j + 1]; //std::cout << index << std::endl; if (index < k) updateNN[index] = currentNN[j]; } //std::cout << "updateNN: " << std::endl; //for (size_t j = 0; j < k; ++j) // std::cout << "(" << updateNN[j].first << ", " << updateNN[j].second << "), "; //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) { if (2 * loc[j] < 2 * k) { index = counterScan[2 * loc[j]] + offset[j]; if (index < k) updateNN[index] = candidateNN[j]; } } currentNN.swap(updateNN); } void mainCycle( const int k, const std::vector<Vortex2D>& vtx, const std::pair<Point2D, Point2D>& gab, const double scale, const int sdvig, std::vector<TParticleCode>& mortonCodeBoth, std::vector<unsigned int>& s, std::vector<TParticleCode>& mortonCodeBoth_temp, std::vector<TParticleCode>& mcdata, std::vector<TParticleCode>& mcquery, std::vector<unsigned int>& iq, std::vector<size_t>& initneig, std::vector<std::vector<std::pair<double, size_t>>>& initdist, double* time ) { //масштабирование координат вихрей в [0;0.75)^2, поиск их мортоновских кодов time[0] = omp_get_wtime(); #pragma omp parallel for for (int i = 0; i < vtx.size(); ++i) { //создание массива data mortonCodeBoth[i].key = Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 }); mortonCodeBoth[i].key <<= 1; //mortonCodeBoth[i].isQuery = 0; mortonCodeBoth[i].originNumber = i; //создание массива query mortonCodeBoth[vtx.size() + i].key = mortonCodeBoth[i].key;// Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 }); mortonCodeBoth[vtx.size() + i].key |= 1; mortonCodeBoth[vtx.size() + i].originNumber = i; //std::cout << "code_data[" << i << "] = " << std::bitset<8>(mortonCodeData[i].key) << std::endl; } time[1] = omp_get_wtime(); //std::cout << "----------" << std::endl; //for (size_t i = 0; i < mortonCodeBoth.size(); ++i) //{ // std::cout << "code[" << i << "] = " << std::bitset<2 * codeLength + 1>(mortonCodeBoth[i].key) << std::endl; //} //сортировка единого массива (data и query) по мортоновским кодам RSort_Parallel(mortonCodeBoth.data(), mortonCodeBoth_temp.data(), (int)mortonCodeBoth.size(), s.data()); time[2] = omp_get_wtime(); //для отладки (рандомно меняем местами data и query с одинаковыми мортоновскими кодами) //for (size_t i = 0; i < mortonCodeBoth.size() / 2; ++i) // if (rand() > RAND_MAX / 2) // std::swap(mortonCodeBoth[2 * i], mortonCodeBoth[2 * i + 1]); //std::cout << "----------" << std::endl; //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные) mcdata.resize(0); mcquery.resize(0); mcdata.reserve(vtx.size()); mcquery.reserve(vtx.size()); time[3] = omp_get_wtime(); iq.resize(0); iq.reserve(vtx.size()); for (size_t i = 0, j = -1; i < mortonCodeBoth.size(); ++i) { if (mortonCodeBoth[i].key & 1) { mcquery.push_back(mortonCodeBoth[i]); iq.push_back((unsigned)(j + 1) == 0 ? 0 : (unsigned int)j); } else { ++j; mcdata.push_back(mortonCodeBoth[i]); } } time[4] = omp_get_wtime(); //std::cout << "----------" << std::endl; // //for (size_t i = 0; i < iq.size(); ++i) // std::cout << "iq[" << i << "] = " << iq[i] << std::endl; // //std::cout << "----------" << std::endl; //for (size_t i = 0; i < mcdata.size(); ++i) // std::cout << "mcdata[" << i << "] = " << std::bitset<8>(mcdata[i].key) << ", isQuery = " << (int)mcdata[i].isQuery << std::endl; //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой) //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей #pragma omp parallel for for (int q = 0; q < mcquery.size(); ++q) { size_t pt = iq[q]; size_t left = pt - k + 1;// , right = pt + k; if (pt < k - 1) { left = 0; //right = 2 * k - 1; } if (pt > mcquery.size() - k - 1) { //right = mcquery.size() - 1; left = (mcquery.size() - 1) - (2 * k - 1); } initneig[q] = left; } time[5] = omp_get_wtime(); //std::cout << "----------" << std::endl; //for (size_t i = 0; i < initneig.size(); ++i) // std::cout << "initneig[" << i << "] = from " << initneig[i] << std::endl; //std::cout << "----------" << std::endl; std::vector<std::pair<double, size_t>> dist(2 * k, { -1.0, -1 }); //std::vector<std::pair<double, size_t>> mergeArray(3 * k); std::vector<size_t> loc(2 * k); std::vector<size_t> counter(2 * k, 0); std::vector<size_t> offset(2 * k, 0); std::vector<size_t> counterScan(2 * k, 0); std::vector<std::pair<double, size_t>> updateNN(k, { 0.0, 0 }); std::vector<std::pair<double, size_t>> dstKeys1(2 * k, { 0.0, 0 }); std::vector<std::pair<double, size_t>> dstKeys(k, { 0.0, 0 }); time[6] = omp_get_wtime(); if (sdvig == 0) #pragma omp parallel for firstprivate(dstKeys1) for (int i = 0; i < initdist.size(); ++i) { for (size_t j = 0; j < 2 * k; ++j) initdist[mcquery[i].originNumber][j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber }; //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; }); newSort(initdist[mcquery[i].originNumber], dstKeys1); initdist[mcquery[i].originNumber].resize(k); } else #pragma omp parallel for firstprivate(dist) firstprivate(loc, counter, offset, counterScan, updateNN, dstKeys)//, firstprivate(mergeArray) for (int i = 0; i < initdist.size(); ++i) { for (size_t j = 0; j < 2 * k; ++j) dist[j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber }; //std::sort(dist.begin(), dist.end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; }); //std::cout << "dist[" << i << "] = "; //for (size_t j = 0; j < k; ++j) // std::cout << "(" << dist[j].first << ", " << dist[j].second << "), "; //std::cout << std::endl; //Mymerge(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), dist.begin(), dist.end(), mergeArray.begin()); newMerge(initdist[mcquery[i].originNumber], dist, loc, counter, offset, counterScan, updateNN); newSort(initdist[mcquery[i].originNumber], dstKeys); //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end()); //std::copy(mergeArray.begin(), mergeArray.begin() + k, initdist[mcquery[i].originNumber].begin()); } time[7] = omp_get_wtime(); } void WakekNN(const std::vector<Vortex2D>& vtx, const size_t k, std::vector<std::vector<std::pair<double, size_t>>>& initdist) { const size_t n = vtx.size(); auto gab = GetGabarit(vtx); double scale = std::max(gab.second[0] - gab.first[0], gab.second[1] - gab.first[1]); //создание единого массива data и query std::vector<TParticleCode> mortonCodeBoth(2 * vtx.size()); //сортировка единого массива (data и query) по мортоновским кодам std::vector<unsigned int> s(256 * omp_get_max_threads()); std::vector<TParticleCode> mortonCodeBoth_temp(2 * vtx.size()); //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные) std::vector<TParticleCode> mcdata, mcquery; std::vector<unsigned int> iq; //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой) //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей std::vector<size_t> initneig(vtx.size(), -1); //выбор из 2*k кандидатов k соседей (выбор по расстоянию от данных до запроса) // initdist = {расстояние; до элемента данных} // сортировка по расстоянию /*std::vector < std::vector < std::pair<double, size_t> >> initdist(vtx.size()); for (auto& d : initdist) d.resize(2 * k, { -1.0, -1 });*/ //Сдвигаем все точки double t1 = omp_get_wtime(); double tStart[5], tFinish[5]; double time[5][8]; for (size_t sdvig = 0; sdvig <= 4; ++sdvig) //4 { tStart[sdvig] = omp_get_wtime(); mainCycle((int)k, vtx, gab, scale, (int)sdvig, mortonCodeBoth, s, mortonCodeBoth_temp, mcdata, mcquery, iq, initneig, initdist, time[sdvig]); tFinish[sdvig] = omp_get_wtime(); //std::cout << "time_knn:" // << " " << time[sdvig][1] - time[sdvig][0] // << " " << time[sdvig][2] - time[sdvig][1] // << " " << time[sdvig][3] - time[sdvig][2] // << " " << time[sdvig][4] - time[sdvig][3] // << " " << time[sdvig][5] - time[sdvig][4] // << " " << time[sdvig][6] - time[sdvig][5] // << " " << time[sdvig][7] - time[sdvig][6] << std::endl; }//for sdvig }
151 
152  if (currentNN[mid].first < x)
153  low = mid + 1;
154  else
155  high = mid - 1;
156  }
157  return mid;
158 }
159 
160 void newSort(std::vector<std::pair<double, size_t>>& mass, std::vector<std::pair<double, size_t>>& dstKeys) {
161  size_t k = mass.size();
162  size_t cnt = 0;
163 
164  for (size_t i = 0; i < k; i++) {
165  double elem = mass[i].first;
166 
167  cnt = 0;
168  for (size_t j = 0; j < k; j++) {
169  cnt += (mass[j].first < elem);
170  //std::cout << cnt;
171  }
172  //if (!AssumeDistinct)
173  for (size_t j = 0; j < i; j++)
174  cnt += (mass[j].first == elem);
175  dstKeys[cnt] = mass[i];
176  }
177  mass.swap(dstKeys);
178 }
179 
180 void newMerge(std::vector<std::pair<double, size_t>>& currentNN, const std::vector<std::pair<double, size_t>>& candidateNN,
181  std::vector<size_t>& loc,
182  std::vector<size_t>& counter,
183  std::vector<size_t>& offset,
184  std::vector<size_t>& counterScan,
185  std::vector<std::pair<double, size_t>>& updateNN
186 )
187 {
188  const size_t k = candidateNN.size() / 2; //количество ближайших соседей //std::cout << "loc: " << std::endl; for (size_t j = 0; j < 2 * k; ++j) { loc[j] = BinSearch(currentNN, candidateNN[j].first, 0, (int)k - 1); //std::cout << loc[j] << " "; } //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) counter[j] = j & 1; //std::cout << "counter: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counter[j] << " "; //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) { if ((loc[j] > 0) && ((loc[j] == k) || (candidateNN[j].first == currentNN[loc[j] - 1].first))) //можно через second, проверить!! offset[j] = k; else offset[j] = counter[2 * loc[j]]++; } //std::cout << "counter: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counter[j] << " "; //std::cout << std::endl; //std::cout << "offset: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << offset[j] << " "; //std::cout << std::endl; counterScan[0] = 0; for (size_t j = 1; j < 2 * k; ++j) counterScan[j] = counterScan[j - 1] + counter[j - 1]; //std::cout << "counterScan: " << std::endl; //for (size_t j = 0; j < 2 * k; ++j) // std::cout << counterScan[j] << " "; //std::cout << std::endl; size_t index; for (size_t j = 0; j < k; ++j) { index = counterScan[2 * j + 1]; //std::cout << index << std::endl; if (index < k) updateNN[index] = currentNN[j]; } //std::cout << "updateNN: " << std::endl; //for (size_t j = 0; j < k; ++j) // std::cout << "(" << updateNN[j].first << ", " << updateNN[j].second << "), "; //std::cout << std::endl; for (size_t j = 0; j < 2 * k; ++j) { if (2 * loc[j] < 2 * k) { index = counterScan[2 * loc[j]] + offset[j]; if (index < k) updateNN[index] = candidateNN[j]; } } currentNN.swap(updateNN); } void mainCycle( const int k, const std::vector<Vortex2D>& vtx, const std::pair<Point2D, Point2D>& gab, const double scale, const int sdvig, std::vector<TParticleCode>& mortonCodeBoth, std::vector<unsigned int>& s, std::vector<TParticleCode>& mortonCodeBoth_temp, std::vector<TParticleCode>& mcdata, std::vector<TParticleCode>& mcquery, std::vector<unsigned int>& iq, std::vector<size_t>& initneig, std::vector<std::vector<std::pair<double, size_t>>>& initdist, double* time ) { //масштабирование координат вихрей в [0;0.75)^2, поиск их мортоновских кодов time[0] = omp_get_wtime(); #pragma omp parallel for for (int i = 0; i < vtx.size(); ++i) { //создание массива data mortonCodeBoth[i].key = Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 }); mortonCodeBoth[i].key <<= 1; //mortonCodeBoth[i].isQuery = 0; mortonCodeBoth[i].originNumber = i; //создание массива query mortonCodeBoth[vtx.size() + i].key = mortonCodeBoth[i].key;// Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 }); mortonCodeBoth[vtx.size() + i].key |= 1; mortonCodeBoth[vtx.size() + i].originNumber = i; //std::cout << "code_data[" << i << "] = " << std::bitset<8>(mortonCodeData[i].key) << std::endl; } time[1] = omp_get_wtime(); //std::cout << "----------" << std::endl; //for (size_t i = 0; i < mortonCodeBoth.size(); ++i) //{ // std::cout << "code[" << i << "] = " << std::bitset<2 * codeLength + 1>(mortonCodeBoth[i].key) << std::endl; //} //сортировка единого массива (data и query) по мортоновским кодам RSort_Parallel(mortonCodeBoth.data(), mortonCodeBoth_temp.data(), (int)mortonCodeBoth.size(), s.data()); time[2] = omp_get_wtime(); //для отладки (рандомно меняем местами data и query с одинаковыми мортоновскими кодами) //for (size_t i = 0; i < mortonCodeBoth.size() / 2; ++i) // if (rand() > RAND_MAX / 2) // std::swap(mortonCodeBoth[2 * i], mortonCodeBoth[2 * i + 1]); //std::cout << "----------" << std::endl; //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные) mcdata.resize(0); mcquery.resize(0); mcdata.reserve(vtx.size()); mcquery.reserve(vtx.size()); time[3] = omp_get_wtime(); iq.resize(0); iq.reserve(vtx.size()); for (size_t i = 0, j = -1; i < mortonCodeBoth.size(); ++i) { if (mortonCodeBoth[i].key & 1) { mcquery.push_back(mortonCodeBoth[i]); iq.push_back((unsigned)(j + 1) == 0 ? 0 : (unsigned int)j); } else { ++j; mcdata.push_back(mortonCodeBoth[i]); } } time[4] = omp_get_wtime(); //std::cout << "----------" << std::endl; // //for (size_t i = 0; i < iq.size(); ++i) // std::cout << "iq[" << i << "] = " << iq[i] << std::endl; // //std::cout << "----------" << std::endl; //for (size_t i = 0; i < mcdata.size(); ++i) // std::cout << "mcdata[" << i << "] = " << std::bitset<8>(mcdata[i].key) << ", isQuery = " << (int)mcdata[i].isQuery << std::endl; //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой) //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей #pragma omp parallel for for (int q = 0; q < mcquery.size(); ++q) { size_t pt = iq[q]; size_t left = pt - k + 1;// , right = pt + k; if (pt < k - 1) { left = 0; //right = 2 * k - 1; } if (pt > mcquery.size() - k - 1) { //right = mcquery.size() - 1; left = (mcquery.size() - 1) - (2 * k - 1); } initneig[q] = left; } time[5] = omp_get_wtime(); //std::cout << "----------" << std::endl; //for (size_t i = 0; i < initneig.size(); ++i) // std::cout << "initneig[" << i << "] = from " << initneig[i] << std::endl; //std::cout << "----------" << std::endl; std::vector<std::pair<double, size_t>> dist(2 * k, { -1.0, -1 }); //std::vector<std::pair<double, size_t>> mergeArray(3 * k); std::vector<size_t> loc(2 * k); std::vector<size_t> counter(2 * k, 0); std::vector<size_t> offset(2 * k, 0); std::vector<size_t> counterScan(2 * k, 0); std::vector<std::pair<double, size_t>> updateNN(k, { 0.0, 0 }); std::vector<std::pair<double, size_t>> dstKeys1(2 * k, { 0.0, 0 }); std::vector<std::pair<double, size_t>> dstKeys(k, { 0.0, 0 }); time[6] = omp_get_wtime(); if (sdvig == 0) #pragma omp parallel for firstprivate(dstKeys1) for (int i = 0; i < initdist.size(); ++i) { for (size_t j = 0; j < 2 * k; ++j) initdist[mcquery[i].originNumber][j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber }; //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; }); newSort(initdist[mcquery[i].originNumber], dstKeys1); initdist[mcquery[i].originNumber].resize(k); } else #pragma omp parallel for firstprivate(dist) firstprivate(loc, counter, offset, counterScan, updateNN, dstKeys)//, firstprivate(mergeArray) for (int i = 0; i < initdist.size(); ++i) { for (size_t j = 0; j < 2 * k; ++j) dist[j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber }; //std::sort(dist.begin(), dist.end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; }); //std::cout << "dist[" << i << "] = "; //for (size_t j = 0; j < k; ++j) // std::cout << "(" << dist[j].first << ", " << dist[j].second << "), "; //std::cout << std::endl; //Mymerge(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), dist.begin(), dist.end(), mergeArray.begin()); newMerge(initdist[mcquery[i].originNumber], dist, loc, counter, offset, counterScan, updateNN); newSort(initdist[mcquery[i].originNumber], dstKeys); //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end()); //std::copy(mergeArray.begin(), mergeArray.begin() + k, initdist[mcquery[i].originNumber].begin()); } time[7] = omp_get_wtime(); } void WakekNN(const std::vector<Vortex2D>& vtx, const size_t k, std::vector<std::vector<std::pair<double, size_t>>>& initdist) { const size_t n = vtx.size(); auto gab = GetGabarit(vtx); double scale = std::max(gab.second[0] - gab.first[0], gab.second[1] - gab.first[1]); //создание единого массива data и query std::vector<TParticleCode> mortonCodeBoth(2 * vtx.size()); //сортировка единого массива (data и query) по мортоновским кодам std::vector<unsigned int> s(256 * omp_get_max_threads()); std::vector<TParticleCode> mortonCodeBoth_temp(2 * vtx.size()); //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные) std::vector<TParticleCode> mcdata, mcquery; std::vector<unsigned int> iq; //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой) //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей std::vector<size_t> initneig(vtx.size(), -1); //выбор из 2*k кандидатов k соседей (выбор по расстоянию от данных до запроса) // initdist = {расстояние; до элемента данных} // сортировка по расстоянию /*std::vector < std::vector < std::pair<double, size_t> >> initdist(vtx.size()); for (auto& d : initdist) d.resize(2 * k, { -1.0, -1 });*/ //Сдвигаем все точки double t1 = omp_get_wtime(); double tStart[5], tFinish[5]; double time[5][8]; for (size_t sdvig = 0; sdvig <= 4; ++sdvig) //4 { tStart[sdvig] = omp_get_wtime(); mainCycle((int)k, vtx, gab, scale, (int)sdvig, mortonCodeBoth, s, mortonCodeBoth_temp, mcdata, mcquery, iq, initneig, initdist, time[sdvig]); tFinish[sdvig] = omp_get_wtime(); //std::cout << "time_knn:" // << " " << time[sdvig][1] - time[sdvig][0] // << " " << time[sdvig][2] - time[sdvig][1] // << " " << time[sdvig][3] - time[sdvig][2] // << " " << time[sdvig][4] - time[sdvig][3] // << " " << time[sdvig][5] - time[sdvig][4] // << " " << time[sdvig][6] - time[sdvig][5] // << " " << time[sdvig][7] - time[sdvig][6] << std::endl; }//for sdvig }
189 
190  //std::cout << "loc: " << std::endl;
191  for (size_t j = 0; j < 2 * k; ++j)
192  {
193  loc[j] = BinSearch(currentNN, candidateNN[j].first, 0, (int)k - 1);
194  //std::cout << loc[j] << " ";
195  }
196  //std::cout << std::endl;
197 
198  for (size_t j = 0; j < 2 * k; ++j)
199  counter[j] = j & 1;
200 
201  //std::cout << "counter: " << std::endl;
202  //for (size_t j = 0; j < 2 * k; ++j)
203  // std::cout << counter[j] << " ";
204  //std::cout << std::endl;
205 
206  for (size_t j = 0; j < 2 * k; ++j)
207  {
208  if ((loc[j] > 0) && ((loc[j] == k) || (candidateNN[j].first == currentNN[loc[j] - 1].first))) //можно через second, проверить!!
209  offset[j] = k;
210  else
211  offset[j] = counter[2 * loc[j]]++;
212  }
213 
214  //std::cout << "counter: " << std::endl;
215  //for (size_t j = 0; j < 2 * k; ++j)
216  // std::cout << counter[j] << " ";
217  //std::cout << std::endl;
218 
219  //std::cout << "offset: " << std::endl;
220  //for (size_t j = 0; j < 2 * k; ++j)
221  // std::cout << offset[j] << " ";
222  //std::cout << std::endl;
223 
224 
225  counterScan[0] = 0;
226  for (size_t j = 1; j < 2 * k; ++j)
227  counterScan[j] = counterScan[j - 1] + counter[j - 1];
228 
229  //std::cout << "counterScan: " << std::endl;
230  //for (size_t j = 0; j < 2 * k; ++j)
231  // std::cout << counterScan[j] << " ";
232  //std::cout << std::endl;
233 
234 
235 
236  size_t index;
237  for (size_t j = 0; j < k; ++j)
238  {
239  index = counterScan[2 * j + 1];
240  //std::cout << index << std::endl;
241  if (index < k)
242  updateNN[index] = currentNN[j];
243  }
244  //std::cout << "updateNN: " << std::endl;
245  //for (size_t j = 0; j < k; ++j)
246  // std::cout << "(" << updateNN[j].first << ", " << updateNN[j].second << "), ";
247  //std::cout << std::endl;
248 
249  for (size_t j = 0; j < 2 * k; ++j)
250  {
251  if (2 * loc[j] < 2 * k)
252  {
253  index = counterScan[2 * loc[j]] + offset[j];
254  if (index < k)
255  updateNN[index] = candidateNN[j];
256  }
257  }
258 
259  currentNN.swap(updateNN);
260 }
261 
262 
263 
264 void mainCycle(
265  const int k,
266  const std::vector<Vortex2D>& vtx, const std::pair<Point2D, Point2D>& gab,
267  const double scale, const int sdvig,
268  std::vector<TParticleCode>& mortonCodeBoth,
269  std::vector<unsigned int>& s, std::vector<TParticleCode>& mortonCodeBoth_temp,
270  std::vector<TParticleCode>& mcdata, std::vector<TParticleCode>& mcquery,
271  std::vector<unsigned int>& iq,
272  std::vector<size_t>& initneig,
273  std::vector<std::vector<std::pair<double, size_t>>>& initdist,
274  double* time
275 )
276 {
277  //масштабирование координат вихрей в [0;0.75)^2, поиск их мортоновских кодов
278 
279  time[0] = omp_get_wtime();
280 #pragma omp parallel for
281  for (int i = 0; i < vtx.size(); ++i)
282  {
283  //создание массива data
284  mortonCodeBoth[i].key = Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 });
285  mortonCodeBoth[i].key <<= 1;
286  //mortonCodeBoth[i].isQuery = 0;
287  mortonCodeBoth[i].originNumber = i;
288 
289  //создание массива query
290  mortonCodeBoth[vtx.size() + i].key = mortonCodeBoth[i].key;// Morton2D((vtx[i].r() - gab.first) * (0.75 / scale) + sdvig * Point2D{ 0.05, 0.05 });
291 
292  mortonCodeBoth[vtx.size() + i].key |= 1;
293  mortonCodeBoth[vtx.size() + i].originNumber = i;
294 
295  //std::cout << "code_data[" << i << "] = " << std::bitset<8>(mortonCodeData[i].key) << std::endl;
296  }
297 
298  time[1] = omp_get_wtime();
299  //std::cout << "----------" << std::endl;
300 
301 
302  //for (size_t i = 0; i < mortonCodeBoth.size(); ++i)
303  //{
304  // std::cout << "code[" << i << "] = " << std::bitset<2 * codeLength + 1>(mortonCodeBoth[i].key) << std::endl;
305  //}
306 
307  //сортировка единого массива (data и query) по мортоновским кодам
308  RSort_Parallel(mortonCodeBoth.data(), mortonCodeBoth_temp.data(), (int)mortonCodeBoth.size(), s.data());
309 
310  time[2] = omp_get_wtime();
311 
312  //для отладки (рандомно меняем местами data и query с одинаковыми мортоновскими кодами)
313  //for (size_t i = 0; i < mortonCodeBoth.size() / 2; ++i)
314  // if (rand() > RAND_MAX / 2)
315  // std::swap(mortonCodeBoth[2 * i], mortonCodeBoth[2 * i + 1]);
316 
317 
318 
319  //std::cout << "----------" << std::endl;
320 
321  //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных
322  //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные)
323  mcdata.resize(0);
324  mcquery.resize(0);
325  mcdata.reserve(vtx.size());
326  mcquery.reserve(vtx.size());
327 
328  time[3] = omp_get_wtime();
329 
330  iq.resize(0);
331  iq.reserve(vtx.size());
332  for (size_t i = 0, j = -1; i < mortonCodeBoth.size(); ++i)
333  {
334  if (mortonCodeBoth[i].key & 1)
335  {
336  mcquery.push_back(mortonCodeBoth[i]);
337  iq.push_back((unsigned)(j + 1) == 0 ? 0 : (unsigned int)j);
338  }
339  else
340  {
341  ++j;
342  mcdata.push_back(mortonCodeBoth[i]);
343  }
344  }
345 
346  time[4] = omp_get_wtime();
347 
348  //std::cout << "----------" << std::endl;
349  //
350  //for (size_t i = 0; i < iq.size(); ++i)
351  // std::cout << "iq[" << i << "] = " << iq[i] << std::endl;
352  //
353  //std::cout << "----------" << std::endl;
354 
355  //for (size_t i = 0; i < mcdata.size(); ++i)
356  // std::cout << "mcdata[" << i << "] = " << std::bitset<8>(mcdata[i].key) << ", isQuery = " << (int)mcdata[i].isQuery << std::endl;
357 
358  //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой)
359  //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей
360 
361 #pragma omp parallel for
362  for (int q = 0; q < mcquery.size(); ++q)
363  {
364  size_t pt = iq[q];
365  size_t left = pt - k + 1;// , right = pt + k;
366  if (pt < k - 1)
367  {
368  left = 0;
369  //right = 2 * k - 1;
370  }
371  if (pt > mcquery.size() - k - 1)
372  {
373  //right = mcquery.size() - 1;
374  left = (mcquery.size() - 1) - (2 * k - 1);
375  }
376 
377  initneig[q] = left;
378  }
379 
380  time[5] = omp_get_wtime();
381 
382  //std::cout << "----------" << std::endl;
383 
384  //for (size_t i = 0; i < initneig.size(); ++i)
385  // std::cout << "initneig[" << i << "] = from " << initneig[i] << std::endl;
386 
387 
388  //std::cout << "----------" << std::endl;
389 
390  std::vector<std::pair<double, size_t>> dist(2 * k, { -1.0, -1 });
391  //std::vector<std::pair<double, size_t>> mergeArray(3 * k);
392 
393  std::vector<size_t> loc(2 * k);
394  std::vector<size_t> counter(2 * k, 0);
395  std::vector<size_t> offset(2 * k, 0);
396  std::vector<size_t> counterScan(2 * k, 0);
397  std::vector<std::pair<double, size_t>> updateNN(k, { 0.0, 0 });
398 
399  std::vector<std::pair<double, size_t>> dstKeys1(2 * k, { 0.0, 0 });
400  std::vector<std::pair<double, size_t>> dstKeys(k, { 0.0, 0 });
401 
402  time[6] = omp_get_wtime();
403 
404  if (sdvig == 0)
405 #pragma omp parallel for firstprivate(dstKeys1)
406  for (int i = 0; i < initdist.size(); ++i)
407  {
408 
409  for (size_t j = 0; j < 2 * k; ++j)
410  initdist[mcquery[i].originNumber][j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber };
411 
412  //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; });
413  newSort(initdist[mcquery[i].originNumber], dstKeys1);
414 
415  initdist[mcquery[i].originNumber].resize(k);
416  }
417  else
418 #pragma omp parallel for firstprivate(dist) firstprivate(loc, counter, offset, counterScan, updateNN, dstKeys)//, firstprivate(mergeArray)
419  for (int i = 0; i < initdist.size(); ++i)
420  {
421  for (size_t j = 0; j < 2 * k; ++j)
422  dist[j] = { (vtx[mcquery[i].originNumber].r() - vtx[mcdata[initneig[i] + j].originNumber].r()).length2(), mcdata[initneig[i] + j].originNumber };
423 
424  //std::sort(dist.begin(), dist.end(), [](const std::pair<double, size_t>& pr1, const std::pair<double, size_t>& pr2) {return pr1.first < pr2.first; });
425 
426  //std::cout << "dist[" << i << "] = ";
427  //for (size_t j = 0; j < k; ++j)
428  // std::cout << "(" << dist[j].first << ", " << dist[j].second << "), ";
429  //std::cout << std::endl;
430 
431  //Mymerge(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end(), dist.begin(), dist.end(), mergeArray.begin());
432  newMerge(initdist[mcquery[i].originNumber], dist, loc, counter, offset, counterScan, updateNN);
433 
434  newSort(initdist[mcquery[i].originNumber], dstKeys);
435 
436  //std::sort(initdist[mcquery[i].originNumber].begin(), initdist[mcquery[i].originNumber].end());
437 
438  //std::copy(mergeArray.begin(), mergeArray.begin() + k, initdist[mcquery[i].originNumber].begin());
439  }
440 
441  time[7] = omp_get_wtime();
442 }
443 
444 
445 
446 
447 void WakekNN(const std::vector<Vortex2D>& vtx, const size_t k, std::vector<std::vector<std::pair<double, size_t>>>& initdist) {
448  const size_t n = vtx.size();
449 
450  auto gab = GetGabarit(vtx);
451  double scale = std::max(gab.second[0] - gab.first[0], gab.second[1] - gab.first[1]);
452 
453  //создание единого массива data и query
454  std::vector<TParticleCode> mortonCodeBoth(2 * vtx.size());
455 
456  //сортировка единого массива (data и query) по мортоновским кодам std::vector<unsigned int> s(256 * omp_get_max_threads()); std::vector<TParticleCode> mortonCodeBoth_temp(2 * vtx.size()); //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные) std::vector<TParticleCode> mcdata, mcquery; std::vector<unsigned int> iq; //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой) //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей std::vector<size_t> initneig(vtx.size(), -1); //выбор из 2*k кандидатов k соседей (выбор по расстоянию от данных до запроса) // initdist = {расстояние; до элемента данных} // сортировка по расстоянию /*std::vector < std::vector < std::pair<double, size_t> >> initdist(vtx.size()); for (auto& d : initdist) d.resize(2 * k, { -1.0, -1 });*/ //Сдвигаем все точки double t1 = omp_get_wtime(); double tStart[5], tFinish[5]; double time[5][8]; for (size_t sdvig = 0; sdvig <= 4; ++sdvig) //4 { tStart[sdvig] = omp_get_wtime(); mainCycle((int)k, vtx, gab, scale, (int)sdvig, mortonCodeBoth, s, mortonCodeBoth_temp, mcdata, mcquery, iq, initneig, initdist, time[sdvig]); tFinish[sdvig] = omp_get_wtime(); //std::cout << "time_knn:" // << " " << time[sdvig][1] - time[sdvig][0] // << " " << time[sdvig][2] - time[sdvig][1] // << " " << time[sdvig][3] - time[sdvig][2] // << " " << time[sdvig][4] - time[sdvig][3] // << " " << time[sdvig][5] - time[sdvig][4] // << " " << time[sdvig][6] - time[sdvig][5] // << " " << time[sdvig][7] - time[sdvig][6] << std::endl; }//for sdvig }
457  std::vector<unsigned int> s(256 * omp_get_max_threads());
458  std::vector<TParticleCode> mortonCodeBoth_temp(2 * vtx.size());
459 
460  //для каждого запроса запоминаем индекс iq -- индекс ближайшего левого элемента данных
461  //отсортиванный массив mortonCodeBoth разделяем на mcquery (только запросы) и mcdata (только данные)
462  std::vector<TParticleCode> mcdata, mcquery;
463 
464  std::vector<unsigned int> iq;
465 
466 
467  //поиск для каждого запроса k точек данных слева и k точек данных справа (при нехватке данных с одной стороны, добавляем с другой)
468  //для каждого запроса находим индекс элемента данных, с которого будут отсчитываться 2*k соседей std::vector<size_t> initneig(vtx.size(), -1); //выбор из 2*k кандидатов k соседей (выбор по расстоянию от данных до запроса) // initdist = {расстояние; до элемента данных} // сортировка по расстоянию /*std::vector < std::vector < std::pair<double, size_t> >> initdist(vtx.size()); for (auto& d : initdist) d.resize(2 * k, { -1.0, -1 });*/ //Сдвигаем все точки double t1 = omp_get_wtime(); double tStart[5], tFinish[5]; double time[5][8]; for (size_t sdvig = 0; sdvig <= 4; ++sdvig) //4 { tStart[sdvig] = omp_get_wtime(); mainCycle((int)k, vtx, gab, scale, (int)sdvig, mortonCodeBoth, s, mortonCodeBoth_temp, mcdata, mcquery, iq, initneig, initdist, time[sdvig]); tFinish[sdvig] = omp_get_wtime(); //std::cout << "time_knn:" // << " " << time[sdvig][1] - time[sdvig][0] // << " " << time[sdvig][2] - time[sdvig][1] // << " " << time[sdvig][3] - time[sdvig][2] // << " " << time[sdvig][4] - time[sdvig][3] // << " " << time[sdvig][5] - time[sdvig][4] // << " " << time[sdvig][6] - time[sdvig][5] // << " " << time[sdvig][7] - time[sdvig][6] << std::endl; }//for sdvig }
469  std::vector<size_t> initneig(vtx.size(), -1);
470 
471  //выбор из 2*k кандидатов k соседей (выбор по расстоянию от данных до запроса)
472  // initdist = {расстояние; до элемента данных}
473  // сортировка по расстоянию
474  /*std::vector < std::vector < std::pair<double, size_t> >> initdist(vtx.size());
475  for (auto& d : initdist)
476  d.resize(2 * k, { -1.0, -1 });*/
477 
478  //Сдвигаем все точки double t1 = omp_get_wtime(); double tStart[5], tFinish[5]; double time[5][8]; for (size_t sdvig = 0; sdvig <= 4; ++sdvig) //4 { tStart[sdvig] = omp_get_wtime(); mainCycle((int)k, vtx, gab, scale, (int)sdvig, mortonCodeBoth, s, mortonCodeBoth_temp, mcdata, mcquery, iq, initneig, initdist, time[sdvig]); tFinish[sdvig] = omp_get_wtime(); //std::cout << "time_knn:" // << " " << time[sdvig][1] - time[sdvig][0] // << " " << time[sdvig][2] - time[sdvig][1] // << " " << time[sdvig][3] - time[sdvig][2] // << " " << time[sdvig][4] - time[sdvig][3] // << " " << time[sdvig][5] - time[sdvig][4] // << " " << time[sdvig][6] - time[sdvig][5] // << " " << time[sdvig][7] - time[sdvig][6] << std::endl; }//for sdvig }
479  double t1 = omp_get_wtime();
480  double tStart[5], tFinish[5];
481  double time[5][8];
482  for (size_t sdvig = 0; sdvig <= 4; ++sdvig) //4
483  {
484 
485  tStart[sdvig] = omp_get_wtime();
486  mainCycle((int)k, vtx, gab, scale, (int)sdvig,
487  mortonCodeBoth,
488  s, mortonCodeBoth_temp, mcdata, mcquery,
489  iq, initneig, initdist, time[sdvig]);
490  tFinish[sdvig] = omp_get_wtime();
491 
492  //std::cout << "time_knn:"
493  // << " " << time[sdvig][1] - time[sdvig][0]
494  // << " " << time[sdvig][2] - time[sdvig][1]
495  // << " " << time[sdvig][3] - time[sdvig][2]
496  // << " " << time[sdvig][4] - time[sdvig][3]
497  // << " " << time[sdvig][5] - time[sdvig][4]
498  // << " " << time[sdvig][6] - time[sdvig][5]
499  // << " " << time[sdvig][7] - time[sdvig][6] << std::endl;
500  }//for sdvig
501 }
VMlib::TParticleCode::key
unsigned int key
Мортоновский код частицы
Definition: Point2D.h:54
ExpandBits
unsigned int ExpandBits(unsigned int v)
Definition: knnCPU.cpp:24
newMerge
void newMerge(std::vector< std::pair< double, size_t >> &currentNN, const std::vector< std::pair< double, size_t >> &candidateNN, std::vector< size_t > &loc, std::vector< size_t > &counter, std::vector< size_t > &offset, std::vector< size_t > &counterScan, std::vector< std::pair< double, size_t >> &updateNN)
Definition: knnCPU.cpp:180
VMlib::Vortex2D::r
HD Point2D & r()
Функция для доступа к радиус-вектору вихря
Definition: Vortex2D.h:85
BinSearch
size_t BinSearch(const std::vector< std::pair< double, size_t >> &currentNN, double x, int low, int high)
Definition: knnCPU.cpp:139
mainCycle
void mainCycle(const int k, const std::vector< Vortex2D > &vtx, const std::pair< Point2D, Point2D > &gab, const double scale, const int sdvig, std::vector< TParticleCode > &mortonCodeBoth, std::vector< unsigned int > &s, std::vector< TParticleCode > &mortonCodeBoth_temp, std::vector< TParticleCode > &mcdata, std::vector< TParticleCode > &mcquery, std::vector< unsigned int > &iq, std::vector< size_t > &initneig, std::vector< std::vector< std::pair< double, size_t >>> &initdist, double *time)
Definition: knnCPU.cpp:264
Morton2D
unsigned int Morton2D(const Point2D &r)
Definition: knnCPU.cpp:56
VMlib::Point2D
Класс, опеделяющий двумерный вектор
Definition: Point2D.h:75
twoPowCodeLength
const int twoPowCodeLength
Definition: knnCPU.cpp:9
VMlib::Vortex2D
Класс, опеделяющий двумерный вихревой элемент
Definition: Vortex2D.h:51
VMlib::dist
R dist(const numvector< T, n > &x, const numvector< P, n > &y)
Вычисление расстояния между двумя точками
Definition: numvector.h:804
RSort_Parallel
void RSort_Parallel(TParticleCode *m, TParticleCode *m_temp, unsigned int n, unsigned int *s)
Сортировка массива из мортоновский кодов
Definition: knnCPU.cpp:66
codeLength
const int codeLength
Definition: knnCPU.cpp:8
GetGabarit
std::pair< Point2D, Point2D > GetGabarit(const std::vector< Vortex2D > &vtx)
Поиск границ габаритного прямоугольника (левый нижний и правый верхний)
Definition: knnCPU.cpp:13
WakekNN
void WakekNN(const std::vector< Vortex2D > &vtx, const size_t k, std::vector< std::vector< std::pair< double, size_t >>> &initdist)
Definition: knnCPU.cpp:447
newSort
void newSort(std::vector< std::pair< double, size_t >> &mass, std::vector< std::pair< double, size_t >> &dstKeys)
Definition: knnCPU.cpp:160