QMCPACK
Functions
SplineSolvers.h File Reference
+ Include dependency graph for SplineSolvers.h:
+ This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions

template<typename T >
void CubicSplineSolve (T const *x, T const *y, int n, T yp0, T ypn, T *y2)
 Solve for second derivatives for cubic splines. More...
 
template<typename Tg , typename T >
void QuinticSplineSolve (int N, const Tg *X, T *Y, T *B, T *C, T *D, T *E, T *F)
 template function: note that the range of data is [0,n) instead of [1,n] solve the linear system for the first and second derivatives Copy of routine from: "One dimensional spline interpolation algorithms" by Helmuth Spath QUINAT routine More...
 

Function Documentation

◆ CubicSplineSolve()

void CubicSplineSolve ( T const *  x,
T const *  y,
int  n,
T  yp0,
T  ypn,
T *  y2 
)

Solve for second derivatives for cubic splines.

Parameters
[in]xlocation of knots
[in]yvalues of function at knots
[in]nnumber of values (size of x,y, and y2)
[in]yp0boundary condition at start of interval - first derivative or natural BC (second derivative equals 0) if yp0 > 1e30)
[in]ypnboundary condition at end of interval - first derivative or natural BC (second derivative equals 0) if ypn > 1e30)
[out]y2second derivative at each knot

Definition at line 33 of file SplineSolvers.h.

References qmcplusplus::n.

Referenced by OneDimCubicSpline< T >::spline(), and qmcplusplus::TEST_CASE().

34 {
35  std::vector<T> u(n);
36  T bc_start1 = ((yp0 < 9.9000000000000002e+29) ? (2 * x[0] - 2 * x[1]) : (1));
37  T bc_start2 = ((yp0 < 9.9000000000000002e+29) ? (x[0] - x[1]) : (0));
38  T bc_end1 = ((ypn < 9.9000000000000002e+29) ? (x[n - 1] - x[n - 2]) : (0));
39  T bc_end2 = ((ypn < 9.9000000000000002e+29) ? (2 * x[n - 1] - 2 * x[n - 2]) : (1));
40  T rhs_start =
41  ((yp0 < 9.9000000000000002e+29) ? (6 * yp0 + 6 * y[0] / (-x[0] + x[1]) - 6 * y[1] / (-x[0] + x[1])) : (0));
42  T rhs_end = ((ypn < 9.9000000000000002e+29)
43  ? (6 * ypn - 6 * y[n - 1] / (x[n - 1] - x[n - 2]) + 6 * y[n - 2] / (x[n - 1] - x[n - 2]))
44  : (0));
45  y2[0] = bc_start2 / bc_start1;
46  u[0] = rhs_start / bc_start1;
47  for (auto i = 1; i < n - 1; i += 1)
48  {
49  T z0 = -x[i];
50  T z1 = z0 + x[i + 1];
51  T z2 = z0 + x[i - 1];
52  T z3 = 2 * x[i + 1];
53  T z4 = 2 * x[i - 1];
54  T z5 = z2 * y2[i - 1];
55  T z6 = -y[i];
56  u[i] = (z1 * z2 * z2 * u[i - 1] - 6 * z1 * (z6 + y[i - 1]) + 6 * z2 * (z6 + y[i + 1])) / (z1 * z2 * (z3 - z4 + z5));
57  y2[i] = (-x[i + 1] + x[i]) / (-z3 + z4 - z5);
58  };
59  y2[n - 1] = (-bc_end1 * u[n - 2] + rhs_end) / (-bc_end1 * y2[n - 2] + bc_end2);
60  for (auto i = n - 2; i >= 0; i += -1)
61  {
62  y2[i] = u[i] - y2[i + 1] * y2[i];
63  };
64 };

◆ QuinticSplineSolve()

void QuinticSplineSolve ( int  N,
const Tg *  X,
T *  Y,
T *  B,
T *  C,
T *  D,
T *  E,
T *  F 
)
inline

template function: note that the range of data is [0,n) instead of [1,n] solve the linear system for the first and second derivatives Copy of routine from: "One dimensional spline interpolation algorithms" by Helmuth Spath QUINAT routine

Definition at line 78 of file SplineSolvers.h.

References qmcplusplus::abs(), B(), qmcplusplus::Units::charge::C, and qmcplusplus::Units::force::N.

Referenced by OneDimQuinticSpline< Td, Tg, CTd, CTg >::spline().

