Add some timers related to FFT (#6537)

* update small places in charge_mixing_residual.cpp

* update charge_mixing_preconditioner

* add timers and remove some PARAM.inp.nspin

* fix problems

* add timers

* small fix

* fix a potential memory leak

* fix bug

Author: mohanchen

source/source_basis/module_pw/test_gpu/pw_basis_C2C.cpp

@@ -44,9 +44,6 @@ class PW_BASIS_K_GPU_TEST : public ::testing::Test
         int distribution_type = 1;
         bool xprime = false;
         const int nks = 1;
-        ModuleBase::Vector3<double>* kvec_d;
-        kvec_d = new ModuleBase::Vector3<double>[nks];
-        kvec_d[0].set(0, 0, 0);
         // init
         const int mypool = 0;
         const int key = 1;

source/source_cell/module_symmetry/symm_rho.cpp

@@ -9,6 +9,10 @@ void Symmetry::rho_symmetry( double *rho,
 {
     ModuleBase::timer::tick("Symmetry","rho_symmetry");
 
+    assert(nr1>0);
+    assert(nr2>0);
+    assert(nr3>0);
+
 	// allocate flag for each FFT grid.
     bool* symflag = new bool[nr1 * nr2 * nr3];
     for (int i=0; i<nr1*nr2*nr3; i++)

source/source_esolver/esolver_of_tool.cpp

@@ -220,7 +220,7 @@ double ESolver_OF::cal_mu(double* pphi, double* pdEdphi, double nelec)
  * @brief Rotate and renormalize the direction |d>,
  * make it orthogonal to phi (<d|phi> = 0), and <d|d> = nelec
  */
-void ESolver_OF::adjust_direction()
+void ESolver_OF::adjust_direction(void)
 {
     // filter the high frequency term in direction if of_full_pw = false
     if (!PARAM.inp.of_full_pw)

source/source_estate/module_charge/charge_mixing_preconditioner.cpp

@@ -6,21 +6,32 @@
 
 void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
 {
-    if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) {
-        return;
-}
+    ModuleBase::TITLE("Charge_Mixing", "Kerker_screen_recip");
+
+	if (this->mixing_gg0 <= 0.0 || this->mixing_beta <= 0.1) 
+	{
+		return;
+	}
+
+    ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_recip");
+
+    const int nspin = PARAM.inp.nspin;
+
     double fac = 0.0;
     double gg0 = 0.0;
     double amin = 0.0;
 
     /// consider a resize for mixing_angle
-    int resize_tmp = 1;
-    if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { resize_tmp = 2;
-}
+	int resize_tmp = 1;
+    if (nspin == 4 && this->mixing_angle > 0) 
+    { 
+    	resize_tmp = 2;
+    }
 
     /// implement Kerker for density and magnetization separately
-    for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is)
+    for (int is = 0; is < nspin / resize_tmp; ++is)
     {
+        const int is_idx = is * this->rhopw->npw;
         /// new mixing method only support nspin=2 not nspin=4
         if (is >= 1)
         {
@@ -29,10 +40,10 @@ void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
 #ifdef __DEBUG
                 assert(is == 1); // make sure break works
 #endif
-                double is_mag = PARAM.inp.nspin - 1;
+                double is_mag = nspin - 1;
                 //for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
                 //{
-                //    drhog[is * this->rhopw->npw + ig] *= 1;
+                //    drhog[is_idx + ig] *= 1;
                 //}
                 break;
             }
@@ -46,32 +57,49 @@ void Charge_Mixing::Kerker_screen_recip(std::complex<double>* drhog)
         }
 
         gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2);
+
+        const double gg0_amin = this->mixing_gg0_min / amin;
+
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 512)
 #endif
         for (int ig = 0; ig < this->rhopw->npw; ++ig)
         {
             double gg = this->rhopw->gg[ig];
-            double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin);
-            drhog[is * this->rhopw->npw + ig] *= filter_g;
+            double filter_g = std::max(gg / (gg + gg0), gg0_amin);
+            drhog[is_idx + ig] *= filter_g;
         }
     }
+
+    ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_recip");
     return;
 }
 
 void Charge_Mixing::Kerker_screen_real(double* drhor)
 {
-    if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) {
-        return;
-}
-    /// consider a resize for mixing_angle
+    ModuleBase::TITLE("Charge_Mixing", "Kerker_screen_real");
+
+	if (this->mixing_gg0 <= 0.0001 || this->mixing_beta <= 0.1) 
+	{
+		return;
+	}
+
+    ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_real");
+
+    const int nspin = PARAM.inp.nspin;
+    assert(nspin==1 || nspin==2 || nspin==4);
+
+	/// consider a resize for mixing_angle
     int resize_tmp = 1;
-    if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { resize_tmp = 2;
-}
+	if (nspin == 4 && this->mixing_angle > 0) 
+	{ 
+		resize_tmp = 2;
+	}
 
