1+ #include " ./math_chebyshev.h"
2+ #include " ./constants.h"
3+ #include " ./blas_connector.h"
4+ #include " ./global_function.h"
5+ namespace ModuleBase
6+ {
7+ // we only have two examples: double and float.
8+ template class Chebyshev <double >;
9+ #ifdef __MIX_PRECISION
10+ template class Chebyshev <float >;
11+ #endif
12+
13+ FFTW<double >::FFTW(const int norder2_in)
14+ {
15+ ccoef = (fftw_complex *) fftw_malloc (sizeof (fftw_complex) * norder2_in);
16+ dcoef = (double *) fftw_malloc (sizeof (double ) * norder2_in);
17+ coef_plan = fftw_plan_dft_r2c_1d (norder2_in, dcoef, ccoef, FFTW_ESTIMATE);
18+ }
19+ FFTW<double >::~FFTW ()
20+ {
21+ fftw_destroy_plan (coef_plan);
22+ fftw_free (ccoef);
23+ fftw_free (dcoef);
24+ }
25+ void FFTW<double >::execute_fftw()
26+ {
27+ fftw_execute (this ->coef_plan );
28+ }
29+
30+ #ifdef __MIX_PRECISION
31+ FFTW<float >::FFTW(const int norder2_in)
32+ {
33+ ccoef = (fftwf_complex *) fftw_malloc (sizeof (fftwf_complex) * norder2_in);
34+ dcoef = (float *) fftw_malloc (sizeof (float ) * norder2_in);
35+ coef_plan = fftwf_plan_dft_r2c_1d (norder2_in, dcoef, ccoef, FFTW_ESTIMATE);
36+ }
37+ FFTW<float >::~FFTW ()
38+ {
39+ fftwf_destroy_plan (coef_plan);
40+ fftw_free (ccoef);
41+ fftw_free (dcoef);
42+ }
43+ void FFTW<float >::execute_fftw()
44+ {
45+ fftwf_execute (this ->coef_plan );
46+ }
47+ #endif
48+
49+ // A number to control the number of grids in C_n integration
50+ #define EXTEND 16
51+
52+ template <typename REAL>
53+ Chebyshev<REAL>::Chebyshev(const int norder_in) : fftw(2 * EXTEND * norder_in)
54+ {
55+ this ->norder = norder_in;
56+ norder2 = 2 * norder * EXTEND;
57+ if (this ->norder < 1 )
58+ {
59+ ModuleBase::WARNING_QUIT (" Stochastic_Chebychev" " The Chebyshev expansion order should be at least 1!"
60+ }
61+ polytrace = new REAL [norder];
62+ coef_real = new REAL [norder];
63+ coef_complex = new std::complex <REAL> [norder];
64+
65+ // ndmin = ndmax = ndmax_in;
66+
67+ getcoef_complex = false ;
68+ getcoef_real = false ;
69+ }
70+
71+ template <typename REAL>
72+ Chebyshev<REAL>::~Chebyshev ()
73+ {
74+ delete [] polytrace;
75+ delete [] coef_real;
76+ delete [] coef_complex;
77+ }
78+
79+ template <typename REAL>
80+ void Chebyshev<REAL>::getpolyval(const REAL x, REAL* polyval, const int N)
81+ {
82+ polyval[0 ] = 1 ;
83+ polyval[1 ] = x;
84+ for (int i = 2 ; i < N; ++i)
85+ {
86+ polyval[i] = 2 * x * polyval[i-1 ] - polyval[i-2 ];
87+ }
88+ }
89+ template <typename REAL>
90+ REAL Chebyshev<REAL>::recurs(const REAL x, const REAL Tn, REAL const Tn_1)
91+ {
92+ return 2 *x*Tn-Tn_1;
93+ }
94+
95+ template <typename REAL>
96+ REAL Chebyshev<REAL>:: ddot_real(
97+ const std::complex <REAL>* psi_L,
98+ const std::complex <REAL>* psi_R,
99+ const int N, const int LDA, const int m)
100+ {
101+ REAL result = 0 ;
102+ if (N == LDA || m==1 )
103+ {
104+ int dim2=2 * N * m;
105+ REAL *pL,*pR;
106+ pL=(REAL *)psi_L;
107+ pR=(REAL *)psi_R;
108+ result=BlasConnector::dot (dim2,pL,1 ,pR,1 );
109+ }
110+ else
111+ {
112+ REAL *pL,*pR;
113+ pL=(REAL *)psi_L;
114+ pR=(REAL *)psi_R;
115+ for (int i = 0 ; i < m ; ++i)
116+ {
117+ int dim2=2 * N;
118+ result += BlasConnector::dot (dim2,pL,1 ,pR,1 );
119+ pL += 2 * LDA;
120+ pR += 2 * LDA;
121+ }
122+ }
123+ return result;
124+ }
125+ }
0 commit comments