1. update annotations about distributeg_method1 and distributeg_method2.

2. correct the heapsort in collect_st in distributeg_method1.cpp.
3. delete the distributeg_method3.cpp.
4. update collect_local_pw in pw_basis.cpp.

Author: sunliang98

source/module_pw/pw_basis.cpp

@@ -98,9 +98,6 @@ void PW_Basis::getstartgr()
 
 //
 // Collect planewaves on current core, and construct gg, gdirect, gcar according to ig2isz and is2ixy.
-// is2ixy contains the x-coordinate and y-coordinate of sticks on current core.
-// ig2isz contains the z-coordinate of planewaves on current core.
-// We will scan the sticks on current core and find the planewaves on them, then store the information into corresponding arrays.
 // known: ig2isz, is2ixy
 // output: gg, gdirect, gcar
 // 
@@ -114,34 +111,24 @@ void PW_Basis::collect_local_pw()
     this->gcar = new ModuleBase::Vector3<double>[this->npw];
 
     ModuleBase::Vector3<double> f;
-    int pw_filled = 0; // how many current core's planewaves have been found.
-    for (int is = 0; is < this->nst; ++is)
+    for(int ig = 0 ; ig < this-> npw ; ++ig)
     {
-        int ix = this->is2ixy[is] / this->ny;
-        int iy = this->is2ixy[is] % this->ny;
+        int isz = this->ig2isz[ig];
+        int iz = isz % this->nz;
+        int is = isz / this->nz;
+        int ixy = this->is2ixy[is];
+        int ix = ixy / this->ny;
+        int iy = ixy % this->ny;
         if (ix >= int(this->nx/2) + 1) ix -= this->nx;
         if (iy >= int(this->bigny/2) + 1) iy -= this->bigny;
-        for (int ig = pw_filled; ig < this->npw; ++ig)
-        {
-            if (this->ig2isz[ig] < (is + 1) * this->nz) // meaning this pw belongs to is^th sticks.
-            {
-                int iz = this->ig2isz[ig] % this->nz;
-                if (iz >= int(this->nz/2) + 1) iz -= this->nz;
-                f.x = ix;
-                f.y = iy;
-                f.z = iz;
-                this->gg[pw_filled] = f * (this->GGT * f);
-                this->gdirect[pw_filled] = f;
-                this->gcar[pw_filled] = f * this->G;
-                pw_filled++;
-            }
-            else
-            {
-                break;
-            }
-        }
+        if (iz >= int(this->nz/2) + 1) iz -= this->nz;
+        f.x = ix;
+        f.y = iy;
+        f.z = iz;
+        this->gg[ig] = f * (this->GGT * f);
+        this->gdirect[ig] = f;
+        this->gcar[ig] = f * this->G;
     }
-    assert(pw_filled == this->npw);
     return;
 }
 

source/module_pw/pw_basis.h

@@ -61,9 +61,9 @@ class PW_Basis
     // on all proc.
     int *ig2isz; // map ig to (is, iz).
     int *istot2bigixy; // istot2bigixy[is]: iy + ix * bigny of is^th stick among all sticks.
-    int *ixy2istot; //useless // ixy2istot[ix + iy * nx]: is of stick on (ix, iy) among all sticks.
-    int *is2ixy; // is2ixy[is]: ix + iy * bignx of is^th stick among sticks on current proc.
-    int *ixy2ip; // useless// ixy2ip[ix + iy * nx]: ip of proc which contains stick on (ix, iy).
+    int *ixy2istot; //useless // ixy2istot[iy + ix * ny]: is of stick on (ix, iy) among all sticks.
+    int *is2ixy; // is2ixy[is]: iy + ix * bigny of is^th stick among sticks on current proc.
+    int *ixy2ip; // useless// ixy2ip[iy + ix * ny]: ip of proc which contains stick on (ix, iy).
     int nst; //num. of sticks in current proc.
     int nstnz; // nst * nz
     int nstot; //num. of sticks in total.
@@ -156,28 +156,6 @@ class PW_Basis
         int* st_length2D,  // the number of planewaves that belong to the stick located on (x, y), stored in 2d x-y plane.
         int* npw_per       // number of planewaves on each core.
     );
