Skip to content

Commit fb51b66

Browse files
pasabanovpasabanov
committed
Modified "Clamped", "FixedSecond" and "FixedThird" cubic splines to use tridiag_solve 
1 parent 4e37fef commit fb51b66

File tree

1 file changed

+41
-81
lines changed

1 file changed

+41
-81
lines changed

‎ALFI/ALFI/spline/cubic.h‎

Lines changed: 41 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -201,40 +201,22 @@ namespace alfi::spline {
201201
},
202202
[&](const typename Types::Clamped& c) {
203203
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
204-
205-
Container<Number> C(n);
206-
if (n == 2) {
207-
C[0] = (3*dY[0]/dX[0] - c.df_n - 2*c.df_1) / dX[0];
208-
C[1] = (c.df_1 + 2*c.df_n - 3*dY[0]/dX[0]) / dX[0];
209-
} else if (n == 3) {
210-
C[1] = (3*(dY[1]/dX[1] - dY[0]/dX[0]) + (c.df_1 - c.df_n)) / (dX[0] + dX[1]);
211-
C[0] = (3*(dY[0]/dX[0] - c.df_1)/dX[0] - C[1]) / 2;
212-
C[2] = (3*(c.df_n - dY[1]/dX[1])/dX[1] - C[1]) / 2;
213-
} else {
214-
Container<Number> diag(n - 2), right(n - 2);
215-
diag[0] = 1.5*dX[0] + 2*dX[1];
216-
right[0] = 3*dY[1]/dX[1] - 4.5*dY[0]/dX[0] + 1.5*c.df_1;
217-
for (SizeT i = 1; i < n - 3; ++i) {
218-
diag[i] = 2 * (dX[i] + dX[i+1]);
219-
right[i] = 3 * (dY[i+1]/dX[i+1] - dY[i]/dX[i]);
220-
}
221-
diag[n-3] = 2*dX[n-3] + 1.5*dX[n-2];
222-
right[n-3] = 4.5*dY[n-2]/dX[n-2] - 3*dY[n-3]/dX[n-3] - 1.5*c.df_n;
223-
224-
for (SizeT i = 1; i < n - 2; ++i) {
225-
const auto m = dX[i] / diag[i-1];
226-
diag[i] -= m * dX[i];
227-
right[i] -= m * right[i-1];
228-
}
229-
230-
C[n-2] = right[n-3] / diag[n-3];
231-
for (SizeT i = n - 3; i >= 1; --i) {
232-
C[i] = (right[i-1] - dX[i] * C[i+1]) / diag[i-1];
233-
}
234-
C[n-1] = 1.5*(c.df_n/dX[n-2] - dY[n-2]/dX[n-2]/dX[n-2]) - 0.5*C[n-2];
235-
C[0] = 1.5*(dY[0]/dX[0]/dX[0] - c.df_1/dX[0]) - 0.5*C[1];
204+
Container<Number> lower(n), diag(n), upper(n), right(n);
205+
for (SizeT i = 1; i < n - 1; ++i) {
206+
lower[i] = dX[i-1];
207+
diag[i] = 2 * (dX[i-1] + dX[i]);
208+
upper[i] = dX[i];
209+
right[i] = 3 * (dY[i]/dX[i] - dY[i-1]/dX[i-1]);
236210
}
237-
211+
lower[0] = NAN;
212+
diag[0] = 2*dX[0];
213+
upper[0] = dX[0];
214+
right[0] = 3 * (dY[0]/dX[0] - c.df_1);
215+
lower[n-1] = dX[n-2];
216+
diag[n-1] = 2*dX[n-2];
217+
upper[n-1] = NAN;
218+
right[n-1] = 3 * (c.df_n - dY[n-2]/dX[n-2]);
219+
const auto C = util::linalg::tridiag_solve(lower, diag, upper, right);
238220
for (SizeT i = 0, index = 0; i < n - 1; ++i) {
239221
coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]);
240222
coeffs[index++] = C[i];
@@ -244,29 +226,22 @@ namespace alfi::spline {
244226
},
245227
[&](const typename Types::FixedSecond& s) {
246228
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
247-
248-
Container<Number> C(n);
249-
C[0] = s.d2f_1 / 2;
250-
C[n-1] = s.d2f_n / 2;
251-
if (n > 2) {
252-
Container<Number> diag(n - 2), right(n - 2);
253-
for (SizeT i = 0; i < n - 2; ++i) {
254-
diag[i] = 2 * (dX[i] + dX[i+1]);
255-
right[i] = 3 * (dY[i+1]/dX[i+1] - dY[i]/dX[i]);
256-
}
257-
right[0] -= C[0] * dX[0];
258-
right[n-3] -= C[n-1] * dX[n-2];
259-
for (SizeT i = 1; i < n - 2; ++i) {
260-
const auto m = dX[i] / diag[i-1];
261-
diag[i] -= m * dX[i];
262-
right[i] -= m * right[i-1];
263-
}
264-
C[n-2] = right[n-3] / diag[n-3];
265-
for (SizeT i = n - 3; i >= 1; --i) {
266-
C[i] = (right[i-1] - dX[i] * C[i+1]) / diag[i-1];
267-
}
229+
Container<Number> lower(n), diag(n), upper(n), right(n);
230+
for (SizeT i = 1; i < n - 1; ++i) {
231+
lower[i] = dX[i-1];
232+
diag[i] = 2 * (dX[i-1] + dX[i]);
233+
upper[i] = dX[i];
234+
right[i] = 3 * (dY[i]/dX[i] - dY[i-1]/dX[i-1]);
268235
}
269-
236+
lower[0] = NAN;
237+
diag[0] = 1;
238+
upper[0] = 0;
239+
right[0] = s.d2f_1 / 2;
240+
lower[n-1] = 0;
241+
diag[n-1] = 1;
242+
upper[n-1] = NAN;
243+
right[n-1] = s.d2f_n / 2;
244+
const auto C = util::linalg::tridiag_solve(lower, diag, upper, right);
270245
for (SizeT i = 0, index = 0; i < n - 1; ++i) {
271246
coeffs[index++] = (C[i+1] - C[i]) / (3*dX[i]);
272247
coeffs[index++] = C[i];
@@ -276,43 +251,28 @@ namespace alfi::spline {
276251
},
277252
[&](const typename Types::FixedThird& t) {
278253
const auto dX = util::arrays::diff(X), dY = util::arrays::diff(Y);
279-
280254
Container<Number> C(n);
281255
if (n == 2) {
282256
C[0] = -dX[0] * (t.d3f_1 + t.d3f_n) / 8;
283257
C[1] = -C[0];
284258
} else {
285-
Container<Number> diag(n), right(n);
286-
diag[0] = 1;
259+
Container<Number> lower(n), diag(n), upper(n), right(n);
287260
for (SizeT i = 1; i < n - 1; ++i) {
261+
lower[i] = dX[i-1];
288262
diag[i] = 2 * (dX[i-1] + dX[i]);
289-
}
290-
diag[n-1] = 1;
291-
right[0] = -dX[0] * t.d3f_1 / 2;
292-
for (SizeT i = 1; i < n - 1; ++i) {
263+
upper[i] = dX[i];
293264
right[i] = 3 * (dY[i]/dX[i] - dY[i-1]/dX[i-1]);
294265
}
266+
lower[0] = NAN;
267+
diag[0] = 1;
268+
upper[0] = -1;
269+
right[0] = -dX[0] * t.d3f_1 / 2;
270+
lower[n-1] = -1;
271+
diag[n-1] = 1;
272+
upper[n-1] = NAN;
295273
right[n-1] = dX[n-2] * t.d3f_n / 2;
296-
297-
const auto m_1 = dX[0] / diag[0];
298-
diag[1] -= m_1 * (-1);
299-
right[1] -= m_1 * right[0];
300-
for (SizeT i = 2; i < n - 1; ++i) {
301-
const auto m = dX[i-1] / diag[i-1];
302-
diag[i] -= m * dX[i-1];
303-
right[i] -= m * right[i-1];
304-
}
305-
const auto m_n_1 = (-1) / diag[n-2];
306-
diag[n-1] -= m_n_1 * dX[n-2];
307-
right[n-1] -= m_n_1 * right[n-2];
308-
309-
C[n-1] = right[n-1] / diag[n-1];
310-
for (SizeT i = n - 2; i >= 1; --i) {
311-
C[i] = (right[i] - dX[i] * C[i+1]) / diag[i];
312-
}
313-
C[0] = right[0] + C[1];
274+
C = util::linalg::tridiag_solve<Number,Container>(lower, diag, upper, right);
314275
}
315-
316276
Container<Number> D(n - 1);
317277
for (SizeT i = 1; i < D.size() - 1; ++i) {
318278
D[i] = (C[i+1] - C[i]) / (3*dX[i]);

0 commit comments

Comments
 (0)