1 //////////////////////////////////////////////////////////////////////////////////////
2 // This file is distributed under the University of Illinois/NCSA Open Source License.
3 // See LICENSE file in top directory for details.
4 //
5 // Copyright (c) 2016 Jeongnim Kim and QMCPACK developers.
6 //
7 // File developed by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
8 // Jeremy McMinnis, jmcminis@gmail.com, University of Illinois at Urbana-Champaign
9 // Mark A. Berrill, berrillma@ornl.gov, Oak Ridge National Laboratory
10 //
11 // File created by: Jeongnim Kim, jeongnim.kim@gmail.com, University of Illinois at Urbana-Champaign
12 //////////////////////////////////////////////////////////////////////////////////////
13
14
15 #ifndef BSPLINE_ONEDIM_SOLVERS
16 #define BSPLINE_ONEDIM_SOLVERS
17 #include <vector>
18 #include <complex>
19 /** dummy declaration to be specialized
20 *
21 * Specializations for double and std::complex<double>
22 */
23 template<typename T>
24 struct SolvePeriodicInterp1D
25 {};
26
27 template<typename T>
FuncSolvePeriodicInterp1D(const std::vector<T> & data,std::vector<T> & p)28 inline void FuncSolvePeriodicInterp1D(const std::vector<T>& data, std::vector<T>& p)
29 {
30 const T ratio = 0.25;
31 int N = data.size();
32 std::vector<T> d(N), gamma(N), mu(N);
33 //d = 1.5*data;
34 for (int i = 0; i < N; i++)
35 d[i] = 1.5 * data[i];
36 // First, eliminate leading coefficients
37 gamma[0] = ratio;
38 mu[0] = ratio;
39 mu[N - 1] = ratio;
40 gamma[N - 1] = 1.0;
41 for (int row = 1; row < (N - 1); row++)
42 {
43 double diag = 1.0 - mu[row - 1] * ratio;
44 double diagInv = 1.0 / diag;
45 gamma[row] = -ratio * gamma[row - 1] * diagInv;
46 mu[row] = diagInv * ratio;
47 d[row] = diagInv * (d[row] - ratio * d[row - 1]);
48 // Last row
49 d[N - 1] -= mu[N - 1] * d[row - 1];
50 gamma[N - 1] -= mu[N - 1] * gamma[row - 1];
51 mu[N - 1] = -mu[N - 1] * mu[row - 1];
52 }
53 // Last row: gamma(N-1) hold diagonal element
54 mu[N - 1] += ratio;
55 gamma[N - 1] -= mu[N - 1] * (mu[N - 2] + gamma[N - 2]);
56 d[N - 1] -= mu[N - 1] * d[N - 2];
57 p[N - 1] = d[N - 1] / gamma[N - 1];
58 // Now go back upward, back substituting
59 for (int row = N - 2; row >= 0; row--)
60 p[row] = d[row] - mu[row] * p[row + 1] - gamma[row] * p[N - 1];
61 }
62
63 template<>
64 struct SolvePeriodicInterp1D<double>
65 {
66 static inline void apply(const std::vector<double>& data, std::vector<double>& p)
67 {
68 FuncSolvePeriodicInterp1D(data, p);
69 }
70
71 template<typename CT>
72 static inline void apply(const CT& data, CT& p, int N)
73 {
74 const double ratio = 0.25;
75 CT d(N), gamma(N), mu(N);
76 for (int i = 0; i < N; i++)
77 d[i] = 1.5 * data[i];
78 // First, eliminate leading coefficients
79 gamma[0] = ratio;
80 mu[0] = ratio;
81 mu[N - 1] = ratio;
82 gamma[N - 1] = 1.0;
83 for (int row = 1; row < (N - 1); row++)
84 {
85 double diag = 1.0 - mu[row - 1] * ratio;
86 double diagInv = 1.0 / diag;
87 gamma[row] = -ratio * gamma[row - 1] * diagInv;
88 mu[row] = diagInv * ratio;
89 d[row] = diagInv * (d[row] - ratio * d[row - 1]);
90 // Last row
91 d[N - 1] -= mu[N - 1] * d[row - 1];
92 gamma[N - 1] -= mu[N - 1] * gamma[row - 1];
93 mu[N - 1] = -mu[N - 1] * mu[row - 1];
94 }
95 // Last row: gamma(N-1) hold diagonal element
96 mu[N - 1] += ratio;
97 gamma[N - 1] -= mu[N - 1] * (mu[N - 2] + gamma[N - 2]);
98 d[N - 1] -= mu[N - 1] * d[N - 2];
99 p[N] = d[N - 1] / gamma[N - 1];
100 // Now go back upward, back substituting
101 for (int row = N - 2; row >= 0; row--)
102 p[row + 1] = d[row] - mu[row] * p[row + 2] - gamma[row] * p[N];
103 }
104 };
105
106 template<>
107 struct SolvePeriodicInterp1D<std::complex<double>>
108 {
109 typedef std::complex<double> value_type;
110 typedef double real_type;
111
112 static inline void apply(const std::vector<value_type>& data, std::vector<value_type>& p)
113 {
114 int N = data.size();
115 std::vector<real_type> dataReal(N), dataImag(N), pReal(N), pImag(N);
116 for (int i = 0; i < N; i++)
117 {
118 dataReal[i] = data[i].real();
119 dataImag[i] = data[i].imag();
120 }
121 SolvePeriodicInterp1D<double>::apply(dataReal, pReal);
122 SolvePeriodicInterp1D<double>::apply(dataImag, pImag);
123 for (int i = 0; i < N; i++)
124 p[i] = value_type(pReal[i], pImag[i]);
125 }
126 };
127
128
129 template<typename T>
130 struct SolveFirstDerivInterp1D
131 {};
132
133 template<>
134 struct SolveFirstDerivInterp1D<double>
135 {
136 template<class CT>
137 static inline void apply(const CT& data, CT& p, int N, double* bcLower, double* bcUpper)
138 {
139 const double ratio = 0.25;
140 CT d(N + 2), mu(N + 2, ratio);
141 for (int i = 1; i <= N; i++)
142 d[i] = 1.5 * data[i - 1];
143 double al = 0.25 * bcLower[0];
144 double bl = 0.25 * bcLower[1];
145 double cl = 0.25 * bcLower[2];
146 double dl = 1.5 * bcLower[3];
147 double ar = 0.25 * bcUpper[0];
148 double br = 0.25 * bcUpper[1];
149 double cr = 0.25 * bcUpper[2];
150 double dr = 1.5 * bcUpper[3];
151 // First, eliminate leading coefficients
152 double alInv = 1.0 / al;
153 bl *= alInv;
154 cl *= alInv;
155 dl *= alInv;
156 d[0] = dl;
157 mu[0] = bl;
158 mu[1] = ratio - ratio * cl;
159 for (int row = 1; row <= N; row++)
160 {
161 double diag = 1.0 - mu[row - 1] * ratio;
162 double diagInv = 1.0 / diag;
163 mu[row] *= diagInv;
164 d[row] = diagInv * (d[row] - ratio * d[row - 1]);
165 }
166 br -= ar * mu[N - 1];
167 dr -= ar * d[N - 1];
168 cr -= br * mu[N];
169 dr -= br * d[N];
170 p[N + 1] = dr / cr;
171 // Now go back upward, back substituting
172 for (int row = N; row >= 1; row--)
173 p[row] = d[row] - mu[row] * p[row + 1];
174 // And do 0th row
175 p[0] = dl - bl * p[1] - cl * p[2];
176 }
177 };
178
179 template<>
180 struct SolveFirstDerivInterp1D<std::complex<double>>
181 {
182 typedef std::complex<double> value_type;
183 typedef double real_type;
184
185 template<class CT>
186 static inline void apply(const CT& data, CT& p, int N, double* bcLower, double* bcUpper)
187 {
188 std::vector<real_type> dataReal(N), dataImag(N), pReal(N), pImag(N);
189 for (int i = 0; i < N; i++)
190 {
191 dataReal[i] = data[i].real();
192 dataImag[i] = data[i].imag();
193 }
194 SolveFirstDerivInterp1D<double>::apply(dataReal, pReal, N, bcLower, bcUpper);
195 SolveFirstDerivInterp1D<double>::apply(dataImag, pImag, N, bcLower, bcUpper);
196 for (int i = 0; i < N; i++)
197 p[i] = value_type(pReal[i], pImag[i]);
198 }
199 };
200
201 template<>
202 struct SolveFirstDerivInterp1D<float>
203 {
204 template<class CT>
205 static inline void apply(const CT& data, CT& p, int N, float* bcLower, float* bcUpper)
206 {
207 std::vector<double> data_d(N), p_d(N);
208 std::copy(data.begin(), data.end(), data_d.begin());
209 double bcLower_d[4];
210 double bcUpper_d[4];
211 std::copy(bcLower, bcLower + 4, bcLower_d);
212 std::copy(bcUpper, bcUpper + 4, bcUpper_d);
213
214 SolveFirstDerivInterp1D<double>::apply(data_d, p_d, N, bcLower_d, bcUpper_d);
215
216 std::copy(p_d.begin(), p_d.end(), p.begin());
217 }
218 };
219
220 template<>
221 struct SolveFirstDerivInterp1D<std::complex<float>>
222 {
223 typedef std::complex<float> value_type;
224 typedef float real_type;
225
226 template<class CT>
227 static inline void apply(const CT& data, CT& p, int N, float* bcLower, float* bcUpper)
228 {
229 std::vector<double> dataReal(N), dataImag(N), pReal(N), pImag(N);
230 for (int i = 0; i < N; i++)
231 {
232 dataReal[i] = data[i].real();
233 dataImag[i] = data[i].imag();
234 }
235
236 double bcLower_d[4];
237 double bcUpper_d[4];
238 std::copy(bcLower, bcLower + 4, bcLower_d);
239 std::copy(bcUpper, bcUpper + 4, bcUpper_d);
240 SolveFirstDerivInterp1D<double>::apply(dataReal, pReal, N, bcLower_d, bcUpper_d);
241 SolveFirstDerivInterp1D<double>::apply(dataImag, pImag, N, bcLower_d, bcUpper_d);
242
243 for (int i = 0; i < N; i++)
244 p[i] = value_type(static_cast<float>(pReal[i]), static_cast<float>(pImag[i]));
245 }
246 };
247 #endif
248