-// for distributeg_method3
-    // void divide_sticks2(
-    //     const int tot_npw,  // total number of planewaves.
-    //     int* st_i,          // x or x + nx (if x < 0) of stick.
-    //     int* st_j,          // y or y + ny (if y < 0) of stick.
-    //     int* st_length,     // the stick on (x, y) consists of st_length[x*ny+y] planewaves.
-    //     int* npw_per,       // number of planewaves on each core.
-    //     int* nst_per,       // number of sticks on each core.
-    //     int* is2ip         // ip of core containing is^th stick, map is to ip.         
-    // );
-    // void get_istot2ixy2(
-    //     int* st_i,          // x or x + nx (if x < 0) of stick.
-    //     int* st_j,          // y or y + ny (if y < 0) of stick.
-    //     int* is2ip          // ip of core containing is^th stick, map is to ip.
-    // );
-    // void get_ig2isz_is2ixy2(
-    //     int* st_i,          // x or x + nx (if x < 0) of stick.
-    //     int* st_j,          // y or y + ny (if y < 0) of stick.
-    //     int* st_bottom,     // minimum z of stick, stored in 1d array with tot_nst elements.
-    //     int* st_length,     // the stick on (x, y) consists of st_length[x*ny+y] planewaves.
-    //     int* is2ip          // ip of core containing is^th stick, map is to ip.
-    // );
 
 //===============================================
 //                  FFT

source/module_pw/pw_distributeg.cpp