-    std::vector<std::complex<double>> drhog(this->rhopw->npw * PARAM.inp.nspin / resize_tmp);
-    std::vector<double> drhor_filter(this->rhopw->nrxx * PARAM.inp.nspin / resize_tmp);
-    for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is)
+    std::vector<std::complex<double>> drhog(this->rhopw->npw * nspin / resize_tmp);
+    std::vector<double> drhor_filter(this->rhopw->nrxx * nspin / resize_tmp);
+
+    for (int is = 0; is < nspin / resize_tmp; ++is)
     {
         // Note after this process some G which is higher than Gmax will be filtered.
         // Thus we cannot use Kerker_screen_recip(drhog.data()) directly after it.
@@ -82,7 +110,7 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
     double gg0 = 0.0;
     double amin = 0.0;
 
-    for (int is = 0; is < PARAM.inp.nspin / resize_tmp; is++)
+    for (int is = 0; is < nspin / resize_tmp; is++)
     {
 
         if (is >= 1)
@@ -92,8 +120,8 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
 #ifdef __DEBUG
                 assert(is == 1); /// make sure break works
 #endif
-                double is_mag = PARAM.inp.nspin - 1;
-                if (PARAM.inp.nspin == 4 && this->mixing_angle > 0) { is_mag = 1;
+                double is_mag = nspin - 1;
+                if (nspin == 4 && this->mixing_angle > 0) { is_mag = 1;
 }
                 for (int ig = 0; ig < this->rhopw->npw * is_mag; ig++)
                 {
@@ -111,6 +139,9 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
         }
 
         gg0 = std::pow(fac * 0.529177 / *this->tpiba, 2);
+
+        const int is_idx = is * this->rhopw->npw;
+        const double gg0_amin = this->mixing_gg0_min / amin;
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 512)
 #endif
@@ -120,24 +151,27 @@ void Charge_Mixing::Kerker_screen_real(double* drhor)
             // I have not decided how to handle gg=0 part, will be changed in future
             //if (gg == 0)
             //{
-            //    drhog[is * this->rhopw->npw + ig] *= 0;
+            //    drhog[is_idx + ig] *= 0;
             //    continue;
             //}
-            double filter_g = std::max(gg / (gg + gg0), this->mixing_gg0_min / amin);
-            drhog[is * this->rhopw->npw + ig] *= (1 - filter_g);
+            double filter_g = std::max(gg / (gg + gg0), gg0_amin);
+            drhog[is_idx + ig] *= (1 - filter_g);
         }
     }
     /// inverse FT
-    for (int is = 0; is < PARAM.inp.nspin / resize_tmp; ++is)
+    for (int is = 0; is < nspin / resize_tmp; ++is)
     {
         this->rhopw->recip2real(drhog.data() + is * this->rhopw->npw, drhor_filter.data() + is * this->rhopw->nrxx);
     }
 
 #ifdef _OPENMP
 #pragma omp parallel for schedule(static, 512)
 #endif
-    for (int ir = 0; ir < this->rhopw->nrxx * PARAM.inp.nspin / resize_tmp; ir++)
+    for (int ir = 0; ir < this->rhopw->nrxx * nspin / resize_tmp; ir++)
     {
         drhor[ir] -= drhor_filter[ir];
     }
+
+    ModuleBase::timer::tick("Charge_Mixing", "Kerker_screen_real");
+    return;
 }

source/source_estate/module_charge/charge_mixing_residual.cpp

@@ -9,11 +9,13 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
 {
     ModuleBase::TITLE("Charge_Mixing", "get_drho");
     ModuleBase::timer::tick("Charge_Mixing", "get_drho");
+    const int nspin = PARAM.inp.nspin;
+    assert(nspin==1 || nspin==2 || nspin==4);
     double drho = 0.0;
 
     if (PARAM.inp.scf_thr_type == 1)
     {
-        for (int is = 0; is < PARAM.inp.nspin; ++is)
+        for (int is = 0; is < nspin; ++is)
         {
             ModuleBase::GlobalFunc::NOTE("Perform FFT on rho(r) to obtain rho(G).");
             chr->rhopw->real2recip(chr->rho[is], chr->rhog[is]);
@@ -23,15 +25,15 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
         }
 
         ModuleBase::GlobalFunc::NOTE("Calculate the charge difference between rho(G) and rho_save(G)");
-        std::vector<std::complex<double>> drhog(PARAM.inp.nspin * this->rhopw->npw);
+        std::vector<std::complex<double>> drhog(nspin * this->rhopw->npw);
 #ifdef _OPENMP
 #pragma omp parallel for collapse(2) schedule(static, 512)
 #endif
-        for (int is = 0; is < PARAM.inp.nspin; ++is)
+        for (int is = 0; is < nspin; ++is)
         {
             for (int ig = 0; ig < this->rhopw->npw; ig++)
             {
-                drhog[is * rhopw->npw + ig] = chr->rhog[is][ig] - chr->rhog_save[is][ig];
+                drhog[is * this->rhopw->npw + ig] = chr->rhog[is][ig] - chr->rhog_save[is][ig];
             }
         }
 
@@ -42,7 +44,7 @@ double Charge_Mixing::get_drho(Charge* chr, const double nelec)
     {
         // Note: Maybe it is wrong.
         //       The inner_product_real function (L1-norm) is different from that (L2-norm) in mixing.
-        for (int is = 0; is < PARAM.inp.nspin; is++)
+        for (int is = 0; is < nspin; is++)
         {
             if (is != 0 && is != 3 && PARAM.globalv.domag_z)
             {