79 {
80  int M;
81  // mmorales: had issues setting std::numeric_limits<T>::epsilon() with gcc, FIX later
82  T eps = 1.0e-9;
83  //eps= std::numeric_limits<double>::epsilon();
84  T Q, R, Q2, Q3, R2, QR, P, P2, PQ, PR, PQQR, B1, V, TT, S, U, P3;
85  M = N - 2;
86  Q = X[1] - X[0];
87  R = X[2] - X[1];
88  Q2 = Q * Q;
89  R2 = R * R;
90  QR = Q + R;
91  D[0] = E[0] = D[1] = 0.0;
92  if (std::abs(Q) > eps)
93  D[1] = 6.0 * Q * Q2 / (QR * QR);
94  for (int i = 1; i < M; i++)
95  {
96  P = Q;
97  Q = R;
98  R = X[i + 2] - X[i + 1];
99  P2 = Q2;
100  Q2 = R2;
101  R2 = R * R;
102  PQ = QR;
103  QR = Q + R;
104  if (std::abs(Q) <= eps)
105  {
106  D[i + 1] = E[i] = F[i - 1] = 0.0;
107  }
108  else
109  {
110  Q3 = Q2 * Q;
111  PR = P * R;
112  PQQR = PQ * QR;
113  D[i + 1] = 6.0 * Q3 / (QR * QR);
114  D[i] += 2.0 * Q *
115  (15.0 * PR * PR + (P + R) * Q * (20.0 * PR + 7.0 * Q2) + Q2 * (8.0 * (P2 + R2) + 21.0 * PR + 2.0 * Q2)) /
116  (PQQR * PQQR);
117  D[i - 1] += 6.0 * Q3 / (PQ * PQ);
118  E[i] = Q2 * (P * QR + 3.0 * PQ * (QR + 2.0 * R)) / (PQQR * QR);
119  E[i - 1] += Q2 * (R * PQ + 3.0 * QR * (PQ + 2.0 * P)) / (PQQR * PQ);
120  F[i - 1] = Q3 / PQQR;
121  }
122  }
123  if (std::abs(R) > eps)
124  D[N - 3] += 6.0 * R * R2 / (QR * QR);
125  for (int i = 1; i < N; i++)
126  {
127  if (std::abs(X[i] - X[i - 1]) < eps)
128  {
129  B[i] = Y[i];
130  Y[i] = Y[i - 1];
131  }
132  else
133  {
134  B[i] = (Y[i] - Y[i - 1]) / (X[i] - X[i - 1]);
135  }
136  }
137  for (int i = 2; i < N; i++)
138  {
139  if (std::abs(X[i] - X[i - 2]) < eps)
140  {
141  C[i] = B[i] * 0.5;
142  B[i] = B[i - 1];
143  }
144  else
145  {
146  C[i] = (B[i] - B[i - 1]) / (X[i] - X[i - 2]);
147  }
148  }
149  // assume N > 5
150  P = C[0] = E[N - 3] = F[0] = F[N - 4] = F[N - 3] = 0.0;
151  C[1] = C[3] - C[2];
152  D[1] = 1.0 / D[1];
153  for (int i = 2; i < M; i++)
154  {
155  Q = D[i - 1] * E[i - 1];
156  D[i] = 1.0 / (D[i] - P * F[i - 2] - Q * E[i - 1]);
157  E[i] -= Q * F[i - 1];
158  C[i] = C[i + 2] - C[i + 1] - P * C[i - 2] - Q * C[i - 1];
159  P = D[i - 1] * F[i - 1];
160  }
161  C[N - 2] = 0.0;
162  C[N - 1] = 0.0;
163  for (int i = N - 3; i > 0; i--)
164  C[i] = (C[i] - E[i] * C[i + 1] - F[i] * C[i + 2]) * D[i];
165  M = N - 1;
166  Q = X[1] - X[0];
167  R = X[2] - X[1];
168  B1 = B[1];
169  Q3 = Q * Q * Q;
170  QR = Q + R;
171  if (std::abs(QR) < eps)
172  {
173  V = TT = 0.0;
174  }
175  else
176  {
177  V = C[1] / QR;
178  TT = V;
179  }
180  F[0] = 0.0;
181  if (std::abs(Q) > eps)
182  F[0] = V / Q;
183  for (int i = 1; i < M; i++)
184  {
185  P = Q;
186  Q = R;
187  R = 0.0;
188  if (i != (M - 1))
189  R = X[i + 2] - X[i + 1];
190  P3 = Q3;
191  Q3 = Q * Q * Q;
192  PQ = QR;
193  QR = Q + R;
194  S = TT;
195  TT = 0;
196  if (std::abs(QR) > eps)
197  TT = (C[i + 1] - C[i]) / QR;
198  U = V;
199  V = TT - S;
200  if (std::abs(PQ) < eps)
201  {
202  C[i] = C[i - 1];
203  D[i] = E[i] = F[i] = 0.0;
204  }
205  else
206  {
207  F[i] = F[i - 1];
208  if (std::abs(Q) > eps)
209  F[i] = V / Q;
210  E[i] = 5.0 * S;
211  D[i] = 10.0 * (C[i] - Q * S);
212  C[i] = D[i] * (P - Q) + (B[i + 1] - B[i] + (U - E[i]) * P3 - (V + E[i]) * Q3) / PQ;
213  B[i] = (P * (B[i + 1] - V * Q3) + Q * (B[i] - U * P3)) / PQ - P * Q * (D[i] + E[i] * (Q - P));
214  }
215  }
216  P = X[1] - X[0];
217  S = F[0] * P * P * P;
218  E[0] = D[0] = E[N - 1] = D[N - 1] = 0.0;
219  C[0] = C[1] - 10.0 * S;
220  B[0] = B1 - (C[0] + S) * P;
221  Q = X[N - 1] - X[N - 2];
222  TT = F[N - 2] * Q * Q * Q;
223  C[N - 1] = C[N - 2] + 10.0 * TT;
224  B[N - 1] = B[N - 1] + (C[N - 1] - TT) * Q;
225 }
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
B
double B(double x, int k, int i, const std::vector< double > &t)
Definition: soecp_eval_reference.cpp:27