forked from abacusmodeling/abacus-develop
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathdiago_lapack.cpp
More file actions
83 lines (71 loc) · 2.47 KB
/
diago_lapack.cpp
File metadata and controls
83 lines (71 loc) · 2.47 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
#include "diago_lapack.h"
#include "../src_pw/global.h"
#include "module_base/global_variable.h"
#include "module_base/lapack_connector.h"
#include "module_base/timer.h"
#include "module_base/tool_quit.h"
typedef hamilt::MatrixBlock<double> matd;
typedef hamilt::MatrixBlock<std::complex<double>> matcd;
namespace hsolver
{
void DiagoLapack::diag(hamilt::Hamilt *phm_in, psi::Psi<std::complex<double>> &psi, double *eigenvalue_in)
{
ModuleBase::TITLE("DiagoLapack", "diag");
assert(GlobalV::NPROC == 1);
ModuleBase::ComplexMatrix Htmp(GlobalV::NLOCAL, GlobalV::NLOCAL);
ModuleBase::ComplexMatrix Stmp(GlobalV::NLOCAL, GlobalV::NLOCAL);
matcd h_mat, s_mat;
phm_in->matrix(h_mat, s_mat);
for (int i = 0; i < GlobalV::NLOCAL; i++)
{
for (int j = 0; j < GlobalV::NLOCAL; j++)
{
Htmp(i, j) = h_mat.p[i * GlobalV::NLOCAL + j];
Stmp(i, j) = s_mat.p[i * GlobalV::NLOCAL + j];
}
}
//----------------------------
// keep this for tests
// out.printcm_norm("Lapack_H", Htmp, 1.0e-5);
// out.printcm_norm("Lapack_S", Stmp, 1.0e-5);
//----------------------------
double *en = new double[GlobalV::NLOCAL];
ModuleBase::GlobalFunc::ZEROS(en, GlobalV::NLOCAL);
ModuleBase::ComplexMatrix hvec(GlobalV::NLOCAL, GlobalV::NBANDS);
GlobalC::hm.diagH_LAPACK(GlobalV::NLOCAL, GlobalV::NBANDS, Htmp, Stmp, GlobalV::NLOCAL, en, hvec);
if (GlobalV::NSPIN != 4)
{
for (int ib = 0; ib < GlobalV::NBANDS; ib++)
{
for (int iw = 0; iw < GlobalV::NLOCAL; iw++)
{
// wfc_k_grid[ib][iw] = hvec(iw, ib);
// wfc_k.c[ib * GlobalV::NLOCAL + iw] = wfc_k_grid[ib][iw];
psi.get_pointer()[ib * GlobalV::NLOCAL + iw] = hvec(iw, ib);
}
}
}
/*
else
{
for (int ib = 0; ib < GlobalV::NBANDS; ib++)
{
for (int iw = 0; iw < GlobalV::NLOCAL / GlobalV::NPOL; iw++)
{
wfc_k_grid[ib][iw] = hvec(iw * GlobalV::NPOL, ib);
wfc_k_grid[ib][iw + GlobalV::NLOCAL / GlobalV::NPOL] = hvec(iw * GlobalV::NPOL + 1, ib);
}
}
}
*/
// energy for k-point ik
for (int ib = 0; ib < GlobalV::NBANDS; ib++)
{
eigenvalue_in[ib] = en[ib];
}
}
void DiagoLapack::diag(hamilt::Hamilt *phm_in, psi::Psi<double> &psi, double *eigenvalue_in)
{
ModuleBase::TITLE("DiagoLapack", "diag");
}
} // namespace hsolver