@@ -41,12 +41,12 @@ void PW_Basis::count_pw_st(
 {
     int ibox[3] = {0, 0, 0};                            // an auxiliary vector, determine the boundary of the scanning area.
     ibox[0] = int(this->nx / 2) + 1;                    // scan x from -ibox[0] to ibox[0].
-    ibox[1] = int(this->ny / 2) + 1;                    // scan y from -ibox[1] to ibox[1].
+    ibox[1] = int(this->ny / 2) + 1;                    // scan y from -ibox[1] to ibox[1], if not gamma-only.
     ibox[2] = int(this->nz / 2) + 1;                    // scan z from -ibox[2] to ibox[2].
 
     ModuleBase::Vector3<double> f;
 
-    int iy_start = -ibox[1]; // determine the scaning area along x-direct, if gamma-only, only positive axis is used.
+    int iy_start = -ibox[1]; // determine the scaning area along y-direct, if gamma-only, only positive axis is used.
     int iy_end = ibox[1];
     if (this->gamma_only)
     {
@@ -69,7 +69,7 @@ void PW_Basis::count_pw_st(
             if (y < 0) y += this->ny;
             int index = x * this->ny + y;
 
-            int length = 0; // number of planewave in stick (x, y).
+            int length = 0; // number of planewave on stick (x, y).
             for (int iz = -ibox[2]; iz <= ibox[2]; ++iz)
             {
                 f.x = ix;
@@ -113,14 +113,14 @@ void PW_Basis::get_ig2isz_is2ixy(
     {
         this->ig2isz = new int[1]; // map ig to the z coordinate of this planewave.
         this->ig2isz[0] = 0;
-        this->is2ixy = new int[1]; // map is (index of sticks) to ixy (ix + iy * nx).
+        this->is2ixy = new int[1]; // map is (index of sticks) to ixy (iy + ix * ny).
         this->is2ixy[0] = -1;
         return;
     }
 
     this->ig2isz = new int[this->npw]; // map ig to the z coordinate of this planewave.
     ModuleBase::GlobalFunc::ZEROS(this->ig2isz, this->npw);
-    this->is2ixy = new int[this->nst]; // map is (index of sticks) to ixy (ix + iy * nx).
+    this->is2ixy = new int[this->nst]; // map is (index of sticks) to ixy (iy + ix * ny).
     for (int is = 0; is < this->nst; ++is) 
     {
         this->is2ixy[is] = -1;

source/module_pw/pw_distributeg_method1.cpp

@@ -9,25 +9,25 @@
 namespace ModulePW
 {
 //
-// Distribute planewaves in reciprocal space to coreors.
+// Distribute planewaves in reciprocal space to cores.
 // Firstly, divide the sphere in reciprocal space into sticks, which are vertical to x-y plane.
+// Secondly, distribute these sticks to coreors.
 // 
 // Example
-//                |  ---- ixy increasing ---> |  ---- ixy increasing --->  |...
-// index of sticks 0, 1, 2, ..., nst_per[0]-1, nst_per[0], ..., nst_per[1], ...
-//                |___________________________|____________________________|___
+//                |  ---- ixy increasing ---> |   ---- ixy increasing --->   |...
+// index of sticks 0, 1, 2, ..., nst_per[0]-1, nst_per[0], ..., nst_per[1]-1, ...
+//                |___________________________|______________________________|___
 // ip                           0                            1              ...
 //                             npw    approximate equal to  npw   approximate equal to...
 // 
-// Secondly, distribute these sticks to coreors.
 //Known: G, GT, GGT, ny, nx, nz, poolnproc, poolrank, ggecut
 //output: ig2isz[ig], istot2bigixy[is], ixy2istot[nxy], is2ixy[is], ixy2ip[ixy], startnsz_per[ip], nst_per[ip], nst
 //
 void PW_Basis::distribution_method1()
 {
     ModuleBase::timer::tick("PW_Basis", "distributeg_method1");
 
-    // initial the variables needed by all proc.
+    // initial the variables needed by all process
     int tot_npw = 0;                     // total number of planewaves.
     this->nstot = 0;                     // total number of sticks.
     int *st_bottom2D = NULL;             // st_bottom2D[ixy], minimum z of stick on (x, y).
@@ -38,10 +38,10 @@ void PW_Basis::distribution_method1()
         // (1) Count the total number of planewaves (tot_npw) and sticks (this->nstot).                  
 
         // Actually we will scan [(2 * ibox[0] + 1) * (2 * ibox[1] + 1)] points on x-y plane,
-        // but we define st_length2D with (ny * nx) points here, because we assume that the diameter
-        // of the sphere is shorter than the sides of the cube.
-        st_length2D = new int[nxy];                    // the number of planewaves that belong to the stick located on (x, y).
-        st_bottom2D = new int[nxy];                    // the z-coordinate of the bottom of stick on (x, y).
+        // but we define st_length2D with (ny * nx) points here, because the diameter
+        // of the sphere should be shorter than the sides of the cube.
+        st_length2D = new int[nxy];
+        st_bottom2D = new int[nxy];
         ModuleBase::GlobalFunc::ZEROS(st_length2D, this->nxy);
         ModuleBase::GlobalFunc::ZEROS(st_bottom2D, this->nxy);
 
@@ -62,7 +62,7 @@ void PW_Basis::distribution_method1()
 #ifdef __MPI
         // Parallel line
 
-        // (2) Collect the x, y indexs, length, bottom of the sticks.
+        // (2) Collect the x, y indexs, and length of the sticks.
         int* st_i = new int[this->nstot];                           // x or x + nx (if x < 0) of stick.
         int* st_j = new int[this->nstot];                           // y or y + ny (if y < 0) of stick.
         int* st_length = new int[this->nstot];                      // number of planewaves in stick.
@@ -73,7 +73,7 @@ void PW_Basis::distribution_method1()
         // std::cout << "\nThe second step done\n";
         // ------------------------------------------------------------
 
-        // (3) Distribute sticks to cores.
+        // (3) Distribute the sticks to cores.
         int *npw_per = new int[this->poolnproc];  // number of planewaves on each core.
         this->nstnz_per = new int[this->poolnproc]; // nz * nst(number of sticks) on each core.
         this->startnsz_per = new int[this->poolnproc];
@@ -82,7 +82,7 @@ void PW_Basis::distribution_method1()
         ModuleBase::GlobalFunc::ZEROS(this->nstnz_per, poolnproc);
         ModuleBase::GlobalFunc::ZEROS(startnsz_per, poolnproc);
 
-        this->ixy2ip = new int[this->nxy];              // ip of core which contains stick on (x, y).
+        this->ixy2ip = new int[this->nxy];              // ip of core which contains the stick on (x, y).
         for (int ixy = 0; ixy < this->nxy; ++ixy)
         {
             this->ixy2ip[ixy] = -1;                 // meaning this stick has not been distributed or there is no stick on (x, y).
@@ -104,7 +104,7 @@ void PW_Basis::distribution_method1()
         // std::cout << "The 3-2 step done\n";
         // --------------------------------------------------------------------------------------
 
-        // (4) Send npw_per, nst[poolrank], st_* to all cores.
+        // (4) Send npw_per[ip], nst_per[ip] to all cores.
         this->npw = npw_per[0];
         this->nst = nst_per[0];
 
@@ -127,7 +127,7 @@ void PW_Basis::distribution_method1()
 
         this->ixy2istot = new int[nxy];
         this->istot2bigixy = new int[this->nstot];
-        this->ixy2ip = new int[nxy];              // ip of core which contains stick on (x, y).
+        this->ixy2ip = new int[nxy];
         int st_move = 0;
         for (int ixy = 0; ixy < nxy; ++ixy)
         {
@@ -140,8 +140,8 @@ void PW_Basis::distribution_method1()
             }
             else
             {
-            ixy2istot[ixy] = -1;
-            ixy2ip[ixy] = -1;
+                ixy2istot[ixy] = -1;
+                ixy2ip[ixy] = -1;
             }
         }
 #endif
@@ -205,7 +205,7 @@ void PW_Basis::collect_st(
     int* st_bottom2D,                               // the z-coordinate of the bottom of stick on (x, y), stored in 2d x-y plane.
     int* st_i,                                      // x or x + nx (if x < 0) of stick.
     int* st_j,                                      // y or y + ny (if y < 0) of stick.
-    int* st_length                                 // number of planewaves in stick, stored in 1d array with this->nstot elements.
+    int* st_length                                  // number of planewaves in stick, stored in 1d array with this->nstot elements.
 )
 {
     int *temp_st_i = new int[this->nstot];                      // x or x + nx (if x < 0) of stick.
@@ -215,7 +215,7 @@ void PW_Basis::collect_st(
 
     int ibox[3] = {0, 0, 0};                            // an auxiliary vector, determine the boundary of the scanning area.
     ibox[0] = int(this->nx / 2) + 1;                    // scan x from -ibox[0] to ibox[0].
-    ibox[1] = int(this->ny / 2) + 1;                    // scan y from -ibox[1] to ibox[1].
+    ibox[1] = int(this->ny / 2) + 1;                    // scan y from -ibox[1] to ibox[1], if not gamma-only.
     ibox[2] = int(this->nz / 2) + 1;                    // scan z from -ibox[2] to ibox[2].
 
     ModuleBase::Vector3<double> f;
@@ -274,13 +274,15 @@ void PW_Basis::collect_st(
 
     int *st_sorted_index = new int[this->nstot]; // indexs in the order of length increasing.
     st_sorted_index[0] = 0;
-    ModuleBase::heapsort(this->nstot, temp_st_length, st_sorted_index); // sort st_* in the order of length decreasing.
+    ModuleBase::heapsort(this->nstot, temp_st_length, st_sorted_index); // sort st_* in the order of length increasing.
 
+    int index = 0;  // indexs in the order of length decreasing.
     for (int istot = 0; istot < this->nstot; ++istot)
     {
-        st_length[istot] = static_cast<int>(temp_st_length[istot]);
-        st_i[istot] = temp_st_i[st_sorted_index[istot]];
-        st_j[istot] = temp_st_j[st_sorted_index[istot]];
+        index = (this->nstot - 1) - istot;
+        st_length[index] = static_cast<int>(temp_st_length[istot]);
+        st_i[index] = temp_st_i[st_sorted_index[istot]];
+        st_j[index] = temp_st_j[st_sorted_index[istot]];
     }
     // std::cout << "st_length    ";
     // for (int is = 0; is < this->nstot; ++is) std::cout << st_length[is] << std::setw(4);
@@ -343,7 +345,7 @@ void PW_Basis::divide_sticks(
                 // }
             }
         }
-        nst_per[ipmin]++;
+        this->nst_per[ipmin]++;
         this->nstnz_per[ipmin] += this->nz;
         npw_per[ipmin] += st_length[is];
         this->ixy2ip[st_i[is] * this->ny + st_j[is]] = ipmin;

source/module_pw/pw_distributeg_method2.cpp

@@ -9,14 +9,14 @@
 namespace ModulePW
 {
 //
-// Distribute planewaves in reciprocal space to coreors.
+// Distribute planewaves in reciprocal space to cores.
 // Firstly, divide the sphere in reciprocal space into sticks, which are vertical to x-y plane.
 // Secondly, distribute these sticks to coreors.
 // 
 // Example
 //                                ---- ixy increasing --->
-// index of sticks 0, 1, 2, ..., nst_per[0]-1, nst_per[0], ..., nst_per[1], ...
-//                |___________________________|____________________________|___
+// index of sticks 0, 1, 2, ..., nst_per[0]-1, nst_per[0], ..., nst_per[1]-1, ...
+//                |___________________________|______________________________|___
 // ip                           0                            1              ...
 // 
 //Known: G, GT, GGT, ny, nx, nz, poolnproc, poolrank, ggecut
@@ -29,7 +29,6 @@ void PW_Basis::distribution_method2()
     // initial the variables needed by all proc.
     int tot_npw = 0;                     // total number of planewaves.
     this->nstot = 0;                     // total number of sticks.
-   // int st_start = 0;                    // index of the first stick on current proc.
     int *st_bottom2D = NULL;             // st_bottom2D[ixy], minimum z of stick on (x, y).
     int *st_length2D = NULL;             // st_length2D[ixy], number of planewaves in stick on (x, y).