Memory Optimization for Hcontainer Class in Parallel Execution

source/module_hamilt_lcao/module_gint/gint.cpp

@@ -23,7 +23,7 @@
 Gint::~Gint() {
 
     delete this->hRGint;
-    delete this->hRGintCd;
+    delete this->hR_tmp;
     // in gamma_only case, DMRGint.size()=0, 
     // in multi-k case, DMRGint.size()=nspin
     for (int is = 0; is < this->DMRGint.size(); is++) {
@@ -33,7 +33,7 @@ Gint::~Gint() {
         delete this->hRGint_tmp[is];
     }
 #ifdef __MPI
-    delete this->DMRGint_full;
+    delete this->DM2D_tmp;
 #endif
 }
 
@@ -155,11 +155,9 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons
         this->hRGint = new hamilt::HContainer<double>(ucell_in.nat);
     } else {
         npol = 2;
-        if (this->hRGintCd != nullptr) {
-            delete this->hRGintCd;
+        if (this->hR_tmp != nullptr) {
+            delete this->hR_tmp;
         }
-        this->hRGintCd
-            = new hamilt::HContainer<std::complex<double>>(ucell_in.nat);
         for (int is = 0; is < nspin; is++) {
             if (this->DMRGint[is] != nullptr) {
                 delete this->DMRGint[is];
@@ -171,10 +169,9 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons
             this->hRGint_tmp[is] = new hamilt::HContainer<double>(ucell_in.nat);
         }
 #ifdef __MPI
-        if (this->DMRGint_full != nullptr) {
-            delete this->DMRGint_full;
+        if (this->DM2D_tmp != nullptr) {
+            delete this->DM2D_tmp;
         }
-        this->DMRGint_full = new hamilt::HContainer<double>(ucell_in.nat);
 #endif
     }
 
@@ -197,8 +194,10 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons
                                    this->DMRGint[0]->get_memory_size()
                                        * this->DMRGint.size());
     } else {
-        this->hRGintCd->insert_ijrs(this->gridt->get_ijr_info(), ucell_in, npol);
-        this->hRGintCd->allocate(nullptr, true);
+        // this->hRGintCd->insert_ijrs(this->gridt->get_ijr_info(), ucell_in, npol);
+        // this->hRGintCd->allocate(nullptr, true);
+        // ModuleBase::Memory::record("Gint::hRGintCd",
+        //                            this->hRGintCd->get_memory_size());
         for(int is = 0; is < nspin; is++) {
             this->hRGint_tmp[is]->insert_ijrs(this->gridt->get_ijr_info(), ucell_in);
             this->DMRGint[is]->insert_ijrs(this->gridt->get_ijr_info(), ucell_in);
@@ -208,14 +207,7 @@ void Gint::initialize_pvpR(const UnitCell& ucell_in, const Grid_Driver* gd, cons
         ModuleBase::Memory::record("Gint::hRGint_tmp",
                                        this->hRGint_tmp[0]->get_memory_size()*nspin);
         ModuleBase::Memory::record("Gint::DMRGint",
-                                       this->DMRGint[0]->get_memory_size()
-                                           * this->DMRGint.size()*nspin);
-#ifdef __MPI
-        this->DMRGint_full->insert_ijrs(this->gridt->get_ijr_info(), ucell_in, npol);
-        this->DMRGint_full->allocate(nullptr, true);
-        ModuleBase::Memory::record("Gint::DMRGint_full",
-                                   this->DMRGint_full->get_memory_size());
-#endif
+                                       this->DMRGint[0]->get_memory_size()*nspin);
     }
 }
 
@@ -231,10 +223,8 @@ void Gint::reset_DMRGint(const int& nspin)
         {
             for (auto& d : this->DMRGint) { d->allocate(nullptr, false); }
 #ifdef __MPI
-            delete this->DMRGint_full;
-            this->DMRGint_full = new hamilt::HContainer<double>(*this->hRGint);
-            this->DMRGint_full->allocate(nullptr, false);
-#endif
+            delete this->DM2D_tmp;
+#endif       
         }
     }
 }
