forked from deepmodeling/abacus-develop
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathsolve_propagation.cpp
More file actions
155 lines (136 loc) · 6.11 KB
/
Copy pathsolve_propagation.cpp
File metadata and controls
155 lines (136 loc) · 6.11 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#include "solve_propagation.h"
#include "source_base/constants.h"
#include "source_base/global_function.h"
#include "source_base/module_external/blas_connector.h"
#include "source_base/module_external/scalapack_connector.h"
#include <iostream>
namespace module_rt
{
#ifdef __MPI
void solve_propagation(const Parallel_Orbitals* pv,
const int nband,
const int nlocal,
const double dt,
const std::complex<double>* Stmp,
const std::complex<double>* Htmp,
const std::complex<double>* psi_k_laststep,
std::complex<double>* psi_k)
{
// (1) init A,B and copy Htmp to A & B
std::complex<double>* operator_A = new std::complex<double>[pv->nloc];
ModuleBase::GlobalFunc::ZEROS(operator_A, pv->nloc);
BlasConnector::copy(pv->nloc, Htmp, 1, operator_A, 1);
std::complex<double>* operator_B = new std::complex<double>[pv->nloc];
ModuleBase::GlobalFunc::ZEROS(operator_B, pv->nloc);
BlasConnector::copy(pv->nloc, Htmp, 1, operator_B, 1);
const double dt_au = dt / ModuleBase::AU_to_FS;
// ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// (2) compute operator_A & operator_B by GEADD
// operator_A = Stmp + i*para * Htmp; beta2 = para = 0.25 * dt
// operator_B = Stmp - i*para * Htmp; beta1 = - para = -0.25 * dt
std::complex<double> alpha = {1.0, 0.0};
std::complex<double> beta1 = {0.0, -0.25 * dt_au};
std::complex<double> beta2 = {0.0, 0.25 * dt_au};
ScalapackConnector::geadd('N', nlocal, nlocal, alpha, Stmp, 1, 1, pv->desc, beta2, operator_A, 1, 1, pv->desc);
ScalapackConnector::geadd('N', nlocal, nlocal, alpha, Stmp, 1, 1, pv->desc, beta1, operator_B, 1, 1, pv->desc);
// ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// (3) b = operator_B @ psi_k_laststep
std::complex<double>* tmp_b = new std::complex<double>[pv->nloc_wfc];
ScalapackConnector::gemm('N',
'N',
nlocal,
nband,
nlocal,
1.0,
operator_B,
1,
1,
pv->desc,
psi_k_laststep,
1,
1,
pv->desc_wfc,
0.0,
tmp_b,
1,
1,
pv->desc_wfc);
//get ipiv
int* ipiv = new int[pv->nloc];
int info = 0;
// (4) solve Ac=b
ScalapackConnector::gesv(nlocal, nband, operator_A, 1, 1, pv->desc, ipiv, tmp_b, 1, 1, pv->desc_wfc, &info);
// copy solution to psi_k
BlasConnector::copy(pv->nloc_wfc, tmp_b, 1, psi_k, 1);
delete[] tmp_b;
delete[] ipiv;
delete[] operator_A;
delete[] operator_B;
}
void solve_propagation(const Parallel_Orbitals* pv,
const int nband,
const int nlocal,
const double dt,
const std::complex<double>* Stmp,
const std::complex<double>* Htmp,
const std::complex<double>* P_k, ///< receives P_k
const std::complex<double>* psi_k_laststep,
std::complex<double>* psi_k)
{
// Print message for debugging, should be removed later
std::cout << "Entering solve_propagation with moving gauge P_k..." << std::endl;
// (1) init A, B and compute HPtmp = Htmp + P_k
std::complex<double>* operator_A = new std::complex<double>[pv->nloc];
std::complex<double>* operator_B = new std::complex<double>[pv->nloc];
// Add up Htmp and P_k to get the effective Hamiltonian matrix for moving spatial gauge
for (int i = 0; i < pv->nloc; ++i)
{
operator_A[i] = Htmp[i] + P_k[i];
operator_B[i] = Htmp[i] + P_k[i];
}
const double dt_au = dt / ModuleBase::AU_to_FS;
// ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// (2) compute operator_A & operator_B by GEADD
// operator_A = Stmp + i*para * (Htmp + P_k);
// operator_B = Stmp - i*para * (Htmp + P_k);
std::complex<double> alpha = {1.0, 0.0};
std::complex<double> beta1 = {0.0, -0.25 * dt_au};
std::complex<double> beta2 = {0.0, 0.25 * dt_au};
ScalapackConnector::geadd('N', nlocal, nlocal, alpha, Stmp, 1, 1, pv->desc, beta2, operator_A, 1, 1, pv->desc);
ScalapackConnector::geadd('N', nlocal, nlocal, alpha, Stmp, 1, 1, pv->desc, beta1, operator_B, 1, 1, pv->desc);
// ->>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
// (3) b = operator_B @ psi_k_laststep
std::complex<double>* tmp_b = new std::complex<double>[pv->nloc_wfc];
ScalapackConnector::gemm('N',
'N',
nlocal,
nband,
nlocal,
1.0,
operator_B,
1,
1,
pv->desc,
psi_k_laststep,
1,
1,
pv->desc_wfc,
0.0,
tmp_b,
1,
1,
pv->desc_wfc);
// get ipiv
int* ipiv = new int[pv->nloc];
int info = 0;
// (4) solve Ac=b
ScalapackConnector::gesv(nlocal, nband, operator_A, 1, 1, pv->desc, ipiv, tmp_b, 1, 1, pv->desc_wfc, &info);
//copy solution to psi_k
BlasConnector::copy(pv->nloc_wfc, tmp_b, 1, psi_k, 1);
delete []tmp_b;
delete []ipiv;
delete []operator_A;
delete []operator_B;
}
#endif // __MPI
} // namespace module_rt