Fix: allocate DiagH only once

Author: dyzheng

source/module_hsolver/diagh.h

@@ -4,6 +4,7 @@
 #include "module_base/complexmatrix.h"
 #include "module_hamilt/hamilt.h"
 #include "module_psi/psi.h"
+#include "string"
 
 namespace hsolver
 {
@@ -12,6 +13,7 @@ class DiagH
 {
   public:
     // virtual void init()=0;
+    std::string method = "none";
 
     virtual void diag(hamilt::Hamilt *phm_in, psi::Psi<std::complex<double>> &psi, double *eigenvalue_in) = 0;
 

source/module_hsolver/hsolver_lcao.cpp

@@ -11,11 +11,43 @@ void HSolverLCAO::solveTemplate(hamilt::Hamilt* pHamilt, psi::Psi<T>& psi, elecs
 {
     // select the method of diagonalization
     if (this->method == "genelpa")
-        pdiagh = new DiagoElpa();
+    {
+        if (pdiagh != nullptr)
+        {
+            if (pdiagh->method != this->method)
+            {
+                delete[] pdiagh;
+                pdiagh = new DiagoElpa();
+                pdiagh->method = this->method;
+            }
+        }
+        else
+        {
+            pdiagh = new DiagoElpa();
+            pdiagh->method = this->method;
+        }
+    }
     else if (this->method == "scalapack_gvx")
-        pdiagh = new DiagoBlas();
+    {
+        if (pdiagh != nullptr)
+        {
+            if (pdiagh->method != this->method)
+            {
+                delete[] pdiagh;
+                pdiagh = new DiagoBlas();
+                pdiagh->method = this->method;
+            }
+        }
+        else
+        {
+            pdiagh = new DiagoBlas();
+            pdiagh->method = this->method;
+        }
+    }
     else
+    {
         ModuleBase::WARNING_QUIT("HSolverLCAO::solve", "This method of DiagH is not supported!");
+    }
 
     /// Loop over k points for solve Hamiltonian to charge density
     for (int ik = 0; ik < psi.get_nk(); ++ik)

source/module_hsolver/hsolver_pw.cpp

@@ -27,11 +27,43 @@ void HSolverPW::solve(hamilt::Hamilt* pHamilt, psi::Psi<std::complex<double>>& p
 
     // select the method of diagonalization
     if (this->method == "cg")
-        pdiagh = new DiagoCG(&(GlobalC::hm.hpw), precondition.data());
+    {
+        if(pdiagh!=nullptr)
+        {
+            if(pdiagh->method != this->method)
+            {
+                delete[] pdiagh;
+                pdiagh = new DiagoCG(&(GlobalC::hm.hpw), precondition.data());
+                pdiagh->method = this->method;
+            }
+        }
+        else
+        {
+            pdiagh = new DiagoCG(&(GlobalC::hm.hpw), precondition.data());
+            pdiagh->method = this->method;
+        }
+    }
     else if (this->method == "david")
-        pdiagh = new DiagoDavid(&(GlobalC::hm.hpw), precondition.data());
+    {
+        if (pdiagh != nullptr)
+        {
+            if (pdiagh->method != this->method)
+            {
+                delete[] pdiagh;
+                pdiagh = new DiagoDavid(&(GlobalC::hm.hpw), precondition.data());
+                pdiagh->method = this->method;
+            }
+        }
+        else
+        {
+            pdiagh = new DiagoDavid(&(GlobalC::hm.hpw), precondition.data());
+            pdiagh->method = this->method;
+        }
+    }
     else
+    {
         ModuleBase::WARNING_QUIT("HSolverPW::solve", "This method of DiagH is not supported!");
+    }
 
     /// Loop over k points for solve Hamiltonian to charge density
     for (int ik = 0; ik < psi.get_nk(); ++ik)