Skip to content

Commit e5f1af2

Browse files
QianruipkuQianruipku
committed
refactor: Add module_pw to cmake
fix: a bug of Chebyshev that it is wrong when using float
1 parent a25b091 commit e5f1af2

File tree

8 files changed

+356
-123
lines changed

8 files changed

+356
-123
lines changed

‎source/module_base/math_chebyshev.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -212,7 +212,7 @@ class FFTW<float>
212212
FFTW(const int norder2_in);
213213
~FFTW();
214214
void execute_fftw();
215-
double* dcoef; //[norder2]
215+
float* dcoef; //[norder2]
216216
fftwf_complex *ccoef;
217217
fftwf_plan coef_plan;
218218
};

‎source/module_base/math_chebyshev_def.h‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -67,7 +67,7 @@ template<typename REAL>
6767
template<class T>
6868
void Chebyshev<REAL>::calcoef_complex(T *ptr, std::complex<REAL> (T::*fun)(std::complex<REAL>))
6969
{
70-
std::complex<double> *pcoef = ( std::complex<double> *)this->fftw.ccoef;
70+
std::complex<REAL> *pcoef = ( std::complex<REAL> *)this->fftw.ccoef;
7171

7272
//three point = 2/3 M + 1/3 T;
7373
//-----------------------------------------------
@@ -80,7 +80,7 @@ void Chebyshev<REAL>::calcoef_complex(T *ptr, std::complex<REAL> (T::*fun)(std::
8080
this->fftw.execute_fftw();
8181
for(int i = 0; i<norder; ++i)
8282
{
83-
double phi=i*ModuleBase::PI/norder2;
83+
REAL phi=i*ModuleBase::PI/norder2;
8484
if(i == 0)
8585
{
8686
coef_complex[i].real( (cos(phi) * pcoef[i].real() + sin(phi) * pcoef[i].imag()) / norder2 * 2 / 3 );
@@ -98,7 +98,7 @@ void Chebyshev<REAL>::calcoef_complex(T *ptr, std::complex<REAL> (T::*fun)(std::
9898
this->fftw.execute_fftw();
9999
for(int i = 0; i<norder; ++i)
100100
{
101-
double phi=i*ModuleBase::PI/norder2;
101+
REAL phi=i*ModuleBase::PI/norder2;
102102
if(i == 0)
103103
{
104104
coef_complex[i].imag( (cos(phi) * pcoef[i].real() + sin(phi) * pcoef[i].imag()) / norder2 * 2 / 3 );
@@ -154,7 +154,7 @@ template<typename REAL>
154154
template<class T>
155155
void Chebyshev<REAL>::calcoef_pair(T *ptr, REAL (T::*fun1)(REAL), REAL (T::*fun2)(REAL))
156156
{
157-
std::complex<double> *pcoef = ( std::complex<double> *)this->fftw.ccoef;
157+
std::complex<REAL> *pcoef = ( std::complex<REAL> *)this->fftw.ccoef;
158158

159159
//three point = 2/3 M + 1/3 T;
160160
//-----------------------------------------------
@@ -167,7 +167,7 @@ void Chebyshev<REAL>::calcoef_pair(T *ptr, REAL (T::*fun1)(REAL), REAL (T::*fun2
167167
this->fftw.execute_fftw();
168168
for(int i = 0; i<norder; ++i)
169169
{
170-
double phi=i*ModuleBase::PI/norder2;
170+
REAL phi=i*ModuleBase::PI/norder2;
171171
if(i == 0)
172172
{
173173
coef_complex[i].real( (cos(phi) * pcoef[i].real() + sin(phi) * pcoef[i].imag()) / norder2 * 2 / 3 );
@@ -185,7 +185,7 @@ void Chebyshev<REAL>::calcoef_pair(T *ptr, REAL (T::*fun1)(REAL), REAL (T::*fun2
185185
this->fftw.execute_fftw();
186186
for(int i = 0; i<norder; ++i)
187187
{
188-
double phi=i*ModuleBase::PI/norder2;
188+
REAL phi=i*ModuleBase::PI/norder2;
189189
if(i == 0)
190190
{
191191
coef_complex[i].imag( (cos(phi) * pcoef[i].real() + sin(phi) * pcoef[i].imag()) / norder2 * 2 / 3 );

‎source/module_base/test/CMakeLists.txt‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -109,6 +109,7 @@ AddTest(
109109
SOURCES container_operator_test.cpp ../container_operator.h
110110
)
111111

112+
#add_definitions(-D__MIX_PRECISION)
112113
AddTest(
113114
TARGET base_math_chebyshev
114115
LIBS ${math_libs}

‎source/module_base/test/math_chebyshev_test.cpp‎

Lines changed: 186 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -45,11 +45,41 @@ class toolfunc{
4545
spin_out[2*i+1] = j*spin_in[2*i];
4646
}
4747
}
48+
#ifdef __MIX_PRECISION
49+
float x7(float x){
50+
return pow(x,7);
51+
}
52+
float x6(float x){
53+
return pow(x,6);
54+
}
55+
float expr(float x){
56+
return exp(x);
57+
}
58+
std::complex<float> expi(std::complex<float> x){
59+
const std::complex<float> j(0.0,1.0);
60+
return exp(j*x);
61+
}
62+
std::complex<float> expi2(std::complex<float> x){
63+
const std::complex<float> j(0.0,1.0);
64+
const float PI = 3.14159265358979323846;
65+
return exp(j*PI/2.0*x);
66+
}
67+
//Pauli matrix: [0,-i;i,0]
68+
void sigma_y(std::complex<float>* spin_in, std::complex<float>* spin_out, const int m = 1){
69+
const std::complex<float> j(0.0,1.0);
70+
for(int i = 0 ; i < m ; ++i)
71+
{
72+
spin_out[2*i] = -j*spin_in[2*i+1];
73+
spin_out[2*i+1] = j*spin_in[2*i];
74+
}
75+
}
76+
#endif
4877
};
4978
class MathChebyshevTest : public testing::Test
5079
{
5180
protected:
5281
ModuleBase::Chebyshev<double> *p_chetest;
82+
ModuleBase::Chebyshev<float> *p_fchetest;
5383
toolfunc fun;
5484
};
5585

@@ -205,4 +235,159 @@ TEST_F(MathChebyshevTest,checkconverge)
205235

206236
delete[] v;
207237
delete p_chetest;
208-
}
238+
}
239+
240+
#ifdef __MIX_PRECISION
241+
TEST_F(MathChebyshevTest,calcoef_real_float)
242+
{
243+
p_fchetest = new ModuleBase::Chebyshev<float>(10);
244+
//x^6 = 1/32*( 10T_0 + 15T_2 + 6T_4 + T_6 )
245+
//x^7 = 1/64*( 35T_1 + 21T_3 + 7T_5 + T_7 )
246+
const float x6ref[10] = {10,0,15,0,6,0,1,0,0,0};
247+
const float x7ref[10] = {0,35,0,21,0,7,0,1,0,0};
248+
p_fchetest->calcoef_real(&fun,&toolfunc::x6);
249+
for(int i = 0;i < 10 ; ++i)
250+
{
251+
EXPECT_NEAR(p_fchetest->coef_real[i]*32.0, x6ref[i] ,1.e-6);
252+
}
253+
p_fchetest->calcoef_real(&fun,&toolfunc::x7);
254+
for(int i = 0;i < 10 ; ++i)
255+
{
256+
EXPECT_NEAR(p_fchetest->coef_real[i]*64.0, x7ref[i] ,1.e-6);
257+
}
258+
delete p_fchetest;
259+
}
260+
261+
TEST_F(MathChebyshevTest,calcoef_pair_float)
262+
{
263+
p_fchetest = new ModuleBase::Chebyshev<float>(10);
264+
//x^6 = 1/32*( 10T_0 + 15T_2 + 6T_4 + T_6 )
265+
//x^7 = 1/64*( 35T_1 + 21T_3 + 7T_5 + T_7 )
266+
const float x6ref[10] = {10,0,15,0,6,0,1,0,0,0};
267+
const float x7ref[10] = {0,35,0,21,0,7,0,1,0,0};
268+
p_fchetest->calcoef_pair(&fun, &toolfunc::x6, &toolfunc::x7);
269+
for(int i = 0;i < 10 ; ++i)
270+
{
271+
EXPECT_NEAR(p_fchetest->coef_complex[i].real()*32.0, x6ref[i] ,1.e-6);
272+
EXPECT_NEAR(p_fchetest->coef_complex[i].imag()*64.0, x7ref[i] ,1.e-6);
273+
}
274+
delete p_fchetest;
275+
}
276+
277+
TEST_F(MathChebyshevTest,calcoef_complex_float)
278+
{
279+
const int norder = 100;
280+
const float PI = 3.14159265358979323846;
281+
p_fchetest = new ModuleBase::Chebyshev<float>(norder);
282+
float *T = new float [norder];
283+
//check exp(i\pi/4) = \sum_n C_n[exp(ix)]T_n(\pi/4) = sqrt(2)/2*(1, i)
284+
p_fchetest->calcoef_complex(&fun, &toolfunc::expi);
285+
p_fchetest->getpolyval(PI/4, T, norder);
286+
std::complex<float> sum(0,0);
287+
for(int i = 0; i < norder ; ++i)
288+
{
289+
sum += p_fchetest->coef_complex[i]*T[i];
290+
}
291+
EXPECT_NEAR(sum.real(), sqrt(2)/2 ,1.e-6);
292+
EXPECT_NEAR(sum.imag(), sqrt(2)/2 ,1.e-6);
293+
delete []T;
294+
delete p_fchetest;
295+
}
296+
297+
TEST_F(MathChebyshevTest,calfinalvec_real_float)
298+
{
299+
const int norder = 100;
300+
const float E = 2.718281828459046;
301+
p_fchetest = new ModuleBase::Chebyshev<float>(norder);
302+
// 1 [ 1/e+e -i(e-1/e) ]
303+
// exp(\sigma_y)= - [ ], where \sigma_y = [0, -i; i, 0]
304+
// 2 [ i(e-1/e) 1/e+e ]
305+
std::complex<float> *v = new std::complex<float> [4];
306+
std::complex<float> *vout = new std::complex<float> [4];
307+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
308+
309+
p_fchetest->calcoef_real(&fun,&toolfunc::expr);
310+
p_fchetest->calfinalvec_real(&fun, &toolfunc::sigma_y, v, vout, 2,2,2);
311+
EXPECT_NEAR(vout[0].real(), 0.5*(E+1/E) ,1.e-6);
312+
EXPECT_NEAR(vout[0].imag(), 0 ,1.e-6);
313+
EXPECT_NEAR(vout[1].real(), 0 ,1.e-6);
314+
EXPECT_NEAR(vout[1].imag(), 0.5*(E-1/E) ,1.e-6);
315+
EXPECT_NEAR(vout[2].real(), 0 ,1.e-6);
316+
EXPECT_NEAR(vout[2].imag(),-0.5*(E-1/E) ,1.e-6);
317+
EXPECT_NEAR(vout[3].real(), 0.5*(E+1/E) ,1.e-6);
318+
EXPECT_NEAR(vout[3].imag(), 0 ,1.e-6);
319+
320+
delete[] v;
321+
delete[] vout;
322+
delete p_fchetest;
323+
}
324+
325+
TEST_F(MathChebyshevTest,calfinalvec_complex_float)
326+
{
327+
const int norder = 100;
328+
const float E = 2.718281828459046;
329+
p_fchetest = new ModuleBase::Chebyshev<float>(norder);
330+
// [ 0 1 ]
331+
// exp(i pi/2*\sigma_y)= [ ], where \sigma_y = [0, -i; i, 0]
332+
// [ -1 0 ]
333+
std::complex<float> *v = new std::complex<float> [4];
334+
std::complex<float> *vout = new std::complex<float> [4];
335+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
336+
337+
p_fchetest->calcoef_complex(&fun,&toolfunc::expi2);
338+
p_fchetest->calfinalvec_complex(&fun, &toolfunc::sigma_y, v, vout, 2,2,2);
339+
EXPECT_NEAR(vout[0].real(), 0 ,1.e-6);
340+
EXPECT_NEAR(vout[0].imag(), 0 ,1.e-6);
341+
EXPECT_NEAR(vout[1].real(),-1 ,1.e-6);
342+
EXPECT_NEAR(vout[1].imag(), 0 ,1.e-6);
343+
EXPECT_NEAR(vout[2].real(), 1 ,1.e-6);
344+
EXPECT_NEAR(vout[2].imag(), 0 ,1.e-6);
345+
EXPECT_NEAR(vout[3].real(), 0 ,1.e-6);
346+
EXPECT_NEAR(vout[3].imag(), 0 ,1.e-6);
347+
348+
delete[] v;
349+
delete[] vout;
350+
delete p_fchetest;
351+
}
352+
353+
TEST_F(MathChebyshevTest,tracepolyA_float)
354+
{
355+
const int norder = 100;
356+
p_fchetest = new ModuleBase::Chebyshev<float>(norder);
357+
358+
std::complex<float> *v = new std::complex<float> [4];
359+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
360+
361+
p_fchetest->tracepolyA(&fun, &toolfunc::sigma_y, v, 2,2,2);
362+
//Trace: even function: 2 ; odd function 0.
363+
for(int i = 0 ; i < norder ; ++i)
364+
{
365+
if(i%2==0) EXPECT_NEAR(p_fchetest->polytrace[i], 2 ,1.e-6);
366+
else EXPECT_NEAR(p_fchetest->polytrace[i], 0 ,1.e-6);
367+
}
368+
369+
delete[] v;
370+
delete p_fchetest;
371+
}
372+
373+
TEST_F(MathChebyshevTest,checkconverge_float)
374+
{
375+
const int norder = 100;
376+
p_fchetest = new ModuleBase::Chebyshev<float>(norder);
377+
378+
std::complex<float> *v = new std::complex<float> [4];
379+
v[0] = 1.0; v[1] = 0.0; v[2] = 0.0; v[3] = 1.0; //[1 0; 0 1]
380+
float tmin = -1.1;
381+
float tmax = 1.1;
382+
bool converge;
383+
converge = p_fchetest->checkconverge(&fun, &toolfunc::sigma_y, v, 2, tmax, tmin, 0.2);
384+
EXPECT_TRUE(converge);
385+
converge = p_fchetest->checkconverge(&fun, &toolfunc::sigma_y, v+2, 2, tmax, tmin, 0.2);
386+
EXPECT_TRUE(converge);
387+
EXPECT_NEAR(tmin, -1.1, 1e-6);
388+
EXPECT_NEAR(tmax, 1.1, 1e-6);
389+
390+
delete[] v;
391+
delete p_fchetest;
392+
}
393+
#endif

‎source/module_esolver/esolver.h‎

Lines changed: 0 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -5,17 +5,9 @@
55
#include "../module_cell/unitcell_pseudo.h"
66
#include "../src_pw/energy.h"
77
#include "../module_base/matrix.h"
8-
//--------------temporary----------------------------
9-
#include "src_lcao/record_adj.h"
10-
#include "src_lcao/local_orbital_charge.h"
11-
#include "src_lcao/local_orbital_wfc.h"
12-
#include "src_lcao/LCAO_hamilt.h"
13-
//--------------\temporary----------------------------
14-
//------It should be moved as fast as possible------
158

169
namespace ModuleESolver
1710
{
18-
1911
class ESolver
2012
{
2113
// protected:

0 commit comments

Comments
 (0)