@@ -262,37 +252,51 @@ void Gint::transfer_DM2DtoGrid(std::vector<hamilt::HContainer<double>*> DM2D) {
     } else // NSPIN=4 case
     {
 #ifdef __MPI
-        hamilt::transferParallels2Serials(*DM2D[0], this->DMRGint_full);
-#else
-        this->DMRGint_full = DM2D[0];
-#endif
-        std::vector<double*> tmp_pointer(4, nullptr);
-        for (int iap = 0; iap < this->DMRGint_full->size_atom_pairs(); ++iap) {
-            auto& ap = this->DMRGint_full->get_atom_pair(iap);
-            int iat1 = ap.get_atom_i();
-            int iat2 = ap.get_atom_j();
-            for (int ir = 0; ir < ap.get_R_size(); ++ir) {
-                const ModuleBase::Vector3<int> r_index = ap.get_R_index(ir);
-                for (int is = 0; is < 4; is++) {
-                    tmp_pointer[is] = this->DMRGint[is]
-                                          ->find_matrix(iat1, iat2, r_index)
-                                          ->get_pointer();
-                }
-                double* data_full = ap.get_pointer(ir);
-                for (int irow = 0; irow < ap.get_row_size(); irow += 2) {
-                    for (int icol = 0; icol < ap.get_col_size(); icol += 2) {
-                        *(tmp_pointer[0])++ = data_full[icol];
-                        *(tmp_pointer[1])++ = data_full[icol + 1];
-                    }
-                    data_full += ap.get_col_size();
-                    for (int icol = 0; icol < ap.get_col_size(); icol += 2) {
-                        *(tmp_pointer[2])++ = data_full[icol];
-                        *(tmp_pointer[3])++ = data_full[icol + 1];
+        // is=0:↑↑, 1:↑↓, 2:↓↑, 3:↓↓
+        const int row_set[4] = {0, 0, 1, 1};
+        const int col_set[4] = {0, 1, 0, 1};
+
+        int mg = DM2D[0]->get_paraV()->get_global_row_size()/2;
+        int ng = DM2D[0]->get_paraV()->get_global_col_size()/2;
+        int nb = DM2D[0]->get_paraV()->get_block_size()/2;
+        int blacs_ctxt = DM2D[0]->get_paraV()->blacs_ctxt;
+        int *iat2iwt = new int[ucell->nat];
+        for (int iat = 0; iat < ucell->nat; iat++) {
+            iat2iwt[iat] = ucell->get_iat2iwt()[iat]/2;
+        }
+        Parallel_Orbitals *pv = new Parallel_Orbitals();
+        pv->set(mg, ng, nb, blacs_ctxt);
+        pv->set_atomic_trace(iat2iwt, ucell->nat, mg);
+        auto ijr_info = DM2D[0]->get_ijr_info();
+        this-> DM2D_tmp = new hamilt::HContainer<double>(pv, nullptr, &ijr_info);
+        ModuleBase::Memory::record("Gint::DM2D_tmp", this->DM2D_tmp->get_memory_size());
+        for (int is = 0; is < 4; is++){
+            this->DM2D_tmp->set_zero();
+            for (int iap = 0; iap < DM2D[0]->size_atom_pairs(); ++iap) {
+                auto& ap = DM2D[0]->get_atom_pair(iap);
+                int iat1 = ap.get_atom_i();
+                int iat2 = ap.get_atom_j();
+                for (int ir = 0; ir < ap.get_R_size(); ++ir) {
+                    const ModuleBase::Vector3<int> r_index = ap.get_R_index(ir);
+                    double* matrix_out = this -> DM2D_tmp -> find_matrix(iat1, iat2, r_index)->get_pointer();
+                    double* matrix_in = ap.get_pointer(ir);
+                    for (int irow = 0; irow < ap.get_row_size()/2; irow ++) {
+                        for (int icol = 0; icol < ap.get_col_size()/2; icol++){
+                            int index_i = irow* ap.get_col_size()/2 + icol;
+                            int index_j = (irow*2+row_set[is]) * ap.get_col_size() + icol*2+col_set[is];
+                            matrix_out[index_i] = matrix_in[index_j];
+                        }
                     }
-                    data_full += ap.get_col_size();
                 }
             }
+            hamilt::transferParallels2Serials( *(this->DM2D_tmp), this->DMRGint[is]);
         }
+        // delete iat2iwt;
+        // delete pv;
+        // delete this-> DM2D_tmp;
+#else
+        //this->DMRGint_full = DM2D[0];
+#endif
     }
     ModuleBase::timer::tick("Gint", "transfer_DMR");
 }

source/module_hamilt_lcao/module_gint/gint.h

@@ -258,13 +258,13 @@ class Gint {
     std::vector<hamilt::HContainer<double>*> hRGint_tmp; 
 
     //! stores Hamiltonian in sparse format
-    hamilt::HContainer<std::complex<double>>* hRGintCd = nullptr; 
+    hamilt::HContainer<std::complex<double>>* hR_tmp = nullptr; 
 
     //! stores DMR in sparse format
     std::vector<hamilt::HContainer<double>*> DMRGint; 
 
     //! tmp tools used in transfer_DM2DtoGrid 
-    hamilt::HContainer<double>* DMRGint_full = nullptr;
+    hamilt::HContainer<double>* DM2D_tmp = nullptr;
 
     std::vector<hamilt::HContainer<double>> pvdpRx_reduced;
     std::vector<hamilt::HContainer<double>> pvdpRy_reduced;

source/module_hamilt_lcao/module_gint/gint_k_pvpr.cpp

@@ -78,84 +78,114 @@ void Gint_k::transfer_pvpR(hamilt::HContainer<std::complex<double>>* hR,
     ModuleBase::TITLE("Gint_k", "transfer_pvpR");
     ModuleBase::timer::tick("Gint_k", "transfer_pvpR");
 
-    this->hRGintCd->set_zero();
-
-    for (int iap = 0; iap < this->hRGintCd->size_atom_pairs(); iap++)
-    {
-        auto* ap = &this->hRGintCd->get_atom_pair(iap);
-        const int iat1 = ap->get_atom_i();
-        const int iat2 = ap->get_atom_j();
-        if (iat1 <= iat2)
-        {
-            hamilt::AtomPair<std::complex<double>>* upper_ap = ap;
-            hamilt::AtomPair<std::complex<double>>* lower_ap = this->hRGintCd->find_pair(iat2, iat1);
-            const hamilt::AtomPair<double>* ap_nspin_0 = this->hRGint_tmp[0]->find_pair(iat1, iat2);
-            const hamilt::AtomPair<double>* ap_nspin_3 = this->hRGint_tmp[3]->find_pair(iat1, iat2);
-            for (int ir = 0; ir < upper_ap->get_R_size(); ir++)
-            {   
-                const auto R_index = upper_ap->get_R_index(ir);
-                auto upper_mat = upper_ap->find_matrix(R_index);
-                auto mat_nspin_0 = ap_nspin_0->find_matrix(R_index);
-                auto mat_nspin_3 = ap_nspin_3->find_matrix(R_index);
+    int mg = hR->get_paraV()->get_global_row_size()/2;
+    int ng = hR->get_paraV()->get_global_col_size()/2;
+    int nb = hR->get_paraV()->get_block_size()/2;
+    hR->set_zero();
+#ifdef __MPI
+    int blacs_ctxt = hR->get_paraV()->blacs_ctxt;
+    int *iat2iwt = new int[ucell_in->nat];
+    for (int iat = 0; iat < ucell_in->nat; iat++) {
+        iat2iwt[iat] = ucell_in->get_iat2iwt()[iat]/2;
+    }
+    Parallel_Orbitals *pv = new Parallel_Orbitals();
+    pv->set(mg, ng, nb, blacs_ctxt);
+    pv->set_atomic_trace(iat2iwt, ucell_in->nat, mg);
+    auto ijr_info = hR->get_ijr_info();
 
-                // The row size and the col size of upper_matrix is double that of matrix_nspin_0
-                for (int irow = 0; irow < mat_nspin_0->get_row_size(); ++irow)
-                {
-                    for (int icol = 0; icol < mat_nspin_0->get_col_size(); ++icol)
-                    {
-                        upper_mat->get_value(2*irow, 2*icol) = mat_nspin_0->get_value(irow, icol) + mat_nspin_3->get_value(irow, icol);
-                        upper_mat->get_value(2*irow+1, 2*icol+1) = mat_nspin_0->get_value(irow, icol) - mat_nspin_3->get_value(irow, icol);
-                    }
-                }
+    this->hR_tmp = new hamilt::HContainer<std::complex<double>>(pv, nullptr, &ijr_info);
+    ModuleBase::Memory::record("Gint::hRGintCd", this->hR_tmp->get_memory_size());
 
-                if (PARAM.globalv.domag)
-                {
-                    const hamilt::AtomPair<double>* ap_nspin_1 = this->hRGint_tmp[1]->find_pair(iat1, iat2);
-                    const hamilt::AtomPair<double>* ap_nspin_2 = this->hRGint_tmp[2]->find_pair(iat1, iat2);
-                    const auto mat_nspin_1 = ap_nspin_1->find_matrix(R_index);
-                    const auto mat_nspin_2 = ap_nspin_2->find_matrix(R_index);
-                    for (int irow = 0; irow < mat_nspin_1->get_row_size(); ++irow)
+    //0,3;1,2;1,2;0,3
+    std::vector<int> first = {0, 1, 1, 0};
+    std::vector<int> second= {3, 2, 2, 3};
+    std::vector<int> row_set = {0, 0, 1, 1};
+    std::vector<int> col_set = {0, 1, 0, 1};
+    std::vector<int> clx_i = {1, 0, 0, -1};
+    std::vector<int> clx_j = {0, 1, -1, 0};
+    for (int is = 0; is < 4; is++){
+        hamilt::HContainer<std::complex<double>>* hRGint_tmpCd = new hamilt::HContainer<std::complex<double>>(this->ucell->nat);
+        hRGint_tmpCd->insert_ijrs(this->gridt->get_ijr_info(), *(this->ucell));
+        hRGint_tmpCd->allocate(nullptr, true);
+        hRGint_tmpCd->set_zero();
+        for (int iap = 0; iap < hRGint_tmpCd->size_atom_pairs(); iap++)
+        {
+            //std::cout<<"iap: "<<iap<<std::endl;
+            auto* ap = &hRGint_tmpCd->get_atom_pair(iap);
+            const int iat1 = ap->get_atom_i();
+            const int iat2 = ap->get_atom_j();
+            if (iat1 <= iat2)
+            {
+                hamilt::AtomPair<std::complex<double>>* upper_ap = ap;
+                hamilt::AtomPair<std::complex<double>>* lower_ap = hRGint_tmpCd->find_pair(iat2, iat1);
+                const hamilt::AtomPair<double>* ap_nspin1 = this->hRGint_tmp[first[is]] ->find_pair(iat1, iat2);
+                const hamilt::AtomPair<double>* ap_nspin2 = this->hRGint_tmp[second[is]] ->find_pair(iat1, iat2);
+                for (int ir = 0; ir < upper_ap->get_R_size(); ir++)
+                {   
+                    //std::cout<<"ir"<<ir<<std::endl;
+                    const auto R_index = upper_ap->get_R_index(ir);
+                    auto upper_mat = upper_ap->find_matrix(R_index);
+                    auto mat_nspin1 = ap_nspin1->find_matrix(R_index);
+                    auto mat_nspin2 = ap_nspin2->find_matrix(R_index);
+                    // The row size and the col size of upper_matrix is double that of matrix_nspin_0
+                    for (int irow = 0; irow < mat_nspin1->get_row_size(); ++irow)
                     {
-                        for (int icol = 0; icol < mat_nspin_1->get_col_size(); ++icol)
+                        for (int icol = 0; icol < mat_nspin1->get_col_size(); ++icol)
+                        {
+                            upper_mat->get_value(irow, icol) = mat_nspin1->get_value(irow, icol) 
+                            + std::complex<double>(clx_i[is], clx_j[is]) * mat_nspin2->get_value(irow, icol);
+                        }
+                    }
+                    //fill the lower triangle matrix
+                    if (PARAM.globalv.domag){
+                        if (iat1 < iat2)
                         {
-                            upper_mat->get_value(2*irow, 2*icol+1) = mat_nspin_1->get_value(irow, icol) +  std::complex<double>(0.0, 1.0) * mat_nspin_2->get_value(irow, icol);
-                            upper_mat->get_value(2*irow+1, 2*icol) = mat_nspin_1->get_value(irow, icol) -  std::complex<double>(0.0, 1.0) * mat_nspin_2->get_value(irow, icol);
+                            auto lower_mat = lower_ap->find_matrix(-R_index);
+                            for (int irow = 0; irow < upper_mat->get_row_size(); ++irow)
+                            {
+                                for (int icol = 0; icol < upper_mat->get_col_size(); ++icol)
+                                {
+                                    lower_mat->get_value(icol, irow) = conj(upper_mat->get_value(irow, icol));
+                                }
+                            }
                         }
                     }
                 }
+            }
+        }
 
-                // fill the lower triangle matrix
-                if (iat1 < iat2)
+        //std::cout<<"success"<<std::endl;
+        this->hR_tmp->set_zero();
+        hamilt::transferSerials2Parallels( *hRGint_tmpCd, this->hR_tmp);
+        for (int iap = 0; iap < hR->size_atom_pairs(); iap++)
+        {
+            //std::cout<<"iap: "<<iap<<std::endl;
+            auto* ap = &hR->get_atom_pair(iap);
+            const int iat1 = ap->get_atom_i();
+            const int iat2 = ap->get_atom_j();
+            auto* ap_nspin = this->hR_tmp ->find_pair(iat1, iat2);
+            for (int ir = 0; ir < ap->get_R_size(); ir++)
+            {   
+                const auto R_index = ap->get_R_index(ir);
+                auto upper_mat = ap->find_matrix(R_index);
+                auto mat_nspin = ap_nspin->find_matrix(R_index);
+
+                // The row size and the col size of upper_matrix is double that of matrix_nspin_0
+                for (int irow = 0; irow < mat_nspin->get_row_size(); ++irow)
                 {
-                    auto lower_mat = lower_ap->find_matrix(-R_index);
-                    for (int irow = 0; irow < upper_mat->get_row_size(); ++irow)
+                    for (int icol = 0; icol < mat_nspin->get_col_size(); ++icol)
                     {
-                        for (int icol = 0; icol < upper_mat->get_col_size(); ++icol)
-                        {
-                            lower_mat->get_value(icol, irow) = conj(upper_mat->get_value(irow, icol));
-                        }
+                        upper_mat->get_value(2*irow+row_set[is], 2*icol+col_set[is]) = 
+                        mat_nspin->get_value(irow, icol);
                     }
                 }
             }
         }
+        delete hRGint_tmpCd;
     }
-
-    // ===================================
-    // transfer HR from Gint to Veff<OperatorLCAO<std::complex<double>, std::complex<double>>>
-    // ===================================
-#ifdef __MPI
-    int size;
-    MPI_Comm_size(MPI_COMM_WORLD, &size);
-    if (size == 1)
-    {
-        hR->add(*this->hRGintCd);
-    }
-    else
-    {
-        hamilt::transferSerials2Parallels<std::complex<double>>(*this->hRGintCd, hR);
-    }
+    delete[] iat2iwt;
 #else
-    hR->add(*this->hRGintCd);
+
 #endif
 
     ModuleBase::timer::tick("Gint_k", "transfer_pvpR");

source/module_lr/utils/gint_move.hpp

@@ -45,8 +45,8 @@ Gint& Gint::operator=(Gint&& rhs)
     // move hR after refactor
     this->hRGint = rhs.hRGint;
     rhs.hRGint = nullptr;
-    this->hRGintCd = rhs.hRGintCd;
-    rhs.hRGintCd = nullptr;
+    this->hR_tmp = rhs.hR_tmp;
+    rhs.hR_tmp = nullptr;
     for (int i = 0; i < this->DMRGint.size(); i++)
     {
         delete this->DMRGint[i];
@@ -60,8 +60,8 @@ Gint& Gint::operator=(Gint&& rhs)
     this->pvdpRz_reduced = std::move(rhs.pvdpRz_reduced);
     this->DMRGint = std::move(rhs.DMRGint);
     this->hRGint_tmp = std::move(rhs.hRGint_tmp);
-    this->DMRGint_full = rhs.DMRGint_full;
-    rhs.DMRGint_full = nullptr;
+    this->DM2D_tmp = rhs.DM2D_tmp;
+    rhs.DM2D_tmp = nullptr;
 
     return *this;
 }