11#include " gint_tools.h"
22#include " module_base/timer.h"
33#include " module_base/ylm.h"
4+
45namespace Gint_Tools {
6+
57void mult_psi_DMR (
68 const Grid_Technique& gt,
79 const int bxyz,
@@ -19,9 +21,11 @@ void mult_psi_DMR(
1921 const UnitCell& ucell = *gt.ucell ;
2022
2123 // parameters for lapack subroutines
22- constexpr char side = ' L' , uplo = ' U' ;
24+ constexpr char side = ' L' ;
25+ constexpr char uplo = ' U' ;
2326 const char trans = ' N' ;
24- const double alpha = 1.0 , beta = 1.0 ;
27+ const double alpha = 1.0 ;
28+ const double beta = 1.0 ;
2529 const double alpha1 = if_symm ? 2.0 : 1.0 ;
2630
2731 for (int ia1 = 0 ; ia1 < na_grid; ia1++)
@@ -31,37 +35,40 @@ void mult_psi_DMR(
3135 const int T1 = ucell.iat2it [iat1];
3236 const int I1 = ucell.iat2ia [iat1];
3337
34- // ~~~~~~~~~~~~~~~~
35- // get cell R1, this step is redundant in gamma_only case.
36- // ~~~~~~~~~~~~~~~~
38+ // ! get cell R1, this step is redundant in gamma_only case.
3739 const int id1 = gt.which_unitcell [bcell1];
3840 const int R1x = gt.ucell_index2x [id1];
3941 const int R1y = gt.ucell_index2y [id1];
4042 const int R1z = gt.ucell_index2z [id1];
4143
42- if (if_symm) // density
44+ // ! density
45+ if (if_symm)
4346 {
44- // ia2==ia1
47+ // ! ia2==ia1
4548 const auto tmp_matrix = DM->find_matrix (iat1, iat1, 0 , 0 , 0 );
46- // maybe the check "tmp_matrix == nullptr" is not necessary
49+
50+ // ! maybe checking "tmp_matrix == nullptr" is not necessary
4751 if (tmp_matrix == nullptr )
4852 {
4953 continue ;
5054 }
55+
5156 const auto cal_info = Gint_Tools::cal_info (bxyz, ia1, ia1, cal_flag);
5257 const int ib_start = cal_info.first ;
5358 const int ib_len = cal_info.second ;
59+
5460 if (ib_len == 0 )
5561 {
5662 continue ;
5763 }
64+
5865 const auto tmp_matrix_ptr = tmp_matrix->get_pointer ();
5966 const int idx1 = block_index[ia1];
6067 dsymm_ (&side, &uplo, &block_size[ia1], &ib_len, &alpha, tmp_matrix_ptr, &block_size[ia1],
6168 &psi[ib_start][idx1], &LD_pool, &beta, &psi_DMR[ib_start][idx1], &LD_pool);
6269 }
6370
64- // get (j,beta,R2)
71+ // ! get (j,beta,R2)
6572 const int start = if_symm ? ia1 + 1 : 0 ;
6673
6774 for (int ia2 = start; ia2 < na_grid; ia2++)
@@ -70,18 +77,14 @@ void mult_psi_DMR(
7077 const int T2 = ucell.iat2it [gt.which_atom [bcell2]];
7178 const int iat2 = gt.which_atom [bcell2];
7279 const int id2 = gt.which_unitcell [bcell2];
73-
74- // ---------------
75- // get cell R2, this step is redundant in gamma_only case.
76- // ---------------
80+
81+ // ! get cell R2, this step is redundant in gamma_only case.
7782 const int R2x = gt.ucell_index2x [id2];
7883 const int R2y = gt.ucell_index2y [id2];
7984 const int R2z = gt.ucell_index2z [id2];
8085
81- // ------------------------------------------------
82- // calculate the 'offset': R2 position relative
83- // to R1 atom, this step is redundant in gamma_only case.
84- // ------------------------------------------------
86+ // ! calculate the 'offset': R2 position relative
87+ // ! to R1 atom, this step is redundant in gamma_only case.
8588 const int dRx = R1x - R2x;
8689 const int dRy = R1y - R2y;
8790 const int dRz = R1z - R2z;
@@ -103,10 +106,12 @@ void mult_psi_DMR(
103106 }
104107 const int idx1 = block_index[ia1];
105108 const int idx2 = block_index[ia2];
109+
106110 dgemm_ (&trans, &trans, &block_size[ia2], &ib_len, &block_size[ia1], &alpha1, tmp_matrix_ptr, &block_size[ia2],
107111 &psi[ib_start][idx1], &LD_pool, &beta, &psi_DMR[ib_start][idx2], &LD_pool);
108112
109- } // ia2
110- } // ia1
111- }
112- }
113+ } // ia2
114+ } // ia1
115+ }// End of mult_psi_DMR
116+
117+ }// End of Gint_Tools
0 commit comments