@@ -10,6 +10,14 @@ namespace ModulePW
1010//
1111// Distribute planewaves in reciprocal space to coreors.
1212// Firstly, divide the sphere in reciprocal space into sticks, which are vertical to x-y plane.
13+ //
14+ // Example
15+ // | ---- ixy increasing ---> | ---- ixy increasing ---> |...
16+ // index of sticks 0, 1, 2, ..., nst_per[0]-1, nst_per[0], ..., nst_per[1], ...
17+ // |___________________________|____________________________|___
18+ // ip 0 1 ...
19+ // npw approximate equal to npw approximate equal to...
20+ //
1321// Secondly, distribute these sticks to coreors.
1422// Known: G, GT, GGT, ny, nx, nz, poolnproc, poolrank, ggecut
1523// output: ig2isz[ig], istot2ixy[is], ixy2istot[nxy], is2ixy[is], ixy2ip[ixy], startnsz_per[ip], nst_per[ip], nst
@@ -50,6 +58,8 @@ void PW_Basis::distribution_method1()
5058 std::cout << ' \n ' ;
5159 }
5260 // ------------------------------------------------------------
61+ #ifdef __MPI
62+ // Parallel line
5363
5464 // (2) Collect the x, y indexs, length, bottom of the sticks.
5565 int * st_i = new int [this ->nstot ]; // x or x + nx (if x < 0) of stick.
@@ -62,8 +72,6 @@ void PW_Basis::distribution_method1()
6272 std::cout << " \n The second step done\n " ;
6373 // ------------------------------------------------------------
6474
65- #ifdef __MPI
66- // Parallel line
6775 // (3) Distribute sticks to cores.
6876 int *npw_per = new int [this ->poolnproc ]; // number of planewaves on each core.
6977 int *nst_per = new int [this ->poolnproc ]; // number of sticks on each core.
@@ -79,7 +87,7 @@ void PW_Basis::distribution_method1()
7987 {
8088 this ->ixy2ip [ixy] = -1 ; // meaning this stick has not been distributed or there is no stick on (x, y).
8189 }
82- this ->divide_sticks (tot_npw, st_i, st_j, st_length, npw_per, nst_per);
90+ this ->divide_sticks (st_i, st_j, st_length, npw_per, nst_per);
8391 delete[] st_length;
8492
8593 // for test -----------------------------------------------------------------------------
@@ -182,72 +190,6 @@ void PW_Basis::distribution_method1()
182190 return ;
183191}
184192
185- //
186- // (1) We count the total number of planewaves (tot_npw) and sticks (this->nstot) here.
187- // Meanwhile, we record the number of planewaves on (x, y) in st_length2D, and store the smallest z-coordinate of each stick in st_bottom2D,
188- // so that we can scan a much smaller area in step(2).
189- // known: nx, ny, nz, ggecut, GGT
190- // output: tot_npw, this->nstot, st_length2D, st_bottom2D
191- //
192- void PW_Basis::count_pw_st (
193- int &tot_npw, // total number of planewaves.
194- int * st_length2D, // the number of planewaves that belong to the stick located on (x, y).
195- int * st_bottom2D // the z-coordinate of the bottom of stick on (x, y).
196- )
197- {
198- int ibox[3 ] = {0 , 0 , 0 }; // an auxiliary vector, determine the boundary of the scanning area.
199- ibox[0 ] = int (this ->nx / 2 ) + 1 ; // scan x from -ibox[0] to ibox[0].
200- ibox[1 ] = int (this ->ny / 2 ) + 1 ; // scan y from -ibox[1] to ibox[1].
201- ibox[2 ] = int (this ->nz / 2 ) + 1 ; // scan z from -ibox[2] to ibox[2].
202-
203- ModuleBase::Vector3<double > f;
204-
205- int iy_start = -ibox[1 ]; // determine the scaning area along x-direct, if gamma-only, only positive axis is used.
206- int iy_end = ibox[1 ];
207- if (this ->gamma_only )
208- {
209- iy_start = 0 ;
210- iy_end = this ->ny ;
211- }
212- for (int ix = -ibox[0 ]; ix <= ibox[0 ]; ++ix)
213- {
214- for (int iy = iy_start; iy <= iy_end; ++iy)
215- {
216- // we shift all sticks to the first quadrant in x-y plane here.
217- // (ix, iy, iz) is the direct coordinates of planewaves.
218- // x and y is the coordinates of shifted sticks in x-y plane.
219- // for example, if ny = nx = 10, we will shift the stick on (-1, 2) to (9, 2),
220- // so that its index in st_length and st_bottom is 9 * 10 + 2 = 92.
221- int x = ix;
222- int y = iy;
223- if (x < 0 ) x += this ->nx ;
224- if (y < 0 ) y += this ->ny ;
225- int index = x * this ->ny + y;
226-
227- int length = 0 ; // number of planewave in stick (x, y).
228- for (int iz = -ibox[2 ]; iz <= ibox[2 ]; ++iz)
229- {
230- f.x = ix;
231- f.y = iy;
232- f.z = iz;
233- double modulus = f * (this ->GGT * f);
234- if (modulus <= this ->ggecut )
235- {
236- if (length == 0 ) st_bottom2D[index] = iz; // length == 0 means this point is the bottom of stick (x, y).
237- ++tot_npw;
238- ++length;
239- }
240- }
241- if (length > 0 )
242- {
243- st_length2D[index] = length;
244- ++this ->nstot ;
245- }
246- }
247- }
248- return ;
249- }
250-
251193//
252194// (2) Collect the x, y indexs, length of the sticks.
253195// Firstly, we scan the area and construct temp_st_*.
@@ -363,7 +305,6 @@ void PW_Basis::collect_st(
363305// output: npw_per, nst_per, this->nstnz_per, this->ixy2ip, this->startnsz_per
364306//
365307void PW_Basis::divide_sticks (
366- const int tot_npw, // total number of planewaves.
367308 int * st_i, // x or x + nx (if x < 0) of stick.
368309 int * st_j, // y or y + ny (if y < 0) of stick.
369310 int * st_length, // the stick on (x, y) consists of st_length[x*ny+y] planewaves.
@@ -462,58 +403,4 @@ void PW_Basis::get_istot2ixy(
462403 return ;
463404}
464405
465- //
466- // (5) Construct ig2isz, and is2ixy.
467- // is2ixy contains the x-coordinate and y-coordinate of sticks on current core.
468- // ig2isz contains the z-coordinate of planewaves on current core.
469- // We will scan all the sticks and find the planewaves on them, then store the information into ig2isz and is2ixy.
470- // known: this->nstot, st_bottom2D, st_length2D
471- // output: ig2isz, is2ixy
472- //
473- void PW_Basis::get_ig2isz_is2ixy (
474- int * st_bottom2D, // minimum z of stick, stored in 1d array with this->nstot elements.
475- int * st_length2D // the stick on (x, y) consists of st_length[x*ny+y] planewaves.
476- )
477- {
478- if (this ->npw == 0 )
479- {
480- this ->ig2isz = new int [1 ]; // map ig to the z coordinate of this planewave.
481- this ->ig2isz [0 ] = 0 ;
482- this ->is2ixy = new int [1 ]; // map is (index of sticks) to ixy (ix + iy * nx).
483- this ->is2ixy [0 ] = -1 ;
484- return ;
485- }
486-
487- this ->ig2isz = new int [this ->npw ]; // map ig to the z coordinate of this planewave.
488- ModuleBase::GlobalFunc::ZEROS (this ->ig2isz , this ->npw );
489- this ->is2ixy = new int [this ->nst ]; // map is (index of sticks) to ixy (ix + iy * nx).
490- for (int is = 0 ; is < this ->nst ; ++is)
491- {
492- this ->is2ixy [is] = -1 ;
493- }
494-
495- int st_move = 0 ; // this is the st_move^th stick on current core.
496- int pw_filled = 0 ; // how many current core's planewaves have been found.
497- for (int ixy = 0 ; ixy < this ->nxy ; ++ixy)
498- {
499- if (this ->ixy2ip [ixy] == this ->poolrank )
500- {
501- int zstart = st_bottom2D[ixy];
502- for (int iz = zstart; iz < zstart + st_length2D[ixy]; ++iz)
503- {
504- int z = iz;
505- if (z < 0 ) z += this ->nz ;
506- this ->ig2isz [pw_filled] = st_move * this ->nz + z;
507- pw_filled++;
508- // assert(pw_filled <= this->npw);
509- }
510- this ->is2ixy [st_move] = ixy;
511- st_move++;
512- // assert(st_move <= this->nst);
513- }
514- if (st_move == this ->nst && pw_filled == this ->npw ) break ;
515- }
516- return ;
517- }
518-
519406}
0 